program snap2d17
***********************************************************************
*     create cleaned maps (512x512 etc.)                              *
*                               ver 6.2 April 19, 1999  Y. Hanaoka    *
***********************************************************************
c
      integer naxes(2),nantst(84),iimg(1000),pir(2000),pil(2000),
     -     plr(2000),pll(2000),pwr(2000),pwl(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(262144),pmat(4),
     -     par(2000),pal(2000)
      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)=npix
      naxes(2)=npix
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,pir,plr,pwr,par,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
      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,npix
      do 400 i=1,npix
  400 outmap((j-1)*npix+i)=cmapr((512-npix)/2+i,(512-npix)/2+j)
      outfft=outf(1:loutf)//'/i'//fform//'a'//datetime(1:ndchar)
      call ftinit(iounit,outfft,2880,istat) 
      if (istat.ne.0) go to 900
      dattyp='cleaned_map'
      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,npix
      do 410 i=1,npix
  410 outmap((j-1)*npix+i)=cmapl((512-npix)/2+i,(512-npix)/2+j)
      outfft=outf(1:loutf)//'/i'//fform//'s'//datetime(1:ndchar)
      call ftinit(iounit,outfft,2880,istat)
      if (istat.ne.0) go to 900
      dattyp='cleaned_map'
      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','snap2d17 v6.2  Y. Hanaoka',
     -     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