program clean_component17
***********************************************************************
*     create clean component (2000x4 etc.)                            *
*                           ver 6.2 October 10, 2007  H. Koshiishi    *
***********************************************************************
c
      integer naxes(2),nantst(84),iimg(1000),pir(2000),pil(2000),
     -     plr(2000),pll(2000),pwr(2000),pwl(2000)
      integer pxr(2000),pyr(2000),pxl(2000),pyl(2000)
      integer*2 icald(336)
      real gbeam(512,512,0:10),flxr(84),flxl(84),
     -     mapra(512,512),maprb(512,512),mapla(512,512),maplb(512,512),
     -     cmapr(512,512),cmapl(512,512),outmap(8000),pmat(4),
     -     par(2000),pal(2000)
      real pacr(2000),pacl(2000),pr(2000,4),pl(2000,4)
      real rlph(84)
      character*80 filehd,infpar,outf,outfft
      character*18 dattyp
      character*17 datetime
      character*1  fform
c
ccc   write(6,'(/'' input directory for output files       : ''$)')
      read(5,'(a)') outf
ccc   write(6,'('' input rawdata filename header          : ''$)')
      read(5,'(a)') filehd
ccc   write(6,'('' input parameter filename               : ''$)')
      read(5,'(a)') infpar
c
ccc   write(6,'('' input pixel-size of the output files   : ''$)')
      read(5,*) npix
c
ccc   write(6,'('' input offset of the image center'')')
ccc   write(6,'('' unit(1:arcsec/2:pixel) & x-/y-offsets  : ''$)')
      read(5,*) nuoff,xoffa,yoffa
      if (nuoff .eq. 1) then
       xoffp=xoffa/4.91106
       yoffp=yoffa/4.91106
      else
       xoffp=xoffa
       yoffp=yoffa
       xoffa=xoffp*4.91106
       yoffa=yoffp*4.91106
      endif
c
ccc   write(6,'('' input first frame number               : ''$)')
      read(5,*) nfrst
      if (nfrst.ne.0) then
ccc    write(6,'('' input last frame number                : ''$)')
       read(5,*) nfrend
ccc    write(6,'('' input frame interval                   : ''$)')
       read(5,*) itvl
       nimg=0
       do 150 i=nfrst,nfrend,itvl
        nimg=nimg+1
        iimg(nimg)=i
  150  continue
      else
ccc    write(6,'('' input number of frames                 : ''$)')
       read(5,*) nimg
ccc    write(6,'('' input frame numbers                    : ''$)')
       read(5,*) (iimg(i),i=1,nimg)
      end if
ccc   write(6,'('' input integration frames'')')
ccc   write(6,'(''      for calibration                   : ''$)')
      read(5,*) ical
ccc   write(6,'(''      for image processing              : ''$)')
      read(5,*) itgr
c
ccc   write(6,'('' input lower limit of criterion(r+l,r-l): ''$)')
      read(5,*) cradd,crsub
ccc   write(6,'('' input factor of criterion              : ''$)')
      read(5,*) crfac
c
***** default parameters *****
* x-/y-offsets on the dirty map
      xoff0=999.
      yoff0=999.
* main beam correction (1)yes (2)no
      mbcyn=1
* disk restoration     (1)yes (2)no
      msrstr=1
c
      loutf=index(outf,' ')-1
c
      call cft512(mapra,maprb,0)
c
      iunit=11
c
      call rdpar(iunit,infpar,
     -     dec1,dec2,dec3,ha1,ha2,ha3,
     -     solr1,solr2,solr3,solp1,solp2,solp3,solb1,solb2,solb3,
     -     nantst)
c
      call msrc(nantst,1.,gbeam)
c
      iounit=2
      nbit=-32
      naxis=2
      naxes(1)=2000
      naxes(2)=4
c
      do 10 n=1,2000
       par(n)=0.
       pal(n)=0.
   10 continue
c
      do 100 iifr=1,nimg
       iframe=iimg(iifr)
c
***** initialize rlph *****
       if (iifr.eq.1) rlph(1)=1.
c
       call snapuv17(iunit,filehd,iframe,iframe+itgr-1,
     -     nantst,ical,icald,rlph,1.,
     -     ndjst,msjst,msjsts,msjste,natt,nfalt,nfreq,nstat,
     -     flxr,flxl,mapra,maprb,mapla,maplb,nfrmst)
       if (nfrmst.ne.0) go to 100
c
       call cft512(mapra,maprb,-1)
c      do 180 il=1,512,64
c  180 write(6,'('' '',8f9.2)') (mapra(i,il),i=1,512,64)
       call cft512(mapla,maplb,-1)
c      do 190 il=1,512,64
c  190 write(6,'('' '',8f9.2)') (mapla(i,il),i=1,512,64)
c
ccc    write(6,'('' time ='',i12)') msjst
       call caleph(dec1,dec2,dec3,ha1,ha2,ha3,
     -     solr1,solr2,solr3,solp1,solp2,solp3,solb1,solb2,solb3,
     -     msjst,dec,ha,solr,solp,solb,az,alt,za,pmat)
ccc    write(6,'('' dec  ='',f10.3,''   ha   ='',f10.3,''   solr ='',
ccc  -     f10.3)') dec,ha,solr
ccc    write(6,'('' p    ='',f10.3,''   b0   ='',f10.3,''   alt  ='',
ccc  -     f10.3,''   az   ='',f10.3,''   za   ='',f10.3)')
ccc  -      solp,solb,alt,az,za
ccc    write(6,'('' matrix ='',4f10.4)') (pmat(i),i=1,4)
c       call msun17(pmat,solr*1.0125,
c     -     nantst,0.,0.,cmapr,cmapl,mstype)
c
***** preclean process for brightness normalization *****
c
      nmaxa=2000
      nmaxs=2000
      do 200 l=1,512
      do 200 i=1,512
  200 maprb(i,l)=mapra(i,l)
ccc   write(6,'('' rcp :'')')
      call clean17(3,maprb,pmat,solr*1.0125,
     -     nantst,gbeam,cradd,crfac,cra,0.,0.,corf,
     -     nmaxa,bgrr,dskr,pir,plr,pwr,par,cmapr)
c
      do 210 l=1,512
      do 210 i=1,512
  210 maprb(i,l)=mapla(i,l)
ccc   write(6,'('' lcp :'')')
      call clean17(3,maprb,pmat,solr*1.0125,
     -     nantst,gbeam,cradd,crfac,cra,0.,0.,corf,
     -     nmaxs,bgrl,dskl,pil,pll,pwl,pal,cmapl)
c
***** add and subtract rcp/lcp *****
c
      do 300 l=1,512
      do 300 i=1,512
  300 maprb(i,l)=mapra(i,l)-bgrr-(mapla(i,l)-bgrl)/dskl*dskr
      do 310 l=1,512
      do 310 i=1,512
  310 mapra(i,l)=(mapra(i,l)+mapla(i,l)/dskl*dskr)/2.
      do 320 l=1,512
      do 320 i=1,512
  320 mapla(i,l)=maprb(i,l)/2.
c
***** clean process *****
c
      if (xoff0.eq.999. .or. yoff0 .eq.999.) then
       pxew=999.
       pxns=999.
      else
       pxew=xoff0
       pxns=yoff0
      end if
      nmaxa=2000
      nmaxs=2000
c
      if (msrstr.eq.1) then
       menu=2
      else
       menu=22
      endif
      call clean17(menu,mapra,pmat,solr*1.0125,
     -     nantst,gbeam,cradd,crfac,cra,pxew,pxns,corf,
     -     nmaxa,bgr,dskbr,pir,plr,pwr,par,cmapr)
c      do 350 il=1,512,64
c  350 write(6,'('' '',8f9.2)') (cmapr(i,il),i=1,512,64)
c
      if (crfac.ne.0.) then
       crsub2=cra
      else
       crsub2=crsub
      endif
      call clean17(-1,mapla,pmat,solr*1.0125,
     -     nantst,gbeam,crsub2,0.,crs,0.,0.,0.,
     -     nmaxs,bgr,dskbr,pil,pll,pwl,pal,cmapl)
c      do 360 il=1,512,64
c  360 write(6,'('' '',8f9.2)') (cmapl(i,il),i=1,512,64)
c
       if (mbcyn.eq.1) then
        pdskbr=-dskbr
       else
        pdskbr=dskbr
       endif
       call proj17(pmat,solp,pdskbr,pxew,pxns,xoffp,yoffp,cmapr)
c      do 370 il=1,512,64
c  370 write(6,'('' '',8f9.2)') (cmapr(i,il),i=1,512,64)
       call proj17(pmat,solp,pdskbr,pxew,pxns,xoffp,yoffp,cmapl)
c      do 380 il=1,512,64
c  380 write(6,'('' '',8f9.2)') (cmapl(i,il),i=1,512,64)
c
      call projxy17(pmat,solp,dskbr,pxew,pxns,pir,plr,par,pxr,pyr,pacr)
      call projxy17(pmat,solp,dskbr,pxew,pxns,pil,pll,pal,pxl,pyl,pacl)
c
      do 390 n=1,2000
         pr(n,1)=real(pxr(n))
         pr(n,2)=real(pyr(n))
         pr(n,3)=real(pwr(n))
         pr(n,4)=pacr(n)/0.02
         pl(n,1)=real(pxl(n))
         pl(n,2)=real(pyl(n))
         pl(n,3)=real(pwl(n))
         pl(n,4)=pacl(n)/0.02
  390 continue
c
      ndut=ndjst
      msut=msjst
      call jst2ut(ndut,msut)
         write(datetime,'(i6.6,''_'',i6.6,''.'',i3.3)')
     -      mod(ndut,1000000),msut/1000,mod(msut,1000)
      if (npix .eq. 512) then
       fform='f'
      else
       fform='p'
      endif
      if (nstat.ne.3) then
       ndchar=13
      else
       ndchar=17
      endif
c
      do 400 j=1,4
      do 400 i=1,2000
  400 outmap((j-1)*2000+i)=pr(i,j)
      outfft=outf(1:loutf)//'/c'//fform//'a'//datetime(1:ndchar)
      call ftinit(iounit,outfft,2880,istat) 
      if (istat.ne.0) go to 900
      dattyp='clean_component'
      call puthdr(iounit,nbit,naxis,naxes,
     -     ndjst,msjst,msjsts,msjste,iframe,iframe+itgr-1,
     -     0,natt,nfalt,nfreq,nstat,dattyp)
c
      call putoep(iounit,solr,solp,solb,dec,ha,az,alt,za,pmat)
      call putcor(iounit,xoffa,yoffa,real(npix)/2.+0.5,
     -  real(npix)/2.+0.5,4.91106,4.91106,1)
      call putpdp(iounit,1.0125,ical,cradd,crfac,mbcyn,msrstr,
     -    pxew,pxns,corf,dskbr,cra,nmaxa)
c
      call ftpdef(iounit,nbit,naxis,naxes,0,0,istat)
      call ftppre(iounit,0,1,naxes(1)*naxes(2),outmap,istat)
      call ftclos(iounit,istat)
      if (istat.ne.0) go to 900
c     write(6,'('' file '',a40,'' created'')') outfft
      write(6,'('' out file '',a70)') outfft
c
      do 410 j=1,4
      do 410 i=1,2000
  410 outmap((j-1)*2000+i)=pl(i,j)
      outfft=outf(1:loutf)//'/c'//fform//'s'//datetime(1:ndchar)
      call ftinit(iounit,outfft,2880,istat)
      if (istat.ne.0) go to 900
      dattyp='clean_component'
      call puthdr(iounit,nbit,naxis,naxes,
     -     ndjst,msjst,msjsts,msjste,iframe,iframe+itgr-1,
     -     2,natt,nfalt,nfreq,nstat,dattyp)
c
      call putoep(iounit,solr,solp,solb,dec,ha,az,alt,za,pmat)
      call putcor(iounit,xoffa,yoffa,real(npix)/2.+0.5,
     -  real(npix)/2.+0.5,4.91106,4.91106,1)
      call putpdp(iounit,1.0125,ical,crsub,crfac,mbcyn,msrstr,
     -    pxew,pxns,corf,dskbr,crs,nmaxs)
c
      call ftpdef(iounit,nbit,naxis,naxes,0,0,istat)
      call ftppre(iounit,0,1,naxes(1)*naxes(2),outmap,istat)
      call ftclos(iounit,istat)
      if (istat.ne.0) go to 900
ccc   write(6,'('' file '',a40,'' created'')') outfft
      write(6,'('' out file '',a70)') outfft
c
  100 continue
c
  900 write(6,'('' istat ='',i8)') istat
c
      end
c
c
***********************************************************************
*     putpdp - write promgram dependent parameters to fits file       *
*       parameters of CLEAN procedure                                 *
*       iunit  : input,integer fortran unit number for fits file      *
*       solrf  : input,real      factor of solar radius correction    *
*       nfrcal : input,integer   number of calibration frames         *
*       crin   : input,real      (lower limit of) clean criterion     *
*       crfac  : input,real      factor for clean criterion           *
*       mbcyn  : input,integer   main beam correction 1=yes, 2=no     *
*       msrstr : input,integer   disk restoration 1=yes, 2=no         *
*       pxew   : input,real      x-offset of the center of dirty disk *
*       pxns   : input,real      y-offset of the center of dirty disk *
*       corf   : input,real      correlation of dirty disk and model  *
*       dskbr  : input,real      brightness on the dirty disk         *
*       cr     : input,real      actual clean criterion               *
*       ncmp   : input,integer   number of clean components           *
*                                       April 19, 1999  Y. Hanaoka    *
***********************************************************************
c
      subroutine putpdp(iunit,solrf,nfrcal,crin,crfac,mbcyn,msrstr,
     -    pxew,pxns,corf,dskbr,cr,ncmp)
c
      character*44 cmnt
c
      cmnt=' '
c
      call ftpkys(iunit,'progname',
     -     'clean_component17 v6.2  H. Koshiishi',cmnt,istat)
c
      call ftpkys(iunit,'bunit','K','disk = 10000 K',istat)
      call ftpkyf(iunit,'solr-fac',solrf,5,'radius correction factor',
     -     istat)
      call ftpkyj(iunit,'nfrcal',nfrcal,'number of calibration frames',
     -     istat)
      call ftpkyf(iunit,'crinput',crin,2,
     -     'lower limit of clean criterion',istat)
      call ftpkyf(iunit,'crfactor',crfac,2,
     -     'factor clean criterion',istat)
      if (mbcyn.eq.1) then
       call ftpkys(iunit,'mbeamc','yes','main beam correction',istat)
      else
       call ftpkys(iunit,'mbeamc','no','main beam correction',istat)
      end if
      if (msrstr.eq.1) then
       call ftpkys(iunit,'diskrstr','yes','disk restoration',istat)
      else
       call ftpkys(iunit,'diskrstr','no','disk restoration',istat)
      end if
c
      call ftpkyf(iunit,'ddoff1',pxew,2,
     -     'x-offset of the dirty disk',istat)
      call ftpkyf(iunit,'ddoff2',pxns,2,
     -     'y-offset of the dirty disk',istat)
      call ftpkyf(iunit,'ddcorr',corf,4,
     -     'correlation between dirty disk and model',istat)
      call ftpkyf(iunit,'dskbr',dskbr,2,
     -     'brightness of the dirty disk',istat)
      call ftpkyf(iunit,'criter',cr,2,'actual clean criterion',
     -     istat)
      call ftpkyj(iunit,'ncompo',ncmp,
     -     'number of clean components',istat)
c
c
      if (istat.ne.0) write(6,'('' putpdp : istat ='',i6)') istat
      return
      end
c
c
***********************************************************************
*     projxy17 - projection correction for clean components 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    disk brightness of the snapshot map     *
*       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                                 *
*       ix    : input,integer(2000) x-position of clean component     *
*       iy    : input,integer(2000) y-position of clean component     *
*       a     : input,real(2000) amplitude of clean component         *
*       ipx   : output,integer(2000) x-position of                    *
*                 projected clean component                           *
*       ipy   : output,integer(2000) y-position of                    *
*                 projected clean component                           *
*       ca    : output,real(2000) corrected amplitude of              *
*                 clean component                                     *
*                                   October 10, 2007  H. Koshiishi    *
***********************************************************************
c
      subroutine projxy17(pmat,solp,dskbr,pxew,pxns,ix,iy,a,ipx,ipy,ca)
c
      integer ix(2000),iy(2000),ipx(2000),ipy(2000)
      real a(2000),ca(2000),pmat(4)
c
      rd=57.29577951
      pi=3.141592653
c
      adskbr=abs(dskbr)
c
*****
c
      do 510 n=1,2000
c
      if (abs(a(n)).ne.0.) then
c
      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)
c
       rx=amod(xx+pxew+256.5+1024.,512.)+0.5
       ry=amod(yy+pxns+256.5+1024.,512.)+0.5
       ii=int(rx+0.5)
       ll=int(ry+0.5)
c
       if (ii.eq.ix(n) .and. ll.eq.iy(n)) then
        ipx(n)=i
        ipy(n)=l
       endif
c
  500 continue
c
* main beam correction factor *
       if (dskbr.lt.0.) then
        rpix2=((256.5-real(ipx(n)))*(256.5-real(ipx(n)))+
     -        (256.5-real(ipy(n)))*(256.5-real(ipy(n))))*
     -        (4.91104/5556.5)*(4.91104/5556.5)
        g=exp(-3.141592*rpix2)
       else
        g=1.
       end if
       ca(n)=a(n)/g
      else
       ca(n)=0.
      end if
c
      ca(n)=ca(n)/adskbr*10000.
c
  510 continue
c
      return
      end