program burst17_3
***********************************************************************
*     create cleaned maps (512x512 etc.)                              *
*                                    October 1, 2008  H. Koshiishi    *
***********************************************************************
c
      integer naxes(2),nantst(84),icimg(1000),iimg(1000),
     -     pir(32767),pil(32767),plr(32767),pll(32767),
     -     pwr(32767),pwl(32767)
      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(32767),pal(32767)
      real rlph(84)
      character*80 filehd,infpar,calf,infcal,outf,outfft
      character*18 dattyp
      character*17 datetime
      character*10 datecal
      character*12 timecal
      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 first frame number               : ''$)')
      read(5,*) ncfrst
      if (ncfrst.ne.0) then
ccc    write(6,'('' input last frame number                : ''$)')
       read(5,*) ncfrend
ccc    write(6,'('' input frame interval                   : ''$)')
       read(5,*) ictvl
       ncimg=0
       do 150 i=ncfrst,ncfrend,ictvl
        ncimg=ncimg+1
        icimg(ncimg)=i
  150  continue
      else
ccc    write(6,'('' input number of frames                 : ''$)')
       read(5,*) ncimg
ccc    write(6,'('' input frame numbers                    : ''$)')
       read(5,*) (icimg(i),i=1,ncimg)
      end if
ccc   write(6,'('' input integration frames'')')
ccc   write(6,'(''      for calibration                   : ''$)')
      read(5,*) iccal
ccc   write(6,'(''      for image processing              : ''$)')
      read(5,*) ictgr
c
ccc   write(6,'('' quiet : 1 / burst : 2               : ''$)')
      read(5,*) iqb
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 160 i=nfrst,nfrend,itvl
        nimg=nimg+1
        iimg(nimg)=i
  160  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 number limit of criterion(r+l,r-l): ''$)')
      read(5,*) inmaxa,inmaxs
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 *****
* main beam correction (1)yes (2)no
      mbcyn=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
      lfilehd=index(filehd,' ')-1      
      if (ncfrst.lt.10)
     -   write(calf,'(''ifd'',a6,''_0000'',i1)')
     -   filehd(lfilehd-5:lfilehd),ncfrst
      if (ncfrst.ge.10 .and. ncfrst.lt.100)
     -   write(calf,'(''ifd'',a6,''_000'',i2)')
     -   filehd(lfilehd-5:lfilehd),ncfrst
      if (ncfrst.ge.100 .and. ncfrst.lt.1000)
     -   write(calf,'(''ifd'',a6,''_00'',i3)')
     -   filehd(lfilehd-5:lfilehd),ncfrst
      if (ncfrst.ge.1000 .and. ncfrst.lt.10000)
     -   write(calf,'(''ifd'',a6,''_0'',i4)')
     -   filehd(lfilehd-5:lfilehd),ncfrst
      if (ncfrst.ge.10000)
     -   write(calf,'(''ifd'',a6,''_'',i5)')
     -   filehd(lfilehd-5:lfilehd),ncfrst
c
      infcal=outf(1:loutf)//'/'//calf
c
      call rdcal(iunit,infcal,datecal,timecal,nfrcal,nfrtgr,icald,
     -     cpxew,cpxns,ccorf)
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
c
      iframe=iimg(iifr)
c
***** initialize rlph *****
c
      if (iifr.eq.1) rlph(1)=1.
c
      call snapuv17(iunit,filehd,iframe,iframe+itgr-1,
     -     nantst,0,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
      call cft512(mapla,maplb,-1)
c
      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)
c
      if (iqb.eq.1) then
c
***** preclean process for brightness normalization *****
c
      nmaxa=inmaxa
      nmaxs=inmaxs
c
      do 200 l=1,512
      do 200 i=1,512
  200 maprb(i,l)=mapra(i,l)
c
      call clean17(3,maprb,pmat,solr*1.0125,
     -     nantst,gbeam,cradd,crfac,cra,0.,0.,0.,
     -     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)
c
      call clean17(3,maprb,pmat,solr*1.0125,
     -     nantst,gbeam,cradd,crfac,cra,0.,0.,0.,
     -     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
      endif
c
***** clean process *****
c
      if (iqb.eq.1) then
c
      pxew=999.
      pxns=999.
      nmaxa=inmaxa
      nmaxs=inmaxs
c
      call clean17(2,mapra,pmat,solr*1.0125,
     -     nantst,gbeam,cradd,crfac,cra,pxew,pxns,corf,
     -     nmaxa,bgr,dskbr,pir,plr,pwr,par,cmapr)
c
      if (crfac.ne.0.) then
       crsub2=cra
      else
       crsub2=crsub
      endif
c
      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
      endif
      if (iqb.eq.2) then
c
      pxew=0.
      pxns=0.
      dskbr=0.
      nmaxa=inmaxa
      nmaxs=inmaxs
      
c
      call clean17(1,mapra,pmat,solr*1.0125,
     -     nantst,gbeam,cradd,crfac,cra,pxew,pxns,corf,
     -     nmaxa,bgr,dskbr,pir,plr,pwr,par,cmapr)
c
      call clean17(1,mapla,pmat,solr*1.0125,
     -     nantst,gbeam,cradd,crfac,crs,pxew,pxns,corf,
     -     nmaxs,bgr,dskbr,pil,pll,pwl,pal,cmapl)
c
      endif
c
      if (iqb.eq.1) pdskbr=-dskbr
      if (iqb.eq.2) pdskbr=-1./10000.
c
      call proj17(pmat,solp,pdskbr,cpxew,cpxns,xoffp,yoffp,cmapr)
c
      call proj17(pmat,solp,pdskbr,cpxew,cpxns,xoffp,yoffp,cmapl)
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)
      ndchar=13
c
      if (npix .eq. 512) then
       fform='f'
      else
       fform='p'
      endif
c
      if (iqb.eq.1) then 
c
      msrstr=1
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)//'/m'//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 putpdf(iounit,iqb,1.0125,0,cradd,crfac,mbcyn,msrstr,
     -    pxew,pxns,corf,dskbr,cra,nmaxa)
      call putcon(iounit,iqb)
      call putcpr(iounit,datecal,timecal,nfrcal,nfrtgr,
     -    cpxew,cpxns,ccorf)
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)
c
      if (istat.ne.0) go to 900
      write(6,'('' out file  :'',a40)') 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)//'/m'//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 putpdf(iounit,iqb,1.0125,0,crsub,crfac,mbcyn,msrstr,
     -    pxew,pxns,corf,dskbr,crs,nmaxs)
      call putcon(iounit,iqb)
      call putcpr(iounit,datecal,timecal,nfrcal,nfrtgr,
     -    cpxew,cpxns,ccorf)
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)
c
      if (istat.ne.0) go to 900
      write(6,'('' out file  :'',a40)') outfft
c
      endif
      if (iqb.eq.2) then
c
      msrstr=2
c
      do 500 j=1,npix
      do 500 i=1,npix
  500 outmap((j-1)*npix+i)=cmapr((512-npix)/2+i,(512-npix)/2+j)
      outfft=outf(1:loutf)//'/m'//fform//'r'//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,
     -     1,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 putpdf(iounit,iqb,1.0125,0,cradd,crfac,mbcyn,msrstr,
     -    pxew,pxns,corf,dskbr,cra,nmaxa)
      call putcon(iounit,iqb)
      call putcpr(iounit,datecal,timecal,nfrcal,nfrtgr,
     -    cpxew,cpxns,ccorf)
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)
c
      if (istat.ne.0) go to 900
      write(6,'('' out file  :'',a40)') outfft
c
      do 510 j=1,npix
      do 510 i=1,npix
  510 outmap((j-1)*npix+i)=cmapl((512-npix)/2+i,(512-npix)/2+j)
      outfft=outf(1:loutf)//'/m'//fform//'l'//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,
     -     -1,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 putpdf(iounit,iqb,1.0125,0,cradd,crfac,mbcyn,msrstr,
     -    pxew,pxns,corf,dskbr,crs,nmaxs)
      call putcon(iounit,iqb)
      call putcpr(iounit,datecal,timecal,nfrcal,nfrtgr,
     -    cpxew,cpxns,ccorf)
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)
c
      if (istat.ne.0) go to 900
      write(6,'('' out file  :'',a40)') outfft
c
      endif
c
  100 continue
c
  900 write(6,'('' istat ='',i8)') istat
c
      end
c
c
***********************************************************************
*     putpdf - write promgram dependent parameters to fits file       *
*       parameters of CLEAN procedure                                 *
*       iunit  : input,integer fortran unit number for fits file      *
*       iqb    : input,integer procedure for quiet(1) or burst(2)     *
*       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           *
*                                    February 14, 2007  H. Koshiishi  *
*                                       April 19, 1999  Y. Hanaoka    *
***********************************************************************
c
      subroutine putpdf(iunit,iqb,solrf,nfrcal,crin,crfac,mbcyn,msrstr,
     -                  pxew,pxns,corf,dskbr,cr,ncmp)
c
      character*44 cmnt
c
      cmnt=' '
c
      call ftpkys(iunit,'progname','burst17_1 H.Koshiishi',
     -     cmnt,istat)
c
      if (iqb.eq.1) then
       call ftpkys(iunit,'bunit','K','disk = 10000 K',istat)
      else
       call ftpkys(iunit,'bunit','K','disk = ***** K',istat)
      end if
      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,'('' putpdf : istat ='',i6)') istat
      return
      end
c
c
***********************************************************************
*     putcon - write procedure condition to fits file                 *
*       iunit  : input,integer fortran unit number for fits file      *
*       iqb    : input,integer procedure for quiet(1) or flare(2)     *
*                                   February 5, 2007  H. Koshiishi    *
***********************************************************************
c
      subroutine putcon(iunit,iqb)
c
      integer iqb
      character*44 cmnt
c
      cmnt=' '
c
      if (iqb.eq.1) call ftpkyj(iunit,'iqb',iqb,
     -                   'procedure for quiet condition',istat)
      if (iqb.eq.2) call ftpkyj(iunit,'iqb',iqb,
     -                   'procedure for flare condition',istat)
c
      if (istat.ne.0) write(6,'('' putcon : istat ='',i6)') istat
c
      return
      end
c
c
***********************************************************************
*     putcpr - write calibration parameters to fits file              *
*       iunit  : input,integer fortran unit number for fits file      *
*       datecal: input,character*10 calibration date in ut            *
*       timecal: input,character*12 center of calibration time in ut  *
*       nfrcal : input,integer number of calibration frames           *
*       nfrtgr : input,integer number of integration frames           *
*       cpxew  : input,real x-offset in calibration                   *
*       cpxns  : input,real y-offset in calibration                   *
*       ccorf  : input,real correlation in calibration                *
*                                   February 5, 2007  H. Koshiishi    *
***********************************************************************
c
      subroutine putcpr(iunit,datecal,timecal,nfrcal,nfrtgr,
     -                  cpxew,cpxns,ccorf)
c
      integer nfrcal,nfrtgr
      real cpxew,cpxns,ccorf
      character*10 datecal
      character*12 timecal
      character*44 cmnt
c
      cmnt=' '
c
      call ftpkys(iunit,'date-cal',datecal,cmnt,istat)
      call ftpkys(iunit,'time-cal',timecal,cmnt,istat)
c
      call ftpkyj(iunit,'cnfrcal',nfrcal,
     -     'calibration frames in calibration',istat)
      call ftpkyj(iunit,'cnfrtgr',nfrtgr,
     -     'integration frames in calibration',istat)
c
      call ftpkyf(iunit,'cddoff1',cpxew,2,
     -     'x-offset in calibration',istat)
      call ftpkyf(iunit,'cddoff2',cpxns,2,
     -     'y-offset in calibration',istat)
      call ftpkyf(iunit,'cddcorr',ccorf,4,
     -     'correlation in calibration',istat)
c
      if (istat.ne.0) write(6,'('' putcpr : istat ='',i6)') istat
c
      return
      end
c
c
***********************************************************************
*     rdcal - read calbiration file                                   *
*       iunit  : input,integer fortran unit number for fits file      *
*       infcal : input,charactor*80 calibration filename              *
*       datecal: output,character*10 calibration date                 *
*       timecal: output,character*12 calibration time                 *
*       nfrcal : output,integer number of calibration frames          *
*       nfrtgr : output,integer number of integration frames          *
*       icald  : output,integer*2(336) calibration data               *
*       pxew   : output,real x-offset of the center of dirty disk     *
*       pxns   : output,real y-offset of the center of dirty disk     *
*       corf   : output,real correlation of dirty disk and model      *
*                                   January 30, 2007  H. Koshiishi    *
***********************************************************************
c
      subroutine rdcal(iunit,infcal,datecal,timecal,
     -                 nfrcal,nfrtgr,icald,pxew,pxns,corf)
c     
      integer nfrcal,nfrtgr
      integer data(336)
      integer*2 icald(336)
      real pxew,pxns,corf
      character*80 infcal
      character*10 datecal
      character*12 timecal
      character*44 cmnt
      logical anyf
c
      call ftopen(iunit,infcal,0,ibsize,istat)
      call ftgkys(iunit,'date-obs',datecal,cmnt,istat)
      call ftgkys(iunit,'time-obs',timecal,cmnt,istat)
      call ftgkyj(iunit,'nfrcal',nfrcal,cmnt,istat)
      call ftgkyj(iunit,'nfrtgr',nfrtgr,cmnt,istat)
      call ftgkye(iunit,'ddoff1',pxew,cmnt,istat)
      call ftgkye(iunit,'ddoff2',pxns,cmnt,istat)
      call ftgkye(iunit,'ddcorr',corf,cmnt,istat)
      call ftgpvj(iunit,0,1,336,0,data,anyf,istat)
      call ftclos(iunit,istat)
c
      if (istat.eq.0) then
         write(6,'('' open file :'',a15)') infcal
         write(6,*) 'pxew , pxns =',pxew,',',pxns
      endif
      if (istat.ne.0) write(6,'('' rdcal : istat ='',i6)') istat
c
      do 100 i=1,336
100   icald(i)=data(i)
c
      return
      end