***********************************************************************
*     lib_clean62 - library for CLEAN subroutines v6.2                *
***********************************************************************
***********************************************************************
*     msrc - dirty beam of gaussian sources                           *
*       nantst : input,integer*(84)  antenna status                   *
*       taper  : input,real taper at 160d                             *
*       gbeam  : output,real(512,512,0:10) dirty beam of              *
*                  gaussian sources                                   *
*                  width : 0(point source),1,...,10                   *
*                                         July 13,1992  Y. Hanaoka    *
*                                         Mar  12,1996     koshix     *
*                                     November 27,1998  Y. Hanaoka    *
***********************************************************************
c
      subroutine msrc(nantst,taper,gbeam)
c
      integer nantst(1)
      real gbeam(512,512,0:10),a(512,512),b(512,512)
c
      do 100 iws=0,10
c
***** clean source *****
c
      do 200 il=1,512
      do 200 i=1,512
       a(i,il)=0.
       a(i,il)=0.
  200 continue
c
      if (iws.eq.0) then
       a(257,257)=100.
      else
       wsr=real(iws)
       sum=0.
       do 410 l=1,512
       do 410 i=1,512
        a(i,l)=a(i,l)+100.*exp(-((257.-real(i))*(257.-real(i))
     -         +(257.-real(l))*(257.-real(l)))/wsr/wsr)
        sum=sum+a(i,l)
  410  continue
       do 420 l=1,512
       do 420 i=1,512
        a(i,l)=a(i,l)*100./sum
  420  continue
      end if
c
c      write(6,'('' '')')
c      do 400 il=254,260
c  400 write(6,'('' '',7f9.2)') (a(i,il),i=254,260)
c
***** dirty map *****
c
      call cft512(a,b,1)
      call uvobs(nantst,taper,a,b)
      call cft512(a,b,-1)
c
      do 300 j=1,512
      do 300 i=1,512
  300 gbeam(i,j,iws)=a(i,j)
c
  100 continue
c
      return
      end
***********************************************************************
*     clhis - calculate background and disk brightness                *
*               of a dirty map                                        *
*       map    : input,real(512,512)  dirty map                       *
*       mstype : input,integer        type of model sun               *
*       itvbin : input,integer        histogram binsize               *
*       bgr    : output,real background(sky) brightness               *
*       dbr    : output,real disk brightness                          *
*                                     November 22,1995  Y. Hanaoka    *
***********************************************************************
c
      subroutine clhis(map,mstype,itvbin,bgr,dbr)
c
      real map(512,512)
      real his0(2000),his1(2000)
      real pargau(6),chis(2000),vw(2006,8)
c
      external hiserr
      common /hist/his0
c
      do 150 i=1,2000
  150 his0(i)=0.
c
***** histogram *****
      do 200 j=1,512
      do 200 i=1,512
       ibin=(int(map(i,j))+itvbin*1000)/itvbin
       if (ibin.ge.1 .and. ibin.le.2000) his0(ibin)=his0(ibin)+1.
  200 continue
c
***** moving mean of histogram *****
      do 250 j=1,1976
      do 250 i=1,24
  250 his0(j)=his0(j)+his0(j+i)
c
***** primary gaussian (initial) *****
      call hisgau(his0,pargau(1))
c
*****
      if (mstype .eq. 2) then
       do 400 i=int(pargau(2)+pargau(3)*2.),2000
  400  his0(i)=0.
      endif
      if (mstype .eq. 3) then
       do 410 i=1,int(pargau(2)-pargau(3)*2.)
  410  his0(i)=0.
      endif
c
***** subtraction of the primary gaussian *****
      do 300 i=1,1976
       his1(i)=his0(i)-pargau(1)*exp(-(real(i)-pargau(2))
     -  *(real(i)-pargau(2))/(pargau(3)*pargau(3)/0.693147))
c       his1(i)=his0(i)-0.6*pargau(1)*exp(-(real(i)-pargau(2))
c     -  *(real(i)-pargau(2))/(0.9*pargau(3)*0.9*pargau(3)))
c     -  -0.4*pargau(1)*exp(-(real(i)-pargau(2))
c     -  *(real(i)-pargau(2))/(0.9*2.5*pargau(3)*0.9*2.5*pargau(3)))
  300 continue
c
***** mainly for mstype=1 *****
c
* centroid of histogram
      cn=0.
      cs=0.
      do 450 i=1,1976
c       if (his0(i).ge.100) then
c        cn=cn+real(his0(i))
c        cs=cs+real(i)*real(his0(i))
       if (his0(i).ge.1000) then
        cn=cn+sqrt(real(his0(i)))
        cs=cs+real(i)*sqrt(real(his0(i)))
       endif
  450 continue
      cs=cs/cn
ccc   write(6,'('' centroid ='',f10.1)') cs
c
*
      if (pargau(2).lt.int(cs)) then
       do 460 i=1,int(cs)
  460  his1(i)=his1(int(cs))
      else
       do 470 i=int(cs),1976
  470  his1(i)=his1(int(cs))
      endif
c
***** secondary gaussian (initial) *****
      call hisgau(his1,pargau(4))
c
ccc   write(6,'('' initial gaussian parameters ='',f10.1,f8.1,f6.1,
ccc  -    f10.1,f8.1,f6.1)') (pargau(i),i=1,6)
c
*****
c
      nmax=6*2000*100
      epsr=1.e-4
      call nolf1(pargau,6,hiserr,2000,epsr,nmax,chis,sums,
     -           vw,2006,icon)
ccc   if (icon.ne.0) write(6,'(''  clhis : icon of nolf1 ='',i6)')
ccc  -     icon
c      write(6,'(''  gaussfit : sums ='',e15.4,''  n ='',i8)')
c     -     sums,nmax
ccc   write(6,'('' fitted gaussian parameters  ='',f10.1,f8.1,f6.1,
ccc  -    f10.1,f8.1,f6.1)') (pargau(i),i=1,6)
c
      if (pargau(2).lt.pargau(5)) then
        db1=(pargau(2)-1000.)*real(itvbin)+12.5*real(itvbin)
        db2=(pargau(5)-1000.)*real(itvbin)+12.5*real(itvbin)
      else
        db1=(pargau(5)-1000.)*real(itvbin)+12.5*real(itvbin)
        db2=(pargau(2)-1000.)*real(itvbin)+12.5*real(itvbin)
      endif
c      write(6,'('' db1,db2 ='',2f10.2)') db1,db2
      if (mstype.eq.1 .or. mstype.eq.2) then
       bgr=db1
       dbr=db2
      else
       bgr=db1-(db2-db1)
       dbr=db1
      endif
c
      return
      end
c
***********************************************************************
***********************************************************************
      subroutine hisgau(his,pargau)
c
      real his(2000),pargau(3)
c
***** peak of histogram *****
      hismax=0
      imax=1
      do 300 i=1,1976
       if (hismax.lt.his(i)) then
        hismax=his(i)
        imax=i
       end if
  300 continue
c
***** half of the peak *****
      ihalf=1
      do 320 i=1,1976
       if ((hismax/2. .gt. his(imax+i)) .or.
     -     (hismax/2. .gt. his(imax-i))) then
        ihalf=i
        goto 330
       end if
  320 continue
  330 continue
c
      pargau(1)=hismax
      pargau(2)=real(imax)
      pargau(3)=real(ihalf)
      return
      end
c
***********************************************************************
***********************************************************************
      subroutine hiserr(pargau,hise)
c
      real pargau(6),hise(2000),his0(2000)
c
      common /hist/his0
c
      do 100 i=1,2000
       if (his0(i) .lt. 2000.) then
        hise(i)=0.
       else
        hise(i)=his0(i)-pargau(1)*exp(-(real(i)-pargau(2))
     -     *(real(i)-pargau(2))/(pargau(3)*pargau(3)/0.693))
     -     -pargau(4)*exp(-(real(i)-pargau(5))
     -     *(real(i)-pargau(5))/(pargau(6)*pargau(6)/0.693))
c        hise(i)=his0(i)-0.6*pargau(1)*exp(-(real(i)-pargau(2))
c     -     *(real(i)-pargau(2))/(0.9*pargau(3)*0.9*pargau(3)))
c     -     -0.4*pargau(1)*exp(-(real(i)-pargau(2))
c     -     *(real(i)-pargau(2))/(0.9*2.5*pargau(3)*0.9*2.5*pargau(3)))
c     -     -0.6*pargau(4)*exp(-(real(i)-pargau(5))
c     -     *(real(i)-pargau(5))/(0.9*pargau(6)*0.9*pargau(6)))
c     -     -0.4*pargau(4)*exp(-(real(i)-pargau(5))
c     -     *(real(i)-pargau(5))/(0.9*2.5*pargau(6)*0.9*2.5*pargau(6)))
        wt=exp(-(real(i)-pargau(2))*(real(i)-pargau(2))
     -    /(pargau(3)*pargau(3)*0.2))
     -    +exp(-(real(i)-pargau(5))*(real(i)-pargau(5))
     -    /(pargau(6)*pargau(6)*0.2))
        hise(i)=hise(i)*wt
       endif
  100 continue
c      write(6,'('' '',10f7.1)') (hise(i),i=801,1200,40)
c
      return
      end
*************************************************************************
*     msun17 - model sun for 17                                         *
*       nfreq  : input,integer observation frequency                    *
*       pmat   : input,real(4) projection matrix                        *
*       solr   : input,real    radius of the model sun (arcsec)         *
*       nantst : input,integer*(84)  antenna status                     *
*       pxew   : input real x-offset of the center                      *
*       pxns   : input,real y-offset of the center                      *
*       taper  : input,real taper at 160d                               *
*       mapd   : output,real(512,512) dirty map                         *
*       mapc   : output,real(512,512) clean map                         *
*       mstype : output,integer  type of model sun                      *
*                                      November 22,1995  Y. Hanaoka     *
*                                      Jun      13,1995     koshix      *
*                                      November 27,1998  Y. Hanaoka     *
*************************************************************************
c
      subroutine msun17(pmat,solr,
     -     nantst,pxew,pxns,taper,mapd,mapc,mstype)
c
      integer nantst(1),mshis(4)
      real mapd(512,512),mapc(512,512),mapa(512,512),pmat(4)
c
      rd=57.29577951
      pi=3.141592654
c
***** model sun at 17GHz *****
c
       rsun=solr/4.64947
c
       do 400 l=1,512
       do 400 i=1,512
        xi=amod(real(i)-pxew+512.,512.)-257.
        xl=amod(real(l)-pxns+512.,512.)-257.
        xc=(pmat(4)*xi-pmat(3)*xl/1.000112494)
     -     /(pmat(1)*pmat(4)-pmat(3)*pmat(2))
        yc=(-pmat(2)*xi+pmat(1)*xl/1.000112494)
     -     /(pmat(1)*pmat(4)-pmat(3)*pmat(2))
        if (xc*xc+yc*yc .gt. rsun*rsun) then
         mapc(i,l)=0.
        else
         mapc(i,l)=100.
        end if
  400  continue
c
***** type of brightness distribution of model sun *****
* histogram of model sun *
c
      do 540 i=1,4
  540 mshis(i)=0
      do 550 l=1,512
      do 550 i=1,512
       j=int(mapc(i,l)+1.)/100+1
       mshis(j)=mshis(j)+1
  550 continue
c
      if (mshis(1) .ge. mshis(2)/3.) mstype=1
      if ((mshis(1) .lt. mshis(2)/3.) .and. (mshis(1).gt.mshis(3)))
     -  mstype=2
      if ((mshis(3) .lt. mshis(2)/3.) .and. (mshis(3).gt.mshis(1)))
     -  mstype=3
      if (mshis(3) .ge. mshis(2)/3.) mstype=4
ccc   write(6,'('' mstype ='',i4)') mstype
c
       call cft512(mapc,mapa,1)
       call filter1(taper,mapc)
       call filter1(taper,mapa)
       call cft512(mapc,mapa,-1)
c
***** dirty disk *****
c
      do 650 l=1,512
      do 650 i=1,512
  650 mapd(i,l)=mapc(i,l)
c
      call cft512(mapd,mapa,1)
      call uvobs(nantst,1.,mapd,mapa)
      call cft512(mapd,mapa,-1)
c
      return
      end
c
***********************************************************************
*     clean17 - clean process for 17 GHz images                       *
*         model source :  gaussian source and flat disk               *
*         adaptive clean-criterion definition available               *
*       menu   : input,integer                                        *
*                  2 : disk/pointsource(+) clean with restoration     *
*                 -2 : disk/pointsource(+/-) clean with restoration   *
*                 22 : disk/pointsource(+) clean with point source    *
*                        restoration, without disk restoration        *
*                -22 : disk/pointsource(+/-) clean with point source  *
*                        restoration, without disk restoration        *
*                 -1 : pointsource(+/-) clean with restore            *
*                  1 : pointsource(+) clean with restore              *
*                  3 : preclean only, criter=dskbr*5.+bgr             *
*       map    : input,real(512,512)  dirty map                       *
*                output,real(512,512) resudual                        *
*       pmat   : input,real(4) projection matrix                      *
*       solr   : input,real    radius of the model sun (arcsec)       *
*       nantst : input,integer*(84)  antenna status                   *
*       gbeam  : input,real(512,512,0:10) dirty beam of gaussian      *
*                   sources                                           *
*       criter : input,real  (lower limit of) clean criterion         *
*                 if criter < 0, -criter*dskbr is actually adopted    *
*       crfac  : input real  factor for clean criterion               *
*                 actual clean criterion is the larger one of         *
*                 criter or crfac*(maximum of the map)                *
*       cr     : output,real actual clean criterion                   *
*       pxew   : input/output,real x-offset of disk center (pixel)    *
*       pxns   : input/output,real y-offset of disk center (pixel)    *
*                  only valid for disk clean menu                     *
*                  else when pxew=999. or pxns=999., x-offset and     *
*                  y-offset is calculated and output pxew and pxns    *
*                  is set to calculated offset.                       *
*       nmax   : input,integer  maximum number of clean components    *
*              : output,integer number of clean components            *
*       bgr    : output,real   background brightness                  *
*       dskbr  : input,real    disk brightness (valid for criter < 0  *
*                              and menu=+/-1,+/-11                    *
*              : output,real   disk brightness                        *
*       pi     : output,integer(nmax) x-coordinate of clean points    *
*       pl     : output,integer(nmax) y-coordinate of clean points    *
*       pw     : output,integer(nmax) width of clean points           *
*       pa     : output,real(nmax)    intensity of clean points       *
*       cmap   : output,real(512,512) clean map                       *
*                                        June 12, 1998  Y. Hanaoka    *
***********************************************************************
c
      subroutine clean17(menu,map,pmat,solr,
     -     nantst,gbeam,criter,crfac,cr,
     -     pxew,pxns,corf,nmax,bgr,dskbr,pi,pl,pw,pa,cmap)
c
      integer nantst(84),pi(nmax),pl(nmax),pw(nmax)
      real map(512,512),cmap(512,512),ddisk(512,512),cdisk(512,512),
     -     cbeam(-30:30,-30:30),gbeam(512,512,0:10),
     -     pmat(4),pa(nmax)
c
      n=0
c
      call mean(1,map,am,amxmap,imx,lmx)
ccc   write(6,'('' maximum of the map :'',f12.1)') amxmap
c
      if (mod(iabs(menu),10) .ge. 2) then
c
*********************************************
*     preclean process for disk clean       *
*********************************************
c
***** model sun *****
       call msun17(pmat,solr,
     -      nantst,0.,0.,1.,ddisk,cdisk,mstype)
c
       call clhis(map,mstype,300,bgr,dbr)
       dskbr=dbr-bgr
ccc    write(6,'('' background/disk(initial) :'',2f10.2)')
ccc  -       bgr,dskbr
       if (criter.lt.0.) then
        cr=amax1(amxmap*crfac,dskbr*abs(criter))
       else
        cr=amax1(amxmap*crfac,criter)
       endif
       cr=amax1(dskbr*5.+bgr,cr*2.)
ccc    write(6,'('' preclean criterion :'',f10.1)') cr
c
       do 100 while(.true.)
c
***** search maximum *****
        call mean(menu,map,am,amx,imx,lmx)
        if (amx.gt.cr .and. n.lt.nmax) then
c
***** subtract point source component *****
         n=n+1
         pi(n)=imx
         pl(n)=lmx
         call cleanp(map,gbeam,0.02,amx,imx,lmx,pw(n),pa(n))
c        write(6,'('' '',i6,'': am,amx,w,a,imx,lmx ='',
c     -       f10.2,f14.2,i3,f10.2,2i5)')
c     -       n,am,amx,pw(n),pa(n),pi(n),pl(n)
c
        else
         go to 110
        end if
c
  100  continue
  110  if (n.ne.0) then
ccc     write(6,'('' n ='',i6,'': am,amx,w,a,imx,lmx ='',
ccc  -       f10.2,f14.2,i3,f10.2,2i5)')
ccc  -       n,am,amx,pw(n),pa(n),pi(n),pl(n)
       else
ccc     write(6,'('' n ='',i6)') n
       end if
c
*********************************************
*     disk clean process                    *
*********************************************
c
       call clhis(map,mstype,300,bgr,dbr)
       dskbr=dbr-bgr
ccc    write(6,'('' background/disk          :'',2f10.2)')
ccc  -       bgr,dskbr
c
c       if (mod(iabs(menu),10).eq.3) return
c
***** subtract disk component *****
c
       do 140 l=1,512
       do 140 i=1,512
  140  map(i,l)=map(i,l)-bgr
c
***** calculate disk center position ***** 
       if (pxew.eq.999. .or. pxns.eq.999.) then
        call possh17(map,dskbr,ddisk,pxew,pxns,corf)
ccc      write(6,'(''  pxew,pxns,corf ='',3f8.3)') pxew,pxns,corf
       end if
c
       if (mod(iabs(menu),10).eq.3) return
c
***** model sun with offset *****
       call msun17(pmat,solr,
     -      nantst,pxew,pxns,1.,ddisk,cdisk,mstype)
c
***** subtract disk component *****
       do 150 l=1,512
       do 150 i=1,512
        map(i,l)=map(i,l)-dskbr*ddisk(i,l)/100.
 150   continue
c
      end if
c
c
*********************************************
*     clean process                         *
*********************************************
c
***** clean criterion *****
c
       if (criter.lt.0.) then
        cr=amax1(amxmap*crfac,dskbr*abs(criter))
       else
        cr=amax1(amxmap*crfac,criter)
       endif
ccc    write(6,'('' clean criterion    :'',f10.1)') cr
c
       do 160 while(.true.)
c
***** search maximum *****
        call mean(menu,map,am,amx,imx,lmx)
        if (abs(amx).gt.cr .and. n.lt.nmax) then
c
***** subtract point source component *****
         n=n+1
         pi(n)=imx
         pl(n)=lmx
         call cleanp(map,gbeam,0.02,amx,imx,lmx,pw(n),pa(n))
c        write(6,'('' '',i6,'': am,amx,w,a,imx,lmx ='',
c     -       f10.2,f14.2,i3,f10.2,2i5)')
c     -       n,am,amx,pw(n),pa(n),pi(n),pl(n)
c
        else
         go to 170
        end if
c
  160  continue
  170  if (n.ne.0) then
ccc     write(6,'('' n ='',i6,'': am,amx,w,a,imx,lmx ='',
ccc  -       f10.2,f14.2,i3,f10.2,2i5)')
ccc  -       n,am,amx,pw(n),pa(n),pi(n),pl(n)
       else
ccc     write(6,'('' n ='',i6)') n
       end if
c
      nmax=n
      if (abs(menu)/10.eq.1) return
c
c
*********************************************
*     restore process                       *
*********************************************
c
      do 240 l=1,512
      do 240 i=1,512
  240 ddisk(i,l)=0.
c
      do 500 iw=0,10
c
***** create clean beam *****
c
       if (iw.eq.0) then
        do 210 l=-30,30
        do 210 i=-30,30
  210   cbeam(i,l)=0.
        cbeam(0,0)=100.
       else
        sum=0.
        do 220 l=-30,30
        do 220 i=-30,30
         cbeam(i,l)=100.*exp(-(real(i)*real(i)+real(l)*real(l))
     -     /real(iw)/real(iw))
         sum=sum+cbeam(i,l)
  220   continue
        do 230 l=-30,30
        do 230 i=-30,30
  230   cbeam(i,l)=100./sum*cbeam(i,l)
       end if   
c
***** restore point sources *****
c
       do 250 ip=1,n
        if (pw(ip).eq.iw) then
         do 260 l=-30,30
         do 260 i=-30,30
c  260    ddisk(pi(ip)+i,pl(ip)+l)=ddisk(pi(ip)+i,pl(ip)+l)+pa(ip)
c     -       *cbeam(i,l)/gbeam(257,257,iw)
          ii=mod(pi(ip)+i+511,512)+1
          ll=mod(pl(ip)+l+511,512)+1
          ddisk(ii,ll)=ddisk(ii,ll)
     -      +pa(ip)*cbeam(i,l)/gbeam(257,257,iw)
  260    continue
c        write(6,'('' '',i4,5f10.3)') iw,(cbeam(ii,0),ii=-10,10,5)
        end if
  250  continue
c
  500 continue
c
***** restore disk component *****
c
      cfact=0.
      if (iabs(menu).eq.2) cfact=dskbr
      do 300 l=1,512
      do 300 i=1,512
       cdisk(i,l)=ddisk(i,l)+cfact*cdisk(i,l)/100.
  300 continue
c
***** convolution by gaussian beam *****
c
       do 340 ll=-10,10
       do 340 ii=-10,10
  340  cbeam(ii,ll)=0.07461*exp(-0.234403
     -        *(real(ll)*real(ll)+real(ii)*real(ii)))
c
      do 310 l=1,512
      do 310 i=1,512
  310 cmap(i,l)=0.
      do 320 l=1,512
      do 320 i=1,512
       do 330 ll=-10,10
       do 330 ii=-10,10
        ib=mod(i+ii+511,512)+1
        lb=mod(l+ll+511,512)+1
        cmap(i,l)=cmap(i,l)+cdisk(ib,lb)*cbeam(ii,ll)
  330  continue
  320 continue
c
***** add residual *****
      do 350 l=1,512
      do 350 i=1,512
       cmap(i,l)=cmap(i,l)+map(i,l)
  350 continue
c
      return
      end
c
***********************************************************************
*      mean - calculate mean and search maximum                       *
*        menu : input,integer  clean menu                             *
*        map  : input,real(512,512)  dirty map                        *
*        am   : output,real    mean of map                            *
*        amx  : output,real    maximum of map                         *
*        imx  : output,integer x-coordinate of the maximum point      *
*        lmx  : output,integer y-coordinate of the maximum point      *
*                                         July 18,1992  Y. Hanaoka    *
***********************************************************************
c
      subroutine mean(menu,map,am,amx,imx,lmx)
c
      real map(512,512)
c
      am=0.
      amx=0.
      if (menu.gt.0) then
       do 100 l=1,512
       do 100 i=1,512
        am=am+map(i,l)
        if (amx.lt.map(i,l)) then
         amx=map(i,l)
         imx=i
         lmx=l
        end if
  100  continue
      else
       do 110 l=1,512
       do 110 i=1,512
        am=am+map(i,l)
        if (abs(amx).lt.abs(map(i,l))) then
         amx=map(i,l)
         imx=i
         lmx=l
        end if
  110  continue
      end if
      am=am/512./512.
c
      return
      end
c
***********************************************************************
*     cleanp - subtract a dirty beam component                        *
*       map    : input,real(512,512)  dirty map                       *
*              : output,real(512,512) dirty map, point source         *
*                  subtracted                                         *
*       gbeam  : input,real(512,512,0:10) dirty beam of gaussian      *
*                  sources                                            *
*       gainp  : input,real  CLEAN loop gain                          *
*       amx    : output,real    maximum of map                        *
*       imx    : output,integer x-coordinate of the maximum point     *
*       lmx    : output,integer y-coordinate of the maximum point     *
*       iwid   : output,integer width of subtracted point source      *
*       ap     : output,real    intensity of subtracted point source  *
*                                         July 17,1992  Y. Hanaoka    *
***********************************************************************
c
      subroutine cleanp(map,gbeam,gainp,amx,imx,lmx,iwid,ap)
c
      real map(512,512),gbeam(512,512,0:10)
c
      do 200 iw=1,10
       do 210 j=-5,5
       do 210 i=-5,5
  210  if (abs(map(imx+i,lmx+j)) .lt.
     -      gbeam(257+i,257+j,iw)/gbeam(257,257,iw)
     -      *abs(map(imx,lmx))) goto 290
  200 continue
  290 iwid=iw-1
c
      ap=amx*gainp
c
      do 300 l=1,512
      do 300 i=1,512
       ib=mod(i-imx+768,512)+1
       lb=mod(l-lmx+768,512)+1
  300 map(i,l)=map(i,l)-ap*gbeam(ib,lb,iwid)/gbeam(257,257,iwid)
c      do 400 il=1,512,64
c  400 write(6,'('' '',8f9.2)') (map(i,il),i=1,512,64)
c
      return
      end
***********************************************************************
*     possh17 - position shift between dirty model disk               *
*              and disk on dirty map                                  *
*       map    : input,real(512,512)  dirty map                       *
*       dskbr  : input,real    brightness of the dirty disk           *
*       disk   : input,real(512,512)  dirty map of model disk         *
*       pxew   : output,real x-offset of disk center                  *
*       pxns   : output,real y-offset of disk center                  *
*       corf   : output,real correlation between dirty map and model  *
*                                        June 12 ,1998  Y. Hanaoka    *
***********************************************************************
c
      subroutine possh17(map,dskbr,disk,pxew,pxns,corf)
c
      real map(512,512),disk(512,512),sumr(512),sumc(512),
     -     corfma(21,21),xp(25),yp(25),zc(25)
c
      pi=3.14159265
c
***** calculate centroid of dirty map *****
c
      sa=0.
      saa=0.
      do 200 i=1,512
       sumr(i)=0.
       sumc(i)=0.
  200 continue
      do 220 l=1,512
      do 220 i=1,512
       a=amin1(dskbr*1.5,amax1(0.,map(i,l)-dskbr/2.)
     -     +dskbr/4.+sign(dskbr/4.,map(i,l)-dskbr/2.))
       sa=sa+a
       saa=saa+a*a
       sumr(l)=sumr(l)+a
       sumc(i)=sumc(i)+a
  220 continue
c
      cs=0.
      sn=0.
      do 250 i=1,512
       cs=cs+sumr(i)*cos(real(i-257)/512.*2.*pi)
       sn=sn+sumr(i)*sin(real(i-257)/512.*2.*pi)
  250 continue
c      write(6,'('' cs/sn ='',2e15.5)') cs,sn
      if (abs(cs).lt.1.e-4) then
       ph=sign(pi/2,sn)
      else
c       ph=acos(cs/sqrt(cs*cs+sn*sn))
c       if (sn.lt.0.) ph=-ph
       ph=atan(sn/cs)
       if (cs.lt.0.) ph=ph+pi
      end if
      sl=ph/pi*256.
      isl=int(sl+0.5)
c
      cs=0.
      sn=0.
      do 260 i=1,512
       cs=cs+sumc(i)*cos(real(i-257)/512.*2.*pi)
       sn=sn+sumc(i)*sin(real(i-257)/512.*2.*pi)
  260 continue
c      write(6,'('' cs/sn ='',2e15.5)') cs,sn
      if (abs(cs).lt.1.e-4) then
       ph=sign(pi/2,sn)
      else
c       ph=acos(cs/sqrt(cs*cs+sn*sn))
c       if (sn.lt.0.) ph=-ph
       ph=atan(sn/cs)
       if (cs.lt.0.) ph=ph+pi
      end if
      sx=ph/pi*256.
      isx=int(sx+0.5)
c
ccc   write(6,'('' disk centroid ='',2f10.3)') sx,sl
c
***** calculate correlation and search maximum correlation *****
c
      sb=0.
      sbb=0.
      do 280 l=1,512
      do 280 i=1,512
       sb=sb+disk(i,l)
       sbb=sbb+disk(i,l)*disk(i,l)
  280 continue
c
***** pass-1 *****
c
      corfmx=0.
      do 300 loff=isl-10,isl+10
      do 300 ioff=isx-10,isx+10
       sab=0.
       do 310 l=1,512
       do 310 i=1,512
        ib=mod(i-ioff+511,512)+1
        lb=mod(l-loff+511,512)+1
        a=amin1(dskbr*1.5,amax1(0.,map(i,l)-dskbr/2.)
     -     +dskbr/4.+sign(dskbr/4.,map(i,l)-dskbr/2.))
        sab=sab+a*disk(ib,lb)
  310  continue
       corf=(512.*512.*sab-sa*sb)/sqrt((512.*512.*saa-sa*sa)
     -      *(512.*512.*sbb-sb*sb))
c       write(6,'('' corf ='',f8.3)') corf
       if (corf.gt.corfmx) then
        corfmx=corf
        ioffmx=ioff
        loffmx=loff
       end if
       corfma(ioff-isx+10+1,loff-isl+10+1)=corf
  300 continue
c
c      write(6,'('' corfm ='')')
c      write(6,'('' '',21f8.5)') corfma
c
      isx=ioffmx
      isl=loffmx
c      write(6,'(''  pass-1: isx  ='',i6,''   isl  ='',i6,
c     -''   correlation ='',f8.4)') isx,isl,corfmx
c
***** pass-2 *****
c
      do 400 loff=isl-2,isl+2
      do 400 ioff=isx-2,isx+2
       sab=0.
       do 410 l=1,512
       do 410 i=1,512
        ib=mod(i-ioff+511,512)+1
        lb=mod(l-loff+511,512)+1
        a=amin1(dskbr*1.5,amax1(0.,map(i,l)-dskbr/2.)
     -     +dskbr/4.+sign(dskbr/4.,map(i,l)-dskbr/2.))
        sab=sab+a*disk(ib,lb)
  410  continue
       corf=(512.*512.*sab-sa*sb)/sqrt((512.*512.*saa-sa*sa)
     -      *(512.*512.*sbb-sb*sb))
c       write(6,'('' corf ='',f8.3)') corf
       corfma(ioff-isx+2+1,loff-isl+2+1)=corf
  400 continue
c
c      write(6,'('' corfm ='')')
c      write(6,'('' '',11f8.5)') corfma
c
      do 450 i=1,25
  450 xp(i)=real(mod(i-1,5)-2)
c      write(6,'('' '',5f10.5)') xp
      do 460 i=1,25
  460 yp(i)=real(((i-1)/5)-2)
c      write(6,'('' '',5f10.5)') yp
      do 470 i=1,25
  470 zc(i)=corfma(mod(i-1,5)+1,(i-1)/5+1)
c      write(6,'('' '',5f10.5)') zc
c
      call getctr(xp,yp,zc,25,pxew,pxns,corf)
c      write(6,'(''  pass-2: pxew ='',f8.3,''   pxns ='',f8.3)')
c     -      pxew,pxns
      pxew=pxew+real(isx)
      pxns=pxns+real(isl)
ccc   write(6,'('' disk center position: pxew ='',f8.3,''   pxns ='',
ccc  -     f8.3,''   correlation ='',f8.4)') pxew,pxns,corf
c
      return
      end
c
c
***********************************************************************
*     getctr - get the summit point of two-dimensional surface        *
*                by least square method                               *
*       x      : input,real(n)  x-positions                           *
*       y      : input,real(n)  y-positions                           *
*       z      : input,real(n)  data value at points [x(i),y(i)]      *
*       n      : input,integer  number of data                        *
*       pxew   : output,real x-coordinate of the summit point         *
*       pxns   : output,real y-coordinate of the summit point         *
*       corf   : output,real correlation at the summit point          *
*                                     November 30,1995  Y. Hanaoka    *
***********************************************************************
c
      subroutine getctr(x,y,z,n,pxew,pxns,corf)
c
      real a(6,6),b(6),vw(6),ip(6)
      real x(n),y(n),z(n)
c
      sx4=0.
      sx3y=0.
      sx2y2=0.
      sxy3=0.
      sy4=0.
      sx3=0.
      sx2y=0.
      sxy2=0.
      sy3=0.
      sx2=0.
      sxy=0.
      sy2=0.
      sx=0.
      sy=0.
      sx2z=0.
      sxyz=0.
      sy2z=0.
      sxz=0.
      syz=0.
      sz=0.
      do 100 i=1,n
       sx4=sx4+x(i)*x(i)*x(i)*x(i)
       sx3y=sx3y+x(i)*x(i)*x(i)*y(i)
       sx2y2=sx2y2+x(i)*x(i)*y(i)*y(i)
       sxy3=sxy3+x(i)*y(i)*y(i)*y(i)
       sy4=sy4+y(i)*y(i)*y(i)*y(i)
       sx3=sx3+x(i)*x(i)*x(i)
       sx2y=sx2y+x(i)*x(i)*y(i)
       sxy2=sxy2+x(i)*y(i)*y(i)
       sy3=sy3+x(i)*y(i)*y(i)
       sx2=sx2+x(i)*x(i)
       sxy=sxy+x(i)*y(i)
       sy2=sy2+y(i)*y(i)
       sx=sx+x(i)
       sy=sy+y(i)
       sx2z=sx2z+x(i)*x(i)*z(i)
       sxyz=sxyz+x(i)*y(i)*z(i)
       sy2z=sy2z+y(i)*y(i)*z(i)
       sxz=sxz+x(i)*z(i)
       syz=syz+y(i)*z(i)
       sz=sz+z(i)
  100 continue
c
      a(1,1)=sx4
      a(2,1)=sx2y2
      a(3,1)=sx3y
      a(4,1)=sx3
      a(5,1)=sx2y
      a(6,1)=sx2
      a(1,2)=sx2y2
      a(2,2)=sy4
      a(3,2)=sxy3
      a(4,2)=sxy2
      a(5,2)=sy3
      a(6,2)=sy2
      a(1,3)=sx3y
      a(2,3)=sxy3
      a(3,3)=sx2y2
      a(4,3)=sx2y
      a(5,3)=sxy2
      a(6,3)=sxy
      a(1,4)=sx3
      a(2,4)=sxy2
      a(3,4)=sx2y
      a(4,4)=sx2
      a(5,4)=sxy
      a(6,4)=sx
      a(1,5)=sx2y
      a(2,5)=sy3
      a(3,5)=sxy2
      a(4,5)=sxy
      a(5,5)=sy2
      a(6,5)=sy
      a(1,6)=sx2
      a(2,6)=sy2
      a(3,6)=sxy
      a(4,6)=sx
      a(5,6)=sy
      a(6,6)=real(n)
      b(1)=sx2z
      b(2)=sy2z
      b(3)=sxyz
      b(4)=sxz
      b(5)=syz
      b(6)=sz
c
      call lax(a,6,6,b,1.e-6,1,is,vw,ip,icon)
c
      pxns=(2.*b(1)*b(2)-sqrt(4.*b(1)*b(1)*b(2)*b(2)-b(3)*b(3)
     -   *(b(3)*b(4)-2.*b(1)*b(5))))/b(3)/b(3)
      pxew=-(b(3)*pxns+b(4))/2./b(1)
      corf=b(1)*pxew*pxew+b(2)*pxns*pxns+b(3)*pxew*pxns+b(4)*pxew
     -     +b(5)*pxns+b(6)
c
      return
      end
***********************************************************************
*     proj17 - projection correction for a snapshot map at 17 GHz     *
*         solar north pole is to the top of the output image          *
*         intensity is normalized to disk brightness = 10000.         *
*         intensity information from nearby 3x3 pixels                *
*         pixel size of the output map = 4.91104 x 4.91104 arcsec^2   *
*           (identical to the SXT half resolution images)             *
*       pmat  : input,real(4) projection matrix                       *
*       solp  : input,real    polar angle of the sun (degree)         *
*       dskbr : input,real    abs(dskbr):disk brightness of           *
*                             the snapshot map                        *
*               if dskbr<0., main beam correction are done            *
*       pxew  : input,real x-offset of the disk center                *
*                 on the original map                                 *
*       pxns  : input,real l-offset of the disk center                *
*                 on the original map                                 *
*       xoffp : input,real  x-offset of the center                    *
*                 of the projected map from the disk center           *
*       yoffp : input,real l-offset of the center                     *
*                 of the projected map from the disk center           *
*       map   : input,real(512,512)  snapshot map                     *
*               output,real(512,512) projected map                    *
*                                        June 11, 1998  Y. Hanaoka    *
***********************************************************************
c
      subroutine proj17(pmat,solp,dskbr,pxew,pxns,xoffp,yoffp,map)
c
      real map(512,512),mapp(512,512),pmat(4)
c
      rd=57.29577951
      pi=3.141592653
c
*****
c
      x=(cos(solp/rd)*xoffp
     -     -sin(solp/rd)*yoffp)*4.91104/4.64947
      y=(sin(solp/rd)*xoffp
     -     +cos(solp/rd)*yoffp)*4.91104/4.64947
      pxx=pmat(1)*x+pmat(3)*y
      pyy=1.000112494*(pmat(2)*x+pmat(4)*y)
      do 500 l=1,512
      do 500 i=1,512
       x=(cos(solp/rd)*(real(i)-256.5)
     -     -sin(solp/rd)*(real(l)-256.5))*4.91104/4.64947
       y=(sin(solp/rd)*(real(i)-256.5)
     -     +cos(solp/rd)*(real(l)-256.5))*4.91104/4.64947
       xx=pmat(1)*x+pmat(3)*y
       yy=1.000112494*(pmat(2)*x+pmat(4)*y)
       ix=int(xx+257.5)
       il=int(yy+257.5)
       if (ix.ge.1 .and. ix.le.512 .and. il.ge.1 .and.
     -     il.le.512) then
        rx=amod(xx+pxx+pxew+256.5+1024.,512.)+0.5
        ry=amod(yy+pyy+pxns+256.5+1024.,512.)+0.5
        ix=int(rx+0.5)
        il=int(ry+0.5)
        sum=0.
        swt=0.
        do 510 ll=il-1,il+1
        do 510 ii=ix-1,ix+1
         wt=exp(-((rx-real(ii))*(rx-real(ii))+(ry-real(ll))
     -       *(ry-real(ll)))/0.5/0.5)
         swt=swt+wt
         sum=sum+wt*map(mod(ii+511,512)+1,
     -        mod(ll+511,512)+1)
  510   continue
* main beam correction factor *
        if (dskbr.lt.0.) then
         rpix2=((256.5-xoffp-real(i))*(256.5-xoffp-real(i))+
     -         (256.5-yoffp-real(l))*(256.5-yoffp-real(l)))*
     -         (4.91104/5556.5)*(4.91104/5556.5)
         g=exp(-3.141592*rpix2)
        else
         g=1.
        end if
        mapp(i,l)=sum/swt/g
       else
        mapp(i,l)=0.
       end if
  500 continue
c
      adskbr=abs(dskbr)
      do 520 l=1,512
      do 520 i=1,512
  520 map(i,l)=mapp(i,l)/adskbr*10000.
c
      return
      end
*************************************************************************
*     msun34 - model sun at 34 GHz                                      *
*       pmat   : input,real(4) projection matrix                        *
*       solr   : input,real    radius of the model sun (arcsec)         *
*       nantst : input,integer*(84)  antenna status                     *
*       pxew   : input real x-offset of the center                      *
*       pxns   : input,real y-offset of the center                      *
*       mapd   : output,real(512,512) dirty map                         *
*       mapc   : output,real(512,512) clean map                         *
*       mstype : output,integer       type of model sun                 *
*                                      November 22,1995  Y. Hanaoka     *
*************************************************************************
c
      subroutine msun34(pmat,solr,
     -     nantst,pxew,pxns,mapd,mapc,mstype)
c
      integer nantst(1),mshis(4)
      real mapd(512,512),mapc(1024,1024),
     -     mapa(512,512),pmat(4)
c
      rd=57.29577951
      pi=3.141592654
c
***** model sun at 34GHz *****
       rsun=solr/4.64947*33.8/17.
c
       do 510 l=1,1024
       do 510 i=1,1024
        xi=amod(real(i)-pxew+1024.,1024.)-513.
        xl=amod(real(l)-pxns+1024.,1024.)-513.
        xc=(pmat(4)*xi-pmat(3)*xl/1.000112494)
     -     /(pmat(1)*pmat(4)-pmat(3)*pmat(2))
        yc=(-pmat(2)*xi+pmat(1)*xl/1.000112494)
     -     /(pmat(1)*pmat(4)-pmat(3)*pmat(2))
        if (xc*xc+yc*yc .le. rsun*rsun) then
         mapc(i,l)=100.
        else
         mapc(i,l)=0.
        end if
  510  continue
c
***** histogram of model sun *****
c
      do 540 i=1,4
  540 mshis(i)=0
      do 550 l=1,512
      do 550 i=1,512
       j=int((mapc(i,l)+mapc(i+512,l)+mapc(i,l+512)+mapc(i+512,l+512))
     -    +1.)/100+1
       mshis(j)=mshis(j)+1
  550 continue
c
      if (mshis(1) .ge. mshis(2)/3.) mstype=1
      if ((mshis(1) .lt. mshis(2)/3.) .and. (mshis(1).gt.mshis(3)))
     -  mstype=2
      if ((mshis(3) .lt. mshis(2)/3.) .and. (mshis(3).gt.mshis(1)))
     -  mstype=3
      if (mshis(3) .ge. mshis(2)/3.) mstype=4
ccc   write(6,'('' mstype ='',i4)') mstype
c
***** dirty disk *****
c
      do 650 l=1,512
      do 650 i=1,512
  650 mapd(i,l)=mapc(i,l)+mapc(i+512,l)
     -          +mapc(i,l+512)+mapc(i+512,l+512)
c
      call cft512(mapd,mapa,1)
      call uvobs(nantst,1.,mapd,mapa)
      call cft512(mapd,mapa,-1)
c
      return
      end
***********************************************************************
*     clean34 - clean process for 34 GHz images                       *
*         model source :  gaussian source and flat disk               *
*         adaptive clean-criterion definition available               *
*       menu   : input,integer                                        *
*                  2 : disk/pointsource(+) clean with restoration     *
*                 -2 : disk/pointsource(+/-) clean with restoration   *
*                 22 : disk/pointsource(+) clean with point source    *
*                        restoration, without disk restoration        *
*                -22 : disk/pointsource(+/-) clean with point source  *
*                        restoration, without disk restoration        *
*                 -1 : pointsource(+/-) clean with restore            *
*                  1 : pointsource(+) clean with restore              *
*                  3 : preclean only, criter=dskbr*5.+bgr             *
*       map    : input,real(512,512)  dirty map                       *
*                output,real(512,512) resudual                        *
*       pmat   : input,real(4) projection matrix                      *
*       solr   : input,real    radius of the model sun (arcsec)       *
*       nantst : input,integer*(84)  antenna status                   *
*       gbeam  : input,real(512,512,0:10) dirty beam of gaussian      *
*                   sources                                           *
*       criter : input,real  (lower limit of) clean criterion         *
*                 if criter < 0, -criter*dskbr is actually adopted    *
*       crfac  : input real  factor for clean criterion               *
*                 actual clean criterion is the larger one of         *
*                 criter or crfac*(maximum of the map)                *
*       cr     : output,real actual clean criterion                   *
*       pxew   : input/output,real x-offset of disk center (pixel)    *
*       pxns   : input/output,real y-offset of disk center (pixel)    *
*                  only valid for disk clean menu                     *
*                  else when pxew=999. or pxns=999., x-offset and     *
*                  y-offset is calculated and output pxew and pxns    *
*                  is set to calculated offset.                       *
*       nmax   : input,integer  maximum number of clean components    *
*              : output,integer number of clean components            *
*       bgr    : output,real   background brightness                  *
*       dskbr  : input,real    disk brightness (valid for criter < 0  *
*                              and menu=+/-1,+/-11                    *
*              : output,real   disk brightness                        *
*       pi     : output,integer(nmax) x-coordinate of clean points    *
*       pl     : output,integer(nmax) y-coordinate of clean points    *
*       pw     : output,integer(nmax) width of clean points           *
*       pa     : output,real(nmax)    intensity of clean points       *
*       cmap   : output,real(1024,1024) clean map                     *
*                                        June 12 ,1998  Y. Hanaoka    *
***********************************************************************
c
      subroutine clean34(menu,map,pmat,solr,
     -     nantst,gbeam,criter,crfac,cr,
     -     pxew,pxns,corf,nmax,bgr,dskbr,pi,pl,pw,pa,cmap)
c
      integer nantst(84),pi(nmax),pl(nmax),pw(nmax)
      real map(512,512),cmap(1024,1024),
     -     ddisk(512,512),cdisk(1024,1024),
     -     cbeam(-30:30,-30:30),gbeam(512,512,0:10),
     -     pmat(4),pa(nmax)
c
      n=0
c
      call mean(1,map,am,amxmap,imx,lmx)
ccc   write(6,'('' maximum of the map :'',f12.1)') amxmap
c
      if (mod(iabs(menu),10) .ge. 2) then
c
*********************************************
*     preclean process for disk clean       *
*********************************************
c
***** model sun *****
       call msun34(pmat,solr,
     -      nantst,0.,0.,ddisk,cdisk,mstype)
c
       call clhis(map,mstype,40,bgr,dbr)
       dskbr=dbr-bgr
ccc    write(6,'('' background/disk(initial) :'',2f10.2)')
ccc  -       bgr,dskbr
       if (criter.lt.0.) then
        cr=amax1(amxmap*crfac,dskbr*abs(criter))
       else
        cr=amax1(amxmap*crfac,criter)
       endif
       cr=amax1(dskbr*5.+bgr,cr*2.)
ccc    write(6,'('' preclean criterion :'',f10.1)') cr
c
       do 100 while(.true.)
c
***** search maximum *****
        call mean(menu,map,am,amx,imx,lmx)
        if (amx.gt.cr .and. n.lt.nmax) then
c
***** subtract point source component *****
         n=n+1
         pi(n)=imx
         pl(n)=lmx
         call cleanp(map,gbeam,0.02,amx,imx,lmx,pw(n),pa(n))
c        write(6,'('' '',i6,'': am,amx,w,a,imx,lmx ='',
c     -       f10.2,f14.2,i3,f10.2,2i5)')
c     -       n,am,amx,pw(n),pa(n),pi(n),pl(n)
c
        else
         go to 110
        end if
c
  100  continue
  110  if (n.ne.0) then
ccc     write(6,'('' n ='',i6,'': am,amx,w,a,imx,lmx ='',
ccc  -       f10.2,f14.2,i3,f10.2,2i5)')
ccc  -       n,am,amx,pw(n),pa(n),pi(n),pl(n)
       else
ccc     write(6,'('' n ='',i6)') n
       end if
c
*********************************************
*     disk clean process                    *
*********************************************
c
       call clhis(map,mstype,40,bgr,dbr)
       dskbr=dbr-bgr
ccc    write(6,'('' background/disk          :'',2f10.2)')
ccc  -       bgr,dskbr
c
c       if (mod(iabs(menu),10).eq.3) return
c
***** subtract disk component *****
c
       do 140 l=1,512
       do 140 i=1,512
  140  map(i,l)=map(i,l)-bgr
c
***** calculate disk center position ***** 
       if (pxew.eq.999. .or. pxns.eq.999.) then
        call possh34(map,dskbr,ddisk,pxew,pxns,corf)
c        write(6,'(''  pxew,pxns,corf ='',3f8.3)') pxew,pxns,corf
       end if
c
       if (mod(iabs(menu),10).eq.3) return
c
***** model sun with offset *****
       call msun34(pmat,solr,
     -      nantst,pxew,pxns,ddisk,cdisk,mstype)
c
***** subtract disk component *****
       do 150 l=1,512
       do 150 i=1,512
        map(i,l)=map(i,l)-dskbr*ddisk(i,l)/100.
 150   continue
c
      end if
c
c
*********************************************
*     clean process                         *
*********************************************
c
       if (criter.lt.0.) then
        cr=amax1(amxmap*crfac,dskbr*abs(criter))
       else
        cr=amax1(amxmap*crfac,criter)
       endif
ccc    write(6,'('' clean criterion    :'',f10.1)') cr
c
       do 160 while(.true.)
c
***** search maximum *****
        call mean(menu,map,am,amx,imx,lmx)
        if (abs(amx).gt.cr .and. n.lt.nmax) then
c
***** subtract point source component *****
         n=n+1
         pi(n)=imx
         pl(n)=lmx
         call cleanp(map,gbeam,0.02,amx,imx,lmx,pw(n),pa(n))
c        write(6,'('' '',i6,'': am,amx,w,a,imx,lmx ='',
c     -       f10.2,f14.2,i3,f10.2,2i5)')
c     -       n,am,amx,pw(n),pa(n),pi(n),pl(n)
c
        else
         go to 170
        end if
c
  160  continue
  170  if (n.ne.0) then
ccc     write(6,'('' n ='',i6,'': am,amx,w,a,imx,lmx ='',
ccc  -       f10.2,f14.2,i3,f10.2,2i5)')
ccc  -       n,am,amx,pw(n),pa(n),pi(n),pl(n)
       else
ccc     write(6,'('' n ='',i6)') n
       end if
c
      nmax=n
      if (abs(menu)/10.eq.1) return
c
c
*********************************************
*     restore process                       *
*********************************************
c
      do 240 l=1,512
      do 240 i=1,512
  240 ddisk(i,l)=0.
c
      do 500 iw=0,10
c
***** create clean beam *****
c
       if (iw.eq.0) then
        do 210 l=-30,30
        do 210 i=-30,30
  210   cbeam(i,l)=0.
        cbeam(0,0)=100.
       else
        sum=0.
        do 220 l=-30,30
        do 220 i=-30,30
         cbeam(i,l)=100.*exp(-(real(i)*real(i)+real(l)*real(l))
     -     /real(iw)/real(iw))
         sum=sum+cbeam(i,l)
  220   continue
        do 230 l=-30,30
        do 230 i=-30,30
  230   cbeam(i,l)=100./sum*cbeam(i,l)
       end if   
c
***** restore point sources *****
c
       do 250 ip=1,n
        if (pw(ip).eq.iw) then
         do 260 l=-30,30
         do 260 i=-30,30
c  260    ddisk(pi(ip)+i,pl(ip)+l)=ddisk(pi(ip)+i,pl(ip)+l)+pa(ip)
c     -       *cbeam(i,l)/gbeam(257,257,iw)
          ii=mod(pi(ip)+i+511,512)+1
          ll=mod(pl(ip)+l+511,512)+1
          ddisk(ii,ll)=ddisk(ii,ll)
     -      +pa(ip)*cbeam(i,l)/gbeam(257,257,iw)
  260    continue
c        write(6,'('' '',i4,5f10.3)') iw,(cbeam(ii,0),ii=-10,10,5)
        end if
  250  continue
c
  500 continue
c
***** restore disk component *****
c
      cfact=0.
      if (iabs(menu).eq.2) cfact=dskbr
      do 300 l=1,1024
      do 300 i=1,1024
       cdisk(i,l)=ddisk(mod(i+511,512)+1,mod(l+511,512)+1)
     -             +cfact*cdisk(i,l)/100.
  300 continue
c
***** convolution by gaussian beam *****
c
       do 340 ll=-10,10
       do 340 ii=-10,10
  340  cbeam(ii,ll)=0.07461*exp(-0.234403
     -        *(real(ll)*real(ll)+real(ii)*real(ii)))
c
      do 310 l=1,1024
      do 310 i=1,1024
  310 cmap(i,l)=0.
      do 320 l=1,1024
      do 320 i=1,1024
       do 330 ll=-10,10
       do 330 ii=-10,10
        ib=mod(i+ii+1023,1024)+1
        lb=mod(l+ll+1023,1024)+1
        cmap(i,l)=cmap(i,l)+cdisk(ib,lb)*cbeam(ii,ll)
  330  continue
  320 continue
c
***** add residual *****
      do 350 l=1,1024
      do 350 i=1,1024
       cmap(i,l)=cmap(i,l)+map(mod(i+511,512)+1,mod(l+511,512)+1)
  350 continue
c
      return
      end
***********************************************************************
*     possh34 - position shift between dirty model disk               *
*              and disk on dirty map                                  *
*       map    : input,real(512,512)  dirty map                       *
*       dskbr  : input,real    brightness of the dirty disk           *
*       disk   : input,real(512,512)  dirty map of model disk         *
*       pxew   : output,real x-offset of disk center                  *
*       pxns   : output,real y-offset of disk center                  *
*       corf   : output,real correlation between dirty map and model  *
*                                        June 12 ,1998  Y. Hanaoka    *
***********************************************************************
c
      subroutine possh34(map,dskbr,disk,pxew,pxns,corf)
c
      real map(512,512),disk(512,512),mapt(512,512),
     -     corfma(16,16),xp(25),yp(25),zc(25)
c
      pi=3.14159265
c
***** calculate correlation and search maximum correlation *****
c
      sa=0.
      saa=0.
      do 220 l=1,512
      do 220 i=1,512
       mapt(i,l)=amin1(dskbr*3.1,amax1(dskbr*(-0.2),map(i,l)))
       sa=sa+mapt(i,l)
       saa=saa+mapt(i,l)*mapt(i,l)
  220 continue
c
      sb=0.
      sbb=0.
      do 280 l=1,512
      do 280 i=1,512
       sb=sb+disk(i,l)
       sbb=sbb+disk(i,l)*disk(i,l)
  280 continue
c
***** pass-1 *****
c
      corfmx=0.
      do 300 loff=-240,240,32
      do 300 ioff=-240,240,32
       sab=0.
       do 310 l=1,512
       do 310 i=1,512
        ib=mod(i-ioff+511,512)+1
        lb=mod(l-loff+511,512)+1
        sab=sab+mapt(i,l)*disk(ib,lb)
  310  continue
       corf=(512.*512.*sab-sa*sb)/sqrt((512.*512.*saa-sa*sa)
     -      *(512.*512.*sbb-sb*sb))
       if (corf.gt.corfmx) then
        corfmx=corf
        isx=ioff
        isl=loff
       end if
       corfma((ioff+240)/32+1,(loff+240)/32+1)=corf
  300 continue
c
c      write(6,'('' corfm ='')')
c      write(6,'('' '',11f8.5)') corfma
c
ccc   write(6,'(''  pass-1: isx  ='',i6,''   isl  ='',i6,
ccc  -''   correlation ='',f8.4)') isx,isl,corfmx
c
***** pass-2 *****
c
      corfmx=0.
      do 350 loff=isl-36,isl+36,6
      do 350 ioff=isx-36,isx+36,6
       sab=0.
       do 360 l=1,512
       do 360 i=1,512
        ib=mod(i-ioff+511,512)+1
        lb=mod(l-loff+511,512)+1
        sab=sab+mapt(i,l)*disk(ib,lb)
  360  continue
       corf=(512.*512.*sab-sa*sb)/sqrt((512.*512.*saa-sa*sa)
     -      *(512.*512.*sbb-sb*sb))
       if (corf.gt.corfmx) then
        corfmx=corf
        ioffmx=ioff
        loffmx=loff
       end if
       corfma((ioff-isx+36)/6+1,(loff-isl+36)/6+1)=corf
  350 continue
c
c      write(6,'('' corfm ='')')
c      write(6,'('' '',11f8.5)') corfma
c
      isx=ioffmx
      isl=loffmx
ccc   write(6,'(''  pass-2: isx  ='',i6,''   isl  ='',i6,
ccc  -''   correlation ='',f8.4)') isx,isl,corfmx
c
***** pass-3 *****
c
      corfmx=0.
      do 370 loff=isl-5,isl+5
      do 370 ioff=isx-5,isx+5
       sab=0.
       do 380 l=1,512
       do 380 i=1,512
        ib=mod(i-ioff+511,512)+1
        lb=mod(l-loff+511,512)+1
        sab=sab+mapt(i,l)*disk(ib,lb)
  380  continue
       corf=(512.*512.*sab-sa*sb)/sqrt((512.*512.*saa-sa*sa)
     -      *(512.*512.*sbb-sb*sb))
       if (corf.gt.corfmx) then
        corfmx=corf
        ioffmx=ioff
        loffmx=loff
       end if
       corfma(ioff-isx+5+1,loff-isl+5+1)=corf
  370 continue
c
c      write(6,'('' corfm ='')')
c      write(6,'('' '',11f8.5)') corfma
c
      isx=ioffmx
      isl=loffmx
ccc   write(6,'(''  pass-3: isx  ='',i6,''   isl  ='',i6,
ccc  -''   correlation ='',f8.4)') isx,isl,corfmx
c
***** pass-4 *****
c
      do 400 loff=isl-2,isl+2
      do 400 ioff=isx-2,isx+2
       sab=0.
       do 410 l=1,512
       do 410 i=1,512
        ib=mod(i-ioff+511,512)+1
        lb=mod(l-loff+511,512)+1
        sab=sab+mapt(i,l)*disk(ib,lb)
  410  continue
       corf=(512.*512.*sab-sa*sb)/sqrt((512.*512.*saa-sa*sa)
     -      *(512.*512.*sbb-sb*sb))
c       write(6,'('' corf ='',f8.3)') corf
       corfma(ioff-isx+2+1,loff-isl+2+1)=corf
  400 continue
c
c      write(6,'('' corfm ='')')
c      write(6,'('' '',11f8.5)') corfma
c
      do 450 i=1,25
  450 xp(i)=real(mod(i-1,5)-2)
c      write(6,'('' '',5f10.5)') xp
      do 460 i=1,25
  460 yp(i)=real(((i-1)/5)-2)
c      write(6,'('' '',5f10.5)') yp
      do 470 i=1,25
  470 zc(i)=corfma(mod(i-1,5)+1,(i-1)/5+1)
c      write(6,'('' '',5f10.5)') zc
c
      call getctr(xp,yp,zc,25,pxew,pxns,corf)
c      write(6,'(''  pass-4: pxew ='',f8.3,''   pxns ='',f8.3)')
c     -      pxew,pxns
      pxew=pxew+real(isx)
      pxns=pxns+real(isl)
ccc   write(6,'('' disk center position: pxew ='',f8.3,''   pxns ='',
ccc  -     f8.3,''   correlation ='',f8.4)') pxew,pxns,corf
c
      return
      end
***********************************************************************
*     proj34 - projection correction for a snapshot map               *
*         solar north pole is to the top of the output image          *
*         intensity is normalized to disk brightness = 10000.         *
*         intensity information from nearby 3x3 pixels                *
*         pixel size of the output map = 4.91104 x 4.91104 arcsec^2   *
*           (identical to the SXT half resolution images)             *
*       pmat  : input,real(4) projection matrix                       *
*       solp  : input,real    polar angle of the sun (degree)         *
*       dskbr : input,real    abs(dskbr):disk brightness of           *
*                             the snapshot map                        *
*               if dskbr<0., main beam correction are done            *
*       pxew  : input,real x-offset of the disk center                *
*                 on the original map                                 *
*       pxns  : input,real l-offset of the disk center                *
*                 on the original map                                 *
*       xoffp : input,real  x-offset of the center                    *
*                 of the projected map from the disk center           *
*       yoffp : input,real l-offset of the center                     *
*                 of the projected map from the disk center           *
*       map   : input,real(512,512)  snapshot map                     *
*                                        June 13, 1998  Y. Hanaoka    *
***********************************************************************
c
      subroutine proj34(pmat,solp,dskbr,pxew,pxns,xoffp,yoffp,map)
c
      real map(1024,1024),mapp(1024,1024),pmat(4)
c
      rd=57.29577951
      pi=3.141592653
c
*****
c
      afact=(4.91104/2.)/(4.64947/33.8*17.)
      x=(cos(solp/rd)*xoffp
     -     -sin(solp/rd)*yoffp)*afact
      y=(sin(solp/rd)*xoffp
     -     +cos(solp/rd)*yoffp)*afact
      pxx=pmat(1)*x+pmat(3)*y
      pyy=1.000112494*(pmat(2)*x+pmat(4)*y)
      do 500 l=1,1024
      do 500 i=1,1024
       x=(cos(solp/rd)*(real(i)-512.5)
     -     -sin(solp/rd)*(real(l)-512.5))*afact
       y=(sin(solp/rd)*(real(i)-512.5)
     -     +cos(solp/rd)*(real(l)-512.5))*afact
       xx=pmat(1)*x+pmat(3)*y
       yy=1.000112494*(pmat(2)*x+pmat(4)*y)
       ix=int(xx+513.5)
       il=int(yy+513.5)
       if (ix.ge.1 .and. ix.le.1024 .and. il.ge.1 .and.
     -     il.le.1024) then
        rx=amod(xx+pxx+pxew+512.5+2048.,1024.)+0.5
        ry=amod(yy+pyy+pxns+512.5+2048.,1024.)+0.5
        ix=int(rx+0.5)
        il=int(ry+0.5)
        sum=0.
        swt=0.
        do 510 ll=il-1,il+1
        do 510 ii=ix-1,ix+1
         wt=exp(-((rx-real(ii))*(rx-real(ii))+(ry-real(ll))
     -       *(ry-real(ll)))/0.5/0.5)
         swt=swt+wt
         sum=sum+wt*map(mod(ii+1023,1024)+1,
     -        mod(ll+1023,1024)+1)
  510   continue
* main beam correction factor *
        if (dskbr.lt.0.) then
         rpix2=((512.5-xoffp-real(i))*(512.5-xoffp-real(i))+
     -         (512.5-yoffp-real(l))*(512.5-yoffp-real(l)))*
     -         (2.45553/5556.5)*(2.45553/5556.5)
         g=exp(-3.141592*rpix2)
        else
         g=1.
        end if
        if (swt.lt.1.e-5) write(6,'('' ix,il ='',2i5)') ix,il
        mapp(i,l)=sum/swt/g
       else
        mapp(i,l)=0.
       end if
  500 continue
c
      adskbr=abs(dskbr)
      do 520 l=1,1024
      do 520 i=1,1024
  520 map(i,l)=mapp(i,l)/adskbr*10000.
c
      return
      end
***********************************************************************
*     cleane  - subtract extended source  Ver.2                       *
*       map    : input,real(512,512)  dirty map                       *
*              : output,real(512,512) dirty map, extended source      *
*                                            subtracted               *
*       nantst : input,integer(84) antenna status                     *
*       amx    : input,real  maximum of map                           *
*       gaine  : input,real  CLEAN loop gain                          *
*       emap   : input,real(512,512)  extended source                 *
*       scale  : output,real  scaling factor for extended source      * 
*                                         Jun  13,1995    koshix      *
***********************************************************************
          subroutine cleane(map,nantst,amx,gaine,emap,scale)

          integer nantst(84)
          real gaine,amx,dmx,scale
          real map(512,512),emap(512,512)
          real a(512,512),b(512,512)


          do 200 i=1,512
             do 100 j=1,512

                a(i,j)=emap(i,j) 
                b(i,j)=0. 

100          continue
200       continue

          call cft512(a,b,1)
          call uvobs(nantst,1.,a,b)
          call cft512(a,b,-1)

          dmx=0.

          do 400 i=1,512
             do 300 j=1,512

                if(abs(a(i,j)) .gt. abs(dmx)) dmx=a(i,j)

300          continue
400       continue                

          do 600 i=1,512
             do 500 j=1,512

                map(i,j)=map(i,j)-a(i,j)/abs(dmx)*abs(amx)*gaine

500          continue
600       continue  

          scale=abs(abs(amx)*gaine/abs(dmx)) 

          return
          end
***********************************************************************
*     cleanf - CLEAN process for extended source                      *
*       model source :  gaussian source and flat disk                 *
*       menu   : input,integer                                        *
*                  2 : disk/pointsource(+)    CLEAN with restore      *
*                 22 : disk/pointsource(+) clean with point source    *
*                  1 : pointsource(+)         CLEAN with restore      *
*                 -1 : pointsource(+/-)       CLEAN with restore      *
*                 -2 : disk/pointsource(+/-)  CLEAN with restore      *
*                -22 : disk/pointsource(+/-) clean with point source  *
*       tmenu  :   1 : search single trim                             *
*              :   2 : search multi trim                              *   
*       map    : input,real(512,512)  dirty map                       *
*       ddisk  : input,real(512,512)  dirty map of model disk         *
*       cdisk  : input,real(512,512)  CLEAN map of model disk         *
*       gaine  : input,real  CLEAN loop gain                          *
*       trim   : input,real  CLEAN trim level                         *
*       criter : input,real  CLEAN criterion                          *
*       nantst : input,integer(84) antenna status                     *
*       nmax   : input,integer  maximum number of CLEAN components    *
*       mstype : input,integer  type of model sun                     *
*       bg     : output,real   background brightness                  *
*       dskbr  : output,real   disk brightness                        *
*       extend : output,real(512,512) CLEAN points                    *
*                                             for extended source     *
*       cmap   : output,real(512,512) CLEAN map                       *
*                                         Jun  13,1995     koshix     *
*                                       January 9, 1999 Y. Hanaoka    *
***********************************************************************

      subroutine cleanf(menu,tmenu,map,ddisk,cdisk,gaine,trim,
     -                  criter,nantst,nmax,mstype,
     -                  bg,dskbr,extend,cmap)

      integer menu,tmenu,nmax,ne,mstype  
      integer nantst(84)
      real gaine,trim,cr,criter,bg,dskbr,cratio
      real scale
      real map(512,512),cmap(512,512),ddisk(512,512),cdisk(512,512)
      real cbeam(-30:30,-30:30)
      real emap(512,512),remap(512,512),map_status(512,512)
      real point(512,512),extend(512,512)


      ne=0


      do 200 l=1,512
         do 100 i=1,512

            remap(i,l)=0.
            extend(i,l)=0.

100      continue
200   continue


      if (mod(iabs(menu),10) .ge. 2) then

*********************************************
*     CLEAN process                         *
*********************************************

         call clhis(map,mstype,300,bg,db)  
       
         dskbr=db-bg

ccc      write(6,'('' background/disk(initial) :'',
ccc  -                2f10.2)') bg,dskbr
 
         if(criter .lt. 0.) then
     
            cr=amax1(dskbr*5.+bg,dskbr*(-criter))
ccc         write(6,'('' preclean criterion :'',f10.1)') cr
   
         else
     
            cr=amax1(dskbr*5.+bg,criter)
    
         endif

         do 500 while(.true.)

***** search maximum *****
            call mean(1,map,am,amx,imx,lmx)
            
            if(amx.gt.cr .and. ne.lt.nmax) then

***** contour trim *****

               if(tmenu .eq. 1) then

                  call contour_trim1(map,gaine,trim,imx,lmx,npix,          
     -                               map_status,emap)          

               endif
               if(tmenu .eq. 2) then

                  call contour_trim2(map,gaine,trim,imx,lmx,npix,          
     -                               map_status,emap)          

               endif

            
***** subtract extended source component *****

               ne=ne+npix
               call cleane(map,nantst,amx,gaine,emap,scale)
c              write(6,'('' '',i6,'': amx,imx,lmx,npix,scale ='',
c     -              f14.2,2i5,i6,f5.2)')ne,amx,imx,lmx,npix,scale

               do 400 l=1,512
                  do 300 i=1,512
 
                     remap(i,l)=remap(i,l)+emap(i,l)*scale

                     extend(i,l)=extend(i,l)+real(map_status(i,l))

300               continue
400            continue 

            else
    
               go to 600
    
            endif

500      continue


600      if(ne.ne.0) then
      
ccc         write(6,'('' ne ='',i6,''   :   amx ='',f11.2)')ne,amx
       
         else
      
ccc         write(6,'('' ne ='',i6)') ne
 
         endif


*********************************************
*     disk CLEAN process                    *
*********************************************

         call clhis(map,mstype,300,bg,db)  
         dskbr=db-bg
      
ccc      write(6,'('' background/disk          :'',
ccc  -                2f10.2)') bg,dskbr

***** subtract disk component *****

         do 800 l=1,512       
            do 700 i=1,512
  
               map(i,l)=map(i,l)-bg

700         continue
800      continue

         call clnd2(map,ddisk,0,0,dskbr)

      endif
      

*********************************************
*     CLEAN process                         *
*********************************************
      cratio = dskbr/10000.
      cratio=1.
      if(criter .lt. 0.) then
       
         cr=dskbr*(-criter)
ccc      write(6,'('' clean criterion    :'',f10.1)') cr
      else
     
c         cr=criter
         cr = criter*cratio
      endif
 
      do 1100 while(.true.)

***** search maximum *****

         call mean(menu,map,am,amx,imx,lmx)
  
         if(abs(amx).gt.cr .and. ne.lt.nmax) then

***** contour trim *****

            if(tmenu .eq. 1) then             

               call contour_trim1(map,gaine,trim,imx,lmx,npix,          
     -                            map_status,emap)          

            endif
            if(tmenu .eq. 2) then

               call contour_trim2(map,gaine,trim,imx,lmx,npix,          
     -                            map_status,emap)          

            endif

 
***** subtract extended source component *****

            ne=ne+npix
            call cleane(map,nantst,amx,gaine,emap,scale)
ccc         write(6,'('' '',i6,'': amx,amx/cratio,
ccc  -            imx,lmx,npix,scale ='',
ccc  -            2f14.2,2i5,i6,f7.2)')ne,amx,amx/cratio,
ccc  -            imx,lmx,npix,scale

            do 1000 l=1,512
               do 900 i=1,512

                  remap(i,l)=remap(i,l)+emap(i,l)*scale

                  extend(i,l)=extend(i,l)+real(map_status(i,l))

900            continue
1000        continue 

         else
    
            go to 1200
    
         endif

1100     continue

1200     if(ne.ne.0) then
      
ccc         write(6,'('' ne ='',i6,''   :   amx ='',f11.2)')ne,amx
       
         else
      
ccc         write(6,'('' ne ='',i6)') ne
 
         endif

         nmax=ne

*********************************************
*     restore process                       *
*********************************************

         do 1400 l=1,512
            do 1300 i=1,512
  
               cmap(i,l)=0.
               point(i,l)=0. 

1300        continue
1400     continue

***** restore disk component *****

      cfact=0.
      if (iabs(menu).eq.2) cfact=dskbr

      do 1600 l=1,512
         do 1500 i=1,512
             
            point(i,l)=cfact*cdisk(i,l)/100.
  
1500     continue
1600  continue 

***** restore extended source *****

      do 1800 l=1,512
         do 1700 i=1,512
      
            point(i,l)=point(i,l)+remap(i,l)
  
1700     continue
1800  continue

***** convolution by gaussian beam *****

      do 2000 ll=-10,10
         do 1900 ii=-10,10
   
            cbeam(ii,ll)=0.07461*exp(-0.234403
     -                  *(real(ll)*real(ll)+real(ii)*real(ii)))

1900     continue
2000  continue 
  
      do 2400 l=1,512
        do 2300 i=1,512
       
           do 2200 ll=-10,10
              do 2100 ii=-10,10
       
                 ib=mod(i+ii+511,512)+1
                 lb=mod(l+ll+511,512)+1
       
                 cmap(i,l)=cmap(i,l)+point(ib,lb)*cbeam(ii,ll)

2100           continue
2200        continue
 
2300     continue 
2400  continue

***** add residual *****

      do 2600 l=1,512
         do 2500 i=1,512
      
            cmap(i,l)=cmap(i,l)+map(i,l)
  
2500     continue
2600  continue

      return
      end
***********************************************************************
*     clnd2 - subtraction of disk component                           *
*       map    : input,real(512,512)  dirty map                       *
*              : output,real(512,512) disty map, dirty disk subtracted*
*       disk   : input,real(512,512)  dirty map of model disk         *
*       ioffm  : input,integer x-offset of disk center                *
*       loffm  : input,integer y-offset of disk center                *
*       dskbr  : output,real   disk brightness                        *
*                                         June 22,1992  Y. Hanaoka    *
***********************************************************************
c
      subroutine clnd2(map,disk,ioffm,loffm,dskbr)
c
      real map(512,512),disk(512,512)
c
      do 400 l=1,512
      do 400 i=1,512
       ib=mod(i-ioffm+511,512)+1
       lb=mod(l-loffm+511,512)+1
c       map(i,l)=map(i,l)-dskbr*disk(ib,lb)/cdisk
       map(i,l)=map(i,l)-dskbr*disk(ib,lb)/100.
 400  continue
c
      return
      end
*****************************************************************
*       contour_trim1 - serach contour trim  Ver.1              *
*                                                               *
*       map    : input,real(512,512)  dirty map                 *
*       gaine  : input,real  CLEAN loop gain                    *
*       trim   : input,real  CLEAN trim level                   *
*       imx    : input,integer                                  *
*                            x-coordinate of the maximum point  *
*       lmx    : input,integer                                  *
*                            y-coordinate of the maximum point  *
*       npix   : output,integer  number of pixels               *
*       map_status  :  output,integer(512,512)                  *
*                            CLEAN point for extended source    *
*       e_map       :  output,real(512,512)  extended component * 
*                                      Jun  13,1995   koshix    *
*****************************************************************

          subroutine contour_trim1(map,gaine,trim,imx,lmx,npix,          
     &                             map_status,emap)

          integer key
          integer i,j,k,l,imx,lmx,ib,lb
          integer npix
          integer n,buffer_num
          integer ibuffer(262144),jbuffer(262144)
          integer map_status(512,512) 
          real gaine,trim
          real map(512,512)
          real emap(512,512)
        

********** Contour Trim ********** 

          do 200 i=1,512
             do 100 j=1,512

                map_status(i,j)=0 
                emap(i,j)=0.

100          continue
200       continue

          do 300 n=1,262144

             ibuffer(n)=-1
             jbuffer(n)=-1

300       continue

          map_status(imx,lmx)=1
          buffer_num=1           
          ibuffer(buffer_num)=imx
          jbuffer(buffer_num)=lmx

          key=1

          do 700 while(key .eq. 1)          

             key=0

             do 600 n=1,buffer_num

                do 500 k=-1,1    
                   do 400 l=-1,1

                      ib=mod(ibuffer(n)+k+511,512)+1
                      lb=mod(jbuffer(n)+l+511,512)+1
   
                      if(abs(map(ib,lb)) .ge.
     &                   abs(map(imx,lmx))*trim  .and.
     &                   map_status(ib,lb) .ne. 1) then          

                         map_status(ib,lb)=1
                         buffer_num=buffer_num+1           
                         ibuffer(buffer_num)=ib
                         jbuffer(buffer_num)=lb 

                         key=1
  
                      endif      

400                continue
500             continue

600          continue

700       continue 

          npix=buffer_num

          do 900 i=1,512
             do 800 j=1,512

                if(map_status(i,j) .eq. 1) then

                   emap(i,j)=map(i,j)*gaine

                endif

800          continue
900       continue


          return
          end
*****************************************************************
*       contour_trim2 - serach contour trim  Ver.2              *
*                                                               *
*       map    : input,real(512,512)  dirty map                 *
*       gaine  : input,real  CLEAN loop gain                    *
*       trim   : input,real  CLEAN trim level                   *
*       imx    : input,integer                                  *
*                            x-coordinate of the maximum point  *
*       lmx    : input,integer                                  *
*                            y-coordinate of the maximum point  *
*       npix   : output,integer  number of pixels               *
*       map_status  :  output,integer(512,512)                  *
*                            CLEAN point for extended source    *
*       e_map       :  output,real(512,512)  extended component * 
*                                      Jun  13,1995   koshix    *
*****************************************************************

          subroutine contour_trim2(map,gaine,trim,imx,lmx,npix,          
     &                             map_status,emap)

          integer i,j,imx,lmx
          integer npix
          integer map_status(512,512) 
          real gaine,trim
          real map(512,512)
          real emap(512,512)
        

********** Contour Trim ********** 

          do 200 i=1,512
             do 100 j=1,512

                map_status(i,j)=0 
                emap(i,j)=0.

100          continue
200       continue

          npix=0          

          do 400 i=1,512
             do 300 j=1,512
   
                if(abs(map(i,j)) .ge.
     &             abs(map(imx,lmx))*trim) then

                   map_status(i,j)=1  
                   npix=npix+1

                endif      

300          continue
400       continue

          do 600 i=1,512
             do 500 j=1,512

                if(map_status(i,j) .eq. 1) then

                   emap(i,j)=map(i,j)*gaine

                endif

500          continue
600       continue

          return
          end
