program burst34_2
***********************************************************************
*     create cleaned maps (1024x1024 etc.)                            *
*                                     March 1, 2007  H. Koshiishi     *
***********************************************************************
c
      integer naxes(2),nantst(84),icimg(1000),iimg(1000),
     -     pi(32767),pl(32767),pw(32767)
      integer*2 icald(168)
      real gbeam(512,512,0:10),flx(84),
     -     mapa(512,512),mapb(512,512),
     -     cmap(1024,1024),outmap(1048576),pmat(4),pa(32767)
      character*80 filehd,infpar,infref,outf,outfft
      character*18 dattyp
      character*12 dateframe
      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 reference filename               : ''$)')
      read(5,'(a)') infref
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/2.45553
       yoffp=yoffa/2.45553
      else
       xoffp=xoffa
       yoffp=yoffa
       xoffa=xoffp*2.45553
       yoffa=yoffp*2.45553
      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        : ''$)')
      read(5,*) inmaxa
ccc   write(6,'('' input lower limit of criterion         : ''$)')
      read(5,*) cradd
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(mapa,mapb,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
c
      iframe=iimg(iifr)
c
      call snapuv34(iunit,filehd,iframe,iframe+itgr-1,
     -     nantst,ical,icald,1.,
     -     ndjst,msjst,msjsts,msjste,natt,nfalt,nfreq,nstat,
     -     flx,mapa,mapb,nfrmst)
      if (nfrmst.ne.0) go to 100
c
      call cft512(mapa,mapb,-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
***** clean process *****
c
      if (iqb.eq.1) then
c
      pxew=999.
      pxns=999.
      nmaxa=inmaxa
c
      call clean34(2,mapa,pmat,solr*1.0125,
     -     nantst,gbeam,cradd,crfac,cra,pxew,pxns,corf,
     -     nmaxa,bgr,dskbr,pi,pl,pw,pa,cmap)
c
      endif
      if (iqb.eq.2) then
c
      pxew=0.
      pxns=0.
      dskbr=0.
      nmaxa=inmaxa
c
      call clean34(1,mapa,pmat,solr*1.0125,
     -     nantst,gbeam,cradd,crfac,cra,pxew,pxns,corf,
     -     nmaxa,bgr,dskbr,pi,pl,pw,pa,cmap)
      endif
c
      lfilehd=index(filehd,' ')-1
      if (iframe.lt.10)
     -    write(dateframe,'(a6,''_0000'',i1)')
     -    filehd(lfilehd-5:lfilehd),iframe
      if (iframe.ge.10 .and. iframe.lt.100)
     -    write(dateframe,'(a6,''_000'',i2)')
     -    filehd(lfilehd-5:lfilehd),iframe
      if (iframe.ge.100 .and. iframe.lt.1000)
     -    write(dateframe,'(a6,''_00'',i3)')
     -    filehd(lfilehd-5:lfilehd),iframe
      if (iframe.ge.1000 .and. iframe.lt.10000)
     -    write(dateframe,'(a6,''_0'',i4)')
     -    filehd(lfilehd-5:lfilehd),iframe
      if (iframe.ge.10000)
     -    write(dateframe,'(a6,''_'',i5)')
     -    filehd(lfilehd-5:lfilehd),iframe
c
      if (npix .eq. 1024) then
       fform='f'
      else
       fform='p'
      endif
c
      if (iqb.eq.1) msrstr=1
      if (iqb.eq.2) msrstr=2
c
      do 400 j=1,npix
      do 400 i=1,npix
  400 outmap((j-1)*npix+i)=cmap((1024-npix)/2+i,(1024-npix)/2+j)
      outfft=outf(1:loutf)//'/l'//fform//'z'//dateframe
      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,0.,0.,real(npix)/2.+0.5,
     -  real(npix)/2.+0.5,2.45553,2.45553,1)
      call putpdf(iounit,iqb,1.0125,ical,cradd,crfac,mbcyn,msrstr,
     -    pxew,pxns,corf,dskbr,cra,nmaxa)
      call putcon(iounit,iqb)
      call putprj(iounit,pmat,solp,dskbr)
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
      if (iqb.eq.1) pdskbr=-dskbr
      if (iqb.eq.2) pdskbr=-1./10000.
c
      call proj34(pmat,solp,pdskbr,0.,0.,0.,0.,cmap)
c
      do 500 j=1,npix
      do 500 i=1,npix
  500 outmap((j-1)*npix+i)=cmap((1024-npix)/2+i,(1024-npix)/2+j)
      outfft=outf(1:loutf)//'/k'//fform//'z'//dateframe
      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,0.,0.,real(npix)/2.+0.5,
     -  real(npix)/2.+0.5,2.45553,2.45553,1)
      call putpdf(iounit,iqb,1.0125,ical,cradd,crfac,mbcyn,msrstr,
     -    pxew,pxns,corf,dskbr,cra,nmaxa)
      call putcon(iounit,iqb)
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
  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_2 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 burst(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 burst condition',istat)
c
      if (istat.ne.0) write(6,'('' putcon : istat ='',i6)') istat
c
      return
      end
c
c
***********************************************************************
*     putprj - write projection parameter to fits file                *
*       iunit  : input,integer fortran unit number for fits file      *
*       pmat   : input,real(4) projection matrix                      *
*       solp   : input,real    polar angle of the sun (degree)        *
*       dskbr  : input,real    disk brightness                        *
*                                      March 5, 2007  H. Koshiishi    *
***********************************************************************
c
      subroutine putprj(iunit,pmat,solp,dskbr)
c
      real solp,dskbr
      real pmat(4)
      character*44 cmnt
c
      cmnt=' '
c
      call ftpkye(iunit,'dsolp',solp,8,
     -     'solar polar angle for projection',istat)
      call ftpkye(iunit,'dpmat1',pmat(1),8,
     -     'projection matrix for projection',istat)
      call ftpkye(iunit,'dpmat2',pmat(2),8,
     -     'projection matrix for projection',istat)
      call ftpkye(iunit,'dpmat3',pmat(3),8,
     -     'projection matrix for projection',istat)
      call ftpkye(iunit,'dpmat4',pmat(4),8,
     -     'projection matrix for projection',istat)
      call ftpkye(iunit,'ddskbr',dskbr,8,
     -     'brightness of the dirty disk for projection',istat)
c
      if (istat.ne.0) write(6,'('' putprj : istat ='',i6)') istat
c
      return
      end