program cor34
***********************************************************************
*     create correlation maps (1024x1024 etc.)                        *
*                                   March 27, 2007  H. Koshiishi      *
***********************************************************************
c
      integer naxes(2),icimg(1000),iimg(1000)
      real map1(1024,1024),map2(1024,1024),map3(1024,1024),
     -     cmap(1024,1024),inmap(1048576),outmap(1048576)
      real map171(512,512),map172(1024,1024),
     -     inmap17(262144)
      real cormap1(41,41),cormap2(5,5),
     -     xarr(5),yarr(5),fmap(5,5),a(9),wk(135)
      real pmat(4)
      character*80 filehd,infpar,infref,infft,outf,outfft
      character*12 dateframe
      character*13 datetime
      character*10 datecal
      character*12 timecal
      character*18 dtobs,tmobs,strobs,endobs,jstdt,jsttm,jsts,jste,
     -             natt,nfalt,nfreq,nfrstat,dattyp
      character*44 cmnt
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 *****
c
      mbcyn=1
      msrstr=0
c
      loutf=index(outf,' ')-1
c
      iunit=11
c
      iounit=2
      nbit=-32
      naxis=2
      npix=512
      naxes(1)=npix
      naxes(2)=npix
c
      infft=infref
      call ftopen(iounit,infft,0,ibsize,istat)
      call ftgpve(iounit,0,1,naxes(1)*naxes(2),0,inmap17,anyf,istat)
      call ftclos(iounit,istat)
      if (istat.eq.0) write(6,'('' open file :'',a16)') infft
      if (istat.ne.0) write(6,'('' cor34 : istat ='',i8)') istat
      do 10 j=1,npix
      do 10 i=1,npix
   10 map171((512-npix)/2+i,(512-npix)/2+j)=inmap17((j-1)*npix+i)
c
      do 20 i=1,512
      do 20 j=1,512
c
         map172(i*2-1,j*2-1)=map171(i,j)
         map172(i*2,j*2-1)=map171(i,j)
         map172(i*2-1,j*2)=map171(i,j)
         map172(i*2,j*2)=map171(i,j)
c
   20 continue
c
      npix=1024
      naxes(1)=npix
      naxes(2)=npix
c
      do 100 iifr=1,nimg
c
      iframe=iimg(iifr)
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
      infft=outf(1:loutf)//'/jfz'//dateframe
      call ftopen(iounit,infft,0,ibsize,istat)
c
      call ftgkye(iounit,'crval1',xoffa,cmnt,istat)
      call ftgkye(iounit,'crval2',yoffa,cmnt,istat)
c
      call ftgkys(iounit,'date-cal',datecal,cmnt,istat)
      call ftgkys(iounit,'time-cal',timecal,cmnt,istat)
      call ftgkyj(iounit,'cnfrcal',iccal,cmnt,istat)
      call ftgkyj(iounit,'cnfrtgr',ictgr,cmnt,istat)
      call ftgkye(iounit,'cddoff1',cpxew,cmnt,istat)
      call ftgkye(iounit,'cddoff2',cpxns,cmnt,istat)
      call ftgkye(iounit,'cddcorr',ccorf,cmnt,istat)
c
      call ftgpve(iounit,0,1,naxes(1)*naxes(2),0,inmap,anyf,istat)
      call ftclos(iounit,istat)
      if (istat.eq.0) write(6,'('' open file :'',a15)') infft
      if (istat.ne.0) write(6,'('' cor34 : istat ='',i8)') istat
      do 200 j=1,npix
      do 200 i=1,npix
  200 map1((1024-npix)/2+i,(1024-npix)/2+j)=inmap((j-1)*npix+i)
c
      infft=outf(1:loutf)//'/kfz'//dateframe
      call ftopen(iounit,infft,0,ibsize,istat)
      call ftgpve(iounit,0,1,naxes(1)*naxes(2),0,inmap,anyf,istat)
      call ftclos(iounit,istat)
      if (istat.eq.0) write(6,'('' open file :'',a15)') infft
      if (istat.ne.0) write(6,'('' cor34 : istat ='',i8)') istat
      do 210 j=1,npix
      do 210 i=1,npix
  210 map2((1024-npix)/2+i,(1024-npix)/2+j)=inmap((j-1)*npix+i)
c
      infft=outf(1:loutf)//'/lfz'//dateframe
      call ftopen(iounit,infft,0,ibsize,istat)
c
      call ftgkys(iounit,'date-obs',dtobs,cmnt,istat)
      call ftgkys(iounit,'time-obs',tmobs,cmnt,istat)
      call ftgkys(iounit,'strt-obs',strobs,cmnt,istat)
      call ftgkys(iounit,'end-obs',endobs,cmnt,istat)
      call ftgkys(iounit,'jstdate',jstdt,cmnt,istat)
      call ftgkys(iounit,'jsttime',jsttm,cmnt,istat)
      call ftgkys(iounit,'jst-strt',jsts,cmnt,istat)
      call ftgkys(iounit,'jst-end',jste,cmnt,istat)
      call ftgkyj(iounit,'startfrm',nfrst,cmnt,istat)
      call ftgkyj(iounit,'endfrm',nfrend,cmnt,istat)
      call ftgkys(iounit,'att-10db',natt,cmnt,istat)
      call ftgkys(iounit,'obs-mode',nfalt,cmnt,istat)
      call ftgkys(iounit,'obs-freq',nfreq,cmnt,istat)
      call ftgkys(iounit,'frm-stat',nfrstat,cmnt,istat)
      call ftgkys(iounit,'data-typ',dattyp,cmnt,istat)
c
      call ftgkye(iounit,'solr',solr,cmnt,istat)
      call ftgkye(iounit,'solb',solb,cmnt,istat)
      call ftgkye(iounit,'dec',dec,cmnt,istat)
      call ftgkye(iounit,'houra',ha,cmnt,istat)
      call ftgkye(iounit,'azimuth',az,cmnt,istat)
      call ftgkye(iounit,'altitude',alt,cmnt,istat)
      call ftgkye(iounit,'zangle',za,cmnt,istat)
c
      call ftgkyj(iounit,'nfrcal',ical,cmnt,istat)
      call ftgkye(iounit,'crinput',cradd,cmnt,istat)
      call ftgkye(iounit,'crfactor',crfac,cmnt,istat)
      call ftgkye(iounit,'ddoff1',dpxew,cmnt,istat)
      call ftgkye(iounit,'ddoff2',dpxns,cmnt,istat)
      call ftgkye(iounit,'ddcorr',dcorf,cmnt,istat)
      call ftgkye(iounit,'criter',cra,cmnt,istat)
      call ftgkyj(iounit,'ncompo',nmaxa,cmnt,istat)
c
      call ftgkyj(iounit,'iqb',iqb,cmnt,istat)
c
      call ftgkye(iounit,'dsolp',solp,cmnt,istat)
      call ftgkye(iounit,'dpmat1',dummy,cmnt,istat)
      pmat(1)=dummy
      call ftgkye(iounit,'dpmat2',dummy,cmnt,istat)
      pmat(2)=dummy
      call ftgkye(iounit,'dpmat3',dummy,cmnt,istat)
      pmat(3)=dummy
      call ftgkye(iounit,'dpmat4',dummy,cmnt,istat)
      pmat(4)=dummy
      call ftgkye(iounit,'ddskbr',dskbr,cmnt,istat)
c
      call ftgpve(iounit,0,1,naxes(1)*naxes(2),0,inmap,anyf,istat)
      call ftclos(iounit,istat)
      if (istat.eq.0) write(6,'('' open file :'',a15)') infft
      if (istat.ne.0) write(6,'('' cor34 : istat ='',i8)') istat
      do 220 j=1,npix
      do 220 i=1,npix
  220 map3((1024-npix)/2+i,(1024-npix)/2+j)=inmap((j-1)*npix+i)
c
      xoffp=xoffa/2.45553
      yoffp=yoffa/2.45553
c
      if (iqb.eq.1) pdskbr=-dskbr
      if (iqb.eq.2) pdskbr=-1./10000.
c
      call mean1024(map172,imx1,lmx1)
      call mean1024(map2,imx2,lmx2)
c
      call contour_trim4(map2,0.1,imx2,lmx2,ipix)
      ipix=int(sqrt(real(ipix)))
      if (ipix.gt.100) ipix=100
      if (ipix.lt.10) ipix=10
      write(6,*) 'ipix =',ipix
c
      write(6,*) 'shift1 =',-1*(imx2-imx1),',',-1*(lmx2-lmx1)
c
      do 310 ix=-20,20
      do 310 iy=-20,20
c
      do 300 i=1,1024
      do 300 j=1,1024
c
         cmap(i,j)=map2(i,j)
c
  300 continue
c
      call shift1024(imx2-imx1+ix,lmx2-lmx1+iy,cmap)
      call calcor1024(map1,cmap,imx1,lmx1,ipix,corf)
      cormap1(ix+21,iy+21)=corf
c
  310 continue
c
      xmx0=imx2-imx1
      ymx0=lmx2-lmx1
c
      call smaxa1(cormap1,ixmx,iymx,corfmx)
      write(6,*) 'shift2 =',
     -           -1*(ixmx-21),',',-1*(iymx-21),',',corfmx
c
      xmx1=real(ixmx)-21.
      ymx1=real(iymx)-21.
c
      do 330 ix=-2,2
      do 330 iy=-2,2
c
      do 320 i=1,1024
      do 320 j=1,1024
c
         cmap(i,j)=map3(i,j)
c
  320 continue
c
      call proj34(pmat,solp,pdskbr,0.,0.,
     -            xmx0+xmx1+real(ix),ymx0+ymx1+real(iy),
     -            cmap)
      call calcor1024(map1,cmap,
     -                imx1+int(xmx1),lmx1+int(ymx1),ipix,corf)
      cormap2(ix+3,iy+3)=corf
c
  330 continue
c
      do 340 i=1,5
         xarr(i)=real(i)-3.
         yarr(i)=real(i)-3.
  340 continue
c
      call LRNGAPL(xarr,5,yarr,5,cormap2,2,2,a,fmap,wk,ierr)
      call smaxa2(a,corfmx,xmx2,ymx2)
      write(6,*) 'shift3 =',
     -           -1.*xmx2,',',-1.*ymx2,',',corfmx
c
      spxew=xmx0+xmx1+xmx2
      spxns=ymx0+ymx1+ymx2
c
      call proj34(pmat,solp,pdskbr,0.,0.,spxew,spxns,map3)
      call calcor1024(map1,map3,
     -                imx1+int(xmx1),lmx1+int(ymx1),ipix,corfmxz)
      write(6,*) 'shift =',-1.*spxew,',',-1.*spxns,',',corfmxz
c
      write(datetime,'(a2,a2,a2,''_'',a2,a2,a2)')
     -      dtobs(3:4),dtobs(6:7),dtobs(9:10),
     -      tmobs(1:2),tmobs(4:5),tmobs(7:8)
c
      do 400 j=1,npix
      do 400 i=1,npix
  400 outmap((j-1)*npix+i)=map3((1024-npix)/2+i,(1024-npix)/2+j)
      outfft=outf(1:loutf)//'/bfz'//datetime
      call ftinit(iounit,outfft,2880,istat) 
      if (istat.ne.0) go to 900
      call puthdf(iounit,nbit,naxis,naxes,
     -      dtobs,tmobs,strobs,endobs,jstdt,jsttm,jsts,jste,
     -      nfrst,nfrend,0,natt,nfalt,nfreq,nfrstat,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,2.45553,2.45553,1)
      call putpdf(iounit,iqb,1.0125,ical,cradd,crfac,mbcyn,msrstr,
     -    dpxew,dpxns,dcorf,dskbr,cra,nmaxa)
      call putcon(iounit,iqb)
      call putcpr(iounit,datecal,timecal,iccal,ictgr,
     -    cpxew,cpxns,ccorf)
      call putsft(iounit,ipix,spxew,spxns,corfmxz)
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
      write(6,'('' out file  :'',a40)') outfft
c
  100 continue
c
  900 write(6,'('' istat ='',i8)') istat
c
      end
c
c
***********************************************************************
*     calcor1024 - correlation calculation                            *
*       map1   : input,real(1024,1024) map1                           *
*       map2   : input,real(1024,1924) map2                           *
*       imx    : input,integer x maximum position                     *
*       lmx    : input,integer y maximum position                     *
*       ipix   : input,integer correlation size                       *
*       corf   : output,real correlation of map1 and map2             *
*                                   February 5, 2007  H. Koshiishi    *
***********************************************************************
c
      subroutine calcor1024(map1,map2,imx,lmx,ipix,corf)
c
      real sum1,sum2,sum,corf
      real map1(1024,1024),map2(1024,1024)
c
      sum1=0.
      sum2=0.
      sum=0.
c
      do 100 i=imx-ipix,imx+ipix
      do 100 j=lmx-ipix,lmx+ipix
c
         sum1=sum1+map1(i,j)*map1(i,j)
         sum2=sum2+map2(i,j)*map2(i,j)
c
         sum=sum+map1(i,j)*map2(i,j)
c
  100 continue
c
      corf=sum/sqrt(sum1)/sqrt(sum2)
c
      return
      end
c
c
***********************************************************************
*      smaxa1 - search maximum position                               *
*                                                                     *
*         map      :  input,real(41,41) map                           *
*         ixmx     :  output,integer x maximum position               *
*         iymx     :  output,integer y maximum position               *
*         amx      :  output,real maximum value                       *
*                                   April 22, 2007  H. Koshiishi      *
***********************************************************************

           subroutine smaxa1(map,ixmx,iymx,amx)

           integer i,j,ixmx,iymx
           real map,amx
           dimension map(41,41)

           amx=0.

           do 100 i=1,41
           do 100 j=1,41

              if (map(i,j).gt.amx) then

                 ixmx=i
                 iymx=j
                 amx=map(i,j)

              endif

100        continue

           return
           end
c
c
***********************************************************************
*      smaxa2 - search maximum position                               *
*                                                                     *
*         a        :  input,real(9)                                   *
*         amx      :  output,real maximum value                       *
*         xmx      :  output,real x maximum position                  *
*         ymx      :  output,real y maximum position                  *
*                                   April 16, 2007  H. Koshiishi      *
***********************************************************************

           subroutine smaxa2(a,amx,xmx,ymx)

           integer i,j
           real a,xmx,ymx,amx,map
           dimension a(9),map(41,41)

           do 100 i=1,41
           do 100 j=1,41

              x=(real(i)-21.)/10.
              y=(real(j)-21.)/10.
              
              map(i,j)=a(1)+
     -                 a(2)*y+
     -                 a(3)*y*y+
     -                 a(4)*x+
     -                 a(5)*x*y+
     -                 a(6)*x*y*y+
     -                 a(7)*x*x+
     -                 a(8)*x*x*y+
     -                 a(9)*x*x*y*y

100        continue

           amx=0.

           do 200 i=1,41
           do 200 j=1,41

              if (map(i,j).gt.amx) then

                 xmx=(real(i)-21.)/10.
                 ymx=(real(j)-21.)/10.
                 amx=map(i,j)

              endif

200        continue

           return
           end
c
c
***********************************************************************
*     puthdf - write standard header information to fits file         *
*       iunit  : input,integer fortran unit number for fits file      *
*       nbit   : input,integer bitpix of fits header                  *
*       naxis  : input,integer naxis of fits header                   *
*       naxes  : input,integer() naxis1,naxis2,... of fits header     *
*       dtobs  : input,character*18 observation date in ut            *
*       tmobs  : input,character*18 center of observation time in ut  *
*       strobs : input,character*18 start of observation time in ut   *
*       endobs : input,character*18 end of observation time in ut     *
*       jstdt  : input,character*18 observation date in jst           *
*       jsttm  : input,character*18 observation date in jst           *
*       jsts   : input,character*18 start of observation time in jst  *
*       jste   : input,character*18 end of observation time in jst    *
*       nfrst  : input,integer first frame number                     *
*       nfrend : input,integer last frame number                      *
*       npol   : input,integer polarization r:1,l:-1,r+l:0,r-l:2,     *
*                              l-r:-2,r|l:10,(r-l)/(r+l):12           *
*       natt   : input,character*18 attenuator                        *
*       nfalt  : input,character*18 frequency alternation             *
*       nfreq  : input,character*18 observation frequency             *
*       nfrstat: input,character*18 frame status                      *
*       dattyp : input,character*18 data type                         *
*                                      March 6, 2007  H. Koshiishi    *
***********************************************************************
c
      subroutine puthdf(iunit,nbit,naxis,naxes,
     -                  dtobs,tmobs,strobs,endobs,
     -                  jstdt,jsttm,jsts,jste,
     -                  nfrst,nfrend,npol,natt,nfalt,
     -                  nfreq,nfrstat,dattyp)
c
      integer naxes(2)
      character*18 dtobs,tmobs,strobs,endobs,jstdt,jsttm,jsts,jste,
     -     pol,natt,nfalt,nfreq,nfrstat,dattyp
      character*44 cmnt
c
*****
c
      cmnt=' '
c
***** write primary header information *****
c
      if (npol.eq.1) then
       write(pol,'(''rcp'')')
      else if (npol.eq.-1) then
       write(pol,'(''lcp'')')
      else if (npol.eq.0) then
       write(pol,'(''r+l'')')
      else if (npol.eq.2) then
       write(pol,'(''r-l'')')
      else if (npol.eq.-2) then
       write(pol,'(''l-r'')')
      else if (npol.eq.10) then
       write(pol,'(''r|l'')')
      else if (npol.eq.12) then
       write(pol,'(''(r-l)/(r+l)'')')
      else
       write(pol,'(''unknown'')')
      end if
c
      call ftpprh(iunit,1,nbit,naxis,naxes,0,0,0,istat)
c
      call ftpkys(iunit,'date-obs',dtobs,cmnt,istat)
      call ftpkys(iunit,'time-obs',tmobs,cmnt,istat)
      call ftpkys(iunit,'strt-obs',strobs,cmnt,istat)
      call ftpkys(iunit,'end-obs',endobs,cmnt,istat)
      call ftpkys(iunit,'jstdate',jstdt,cmnt,istat)
      call ftpkys(iunit,'jsttime',jsttm,cmnt,istat)
      call ftpkys(iunit,'jst-strt',jsts,cmnt,istat)
      call ftpkys(iunit,'jst-end',jste,cmnt,istat)
      call ftpkyj(iunit,'startfrm',nfrst,cmnt,istat)
      call ftpkyj(iunit,'endfrm',nfrend,cmnt,istat)
c
      call ftpkys(iunit,'polariz',pol,cmnt,istat)
      call ftpkys(iunit,'att-10db',natt,cmnt,istat)
      call ftpkys(iunit,'obs-mode',nfalt,cmnt,istat)
      call ftpkys(iunit,'obs-freq',nfreq,cmnt,istat)
      call ftpkys(iunit,'frm-stat',nfrstat,cmnt,istat)
      call ftpkys(iunit,'data-typ',dattyp,cmnt,istat)
c
      call ftpkys(iunit,'object','sun',cmnt,istat)
      call ftpkys(iunit,'telescop','radioheliograph',cmnt,istat)
      call ftpkys(iunit,'origin','nobeyama radio obs',cmnt,istat)
c
      call ftpkys(iunit,'hdrident','HeliogFITS 2.0',cmnt,istat)
c
      if (istat.ne.0) write(6,'('' puthdf : istat ='',i6)') istat
      return
      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','burst34 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
***********************************************************************
*     putsft - write shift parameters to fits file                    *
*       iunit  : input,integer fortran unit number for fits file      *
*       ipix   : input,integer correlation size                       *
*       spxew  : input,real x-offset in shift                         *
*       spxns  : input,real y-offset in shift                         *
*       scorf  : input,real correlation in shift                      *
*                                      March 6, 2007  H. Koshiishi    *
***********************************************************************
c
      subroutine putsft(iunit,ipix,spxew,spxns,scorf)
c
      real spxew,spxns,scorf
      character*44 cmnt
c
      cmnt=' '
c
      call ftpkyj(iunit,'ipix',ipix,
     -     'correlation size in shift',istat)
      call ftpkyf(iunit,'sddoff1',spxew,2,
     -     'x-offset in shift',istat)
      call ftpkyf(iunit,'sddoff2',spxns,2,
     -     'y-offset in shift',istat)
      call ftpkyf(iunit,'sddcorr',scorf,4,
     -     'correlation in shift',istat)
c
      if (istat.ne.0) write(6,'('' putsft : istat ='',i6)') istat
c
      return
      end
c
c
***********************************************************************
*     mean1024 - search maximum point                                 *
*       map    : input,real(1024,1024)                                *
*       imax   : output,integer                                       *
*       jmax   : output,integer                                       *
*                                      March 27, 2007  H. Koshiishi   *
***********************************************************************
c
      subroutine mean1024(map,imax,jmax)
c
      real map(1024,1024)
c
      amx=0.
c
      do 100 i=1,1024
      do 100 j=1,1024
c
         if (map(i,j).gt.amx) then
c
            imax=i
            jmax=j
            amx=map(i,j)
c
         endif
c
  100 continue
c
      return
      end
c
c
***********************************************************************
*       contour_trim4 - serach contour trim                           *
*                                                                     *
*       map    : input,real(1024,1024)  dirty map                     *
*       trim   : input,real trim level                                *
*       imx    : input,integer x-coordinate of the maximum point      *
*       lmx    : input,integer y-coordinate of the maximum point      *
*       ipix   : output,integer number of pixels                      *
*                                      April 19, 2007  H. Koshiishi   *
***********************************************************************

          subroutine contour_trim4(map,trim,imx,lmx,ipix)

          integer key
          integer i,j,k,l,imx,lmx,ib,lb
          integer ipix
          integer n,buffer_num
          integer ibuffer(1048576),jbuffer(1048576)
          integer map_status(1024,1024) 
          real trim
          real map(1024,1024)

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

                map_status(i,j)=0

100          continue
200       continue

          do 300 n=1,1048576

             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+1023,1024)+1
                      lb=mod(jbuffer(n)+l+1023,1024)+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

          ipix=buffer_num

          return
          end
c
c
***********************************************************************
*      shift1024 - position shift                                     *
*                                                                     *
*         ioffm    :  input,integer                                   *
*         loffm    :  input,integer                                   *
*         map      :  input,real(1024,1024)                           *
*                     output,real(1024,1024)                          *
*                                      March 2, 2007  H. Koshiishi    *
***********************************************************************

           subroutine shift1024(ioffm,loffm,map)

           integer i,j,ib,jb,ioffm,loffm
           real map,shift_map
           dimension map(1024,1024),shift_map(1024,1024)


           do 200 i=1,1024

              do 100 j=1,1024

                 ib=mod(i+ioffm+1023,1024)+1
                 jb=mod(j+loffm+1023,1024)+1

                 shift_map(i,j)=map(ib,jb)

100           continue

200        continue

           do 400 i=1,1024

              do 300 j=1,1024

                 map(i,j)=shift_map(i,j)

300           continue

400        continue

           return
           end