program snap2d17suw
**************************************************************
*     create cleaned maps (512x512)                          *
**************************************************************

      integer naxes(2),nantst(84),iimg(1000),pir(2000),pil(2000),
     -     plr(2000),pll(2000),pwr(2000),pwl(2000)
      integer*2 icald(336),icald0(336)
      real gbeam(512,512,0:3),flxr(84),flxl(84),
     -     mapra(512,512),maprb(512,512),mapla(512,512),maplb(512,512),
     -     cmapr(512,512),cmapl(512,512),pmat(4),
     -     par(2000),pal(2000)
      real rlph(84)
      character*80 filehd,infpar,outf,outfft
      character*18 dattyp
      character*5 fnum
      character*17 datetime
c
c
      integer naxes1(2),isam
      real mapra1(512,512),dmy512(512,512),gaine,
     -     mapla1(512,512),mask_all(512,512),
     -     xd,taper,mask_super(512,512),gbeam2(512,512,0:5),
     -     mapra3(1024,1024),cbeam2(256,256,0:5),
     -     dummy(1024,1024),
     -     mapra4(256,256),maprb4(256,256),
     -     mapla4(256,256),maplb4(256,256),
     -     mask_pfi(512,512),
     -     mapra2(128,128),mapla2(128,128),
     -     suma(128,128),sumb(128,128)

      character*80 outfft2
      character*60 comment
      rd=57.29577951
      pi=3.141592653
c
c
ccC   write(6,*)'*****************************************'
ccc   write(6,*)' OUTPUT FILES'
ccc   write(6,*)'*****************************************'
ccc   write(6,'(/'' input output filename                  : ''$)')
      read(5,'(a)') outf
c
ccc   write(6,*)'*****************************************'
ccc   write(6,*)' DATA AND FREME PART'
ccc   write(6,*)'*****************************************'
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,*)'*****************************************'
ccc   write(6,*)' IMAGE POSITION AND SIZE PART'
ccc   write(6,*)'*****************************************'
ccc   write(6,*) 'PFI center from sun center'
c     write(6,'(''0(deg.)/1(pixel(4.911)),x,y): ''$)')
ccc   write(6,'(''1(arcsec)/2(pixel(4.911)),x,y): ''$)')
      read(5,*) ictype,cenx,ceny
      if(cenx.ne.999.and.ceny.ne.999) then
       if(ictype.eq.1) then
         cenx=cenx/4.911+257.
         ceny=ceny/4.911+257.
       else
cvvvvvvvvvvv T. Yokoyama 1999-2-4
c        cenx=cenx+257.
c        ceny=ceny+257.
         cenx=0.5*cenx+257.
         ceny=0.5*ceny+257.
c^^^^^^^^^^^ T. Yokoyama 1999-2-4
       endif
      endif
c      write(6,*)'cenx,ceny',cenx,ceny
c
ccc   write(6,'('' input first frame number               : '')')
ccc   write(6,'('' for irregular interval, you set 0      : ''$)')
      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=1
       do 150 i=nfrst,nfrend,itvl
        nimg=nimg+1
        iimg(nimg)=i
  150  continue
      else
ccc    write(6,'('' input total number of frames            : ''$)')
       read(5,*) nimg
       nimg=nimg+1
ccc    write(6,'('' input frame numbers                    : ''$)')
       read(5,*) (iimg(i),i=2,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 reference frame number           : '')')
      read(5,*) nfrq
ccc   write(6,'('' input integration frames for reference'')')
ccc   write(6,'(''      for calibration                   : '')')
      read(5,*) icalq
ccc   write(6,'(''      for image processing              : '')')
      read(5,*) itgrq
      iimg(1)=nfrq
c
ccc   write(6,*)'*****************************************'
ccc   write(6,*)' CLEAN PARAMETER'
ccc   write(6,*)'*****************************************'
ccc   write(6,'('' CLEAN criterion for pre CLEAN(acp,scp)     : ''$)')
      read(5,*) pcrit_r,pcrit_l
ccc   write(6,'('' CLEAN Criterion for ACP,SCP                : ''$)')
      read(5,*) cradd_r,cradd_l
ccc   write(6,'('' CLEAN Criterion ratio                      : ''$)')
      read(5,*) c_ratio
c
c
***** default parameters *****
      cradd=30000.
      crsub=7500.
      xoff0=999.
      yoff0=999.
      ioffp=0
      loffp=0
      mbcyn=1
c
c      write(6,'('' CLEAN PFI dirty map (Y:1 / N:0)  : ''$)')
c      read(5,*) i_clean

* CLEAN negative components (Y:1 / N:0)
      negative=1

      gaine=0.2
      c_ratio=0.005

* integration frames after restoration
      itgra=1
* image type(RCP/LCP:1, ACP/SCP:2)
      image_type=2
c
      i_alias=1
* Use Steer Algorithm(ON:1 / OFF:0)
      istr_r=1
      istr_l=0
      taper=1.0
      xd=160.
* LOOP GAIN FOR CLEAN
      gain=0.2
* output intrium map (yes=1/no=0)
      intr=0
* CELL SIZE FOR SUPER-UNIFORM WEIGHTING
      icell=5

c
      loutf=index(outf,' ')-1
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
*****      Initialize FFT table
      call cft1024(mapra3,dummy,0)
      call cft512(mapra,maprb,0)
      call full_uv(nantst,mask_all)
      call make_super(mask_super,mask_all,icell)
c      outfft2=outf(1:loutf)//'mask_super.fits'
c      call wt_fits(512,mask_super,outfft2)
      call msrc512(nantst,gbeam)
      call msrc_super(nantst,gbeam2,cbeam2,mask_super)
c
      iounit=2
      iounit2=2

      nbit=-32
      naxis=2
      naxes(1)=512
      naxes(2)=512
      naxes1(1)=128
      naxes1(2)=128
c
ccc   write(6,*)'start imaging'
      do 100 iifr=1,nimg
       iframe=iimg(iifr)
c    
       do j=1,128
       do i=1,128
          suma(i,j)=0.
          sumb(i,j)=0.
       enddo
       enddo

***** loop for integration of restored images
       if(iifr.eq.1) then 
         itgrav=1
       else
         itgrav=itgra
       end if

       do 99 isum=1,itgrav
ccc    write(6,*) isum,' th image / ',itgrav
       if (iifr.eq.1) then
         icalv=icalq
         itgrv=itgrq
         itgrb=1
         if(itgr.lt.0) then
           itgrb=-itgrq
           itgrv=1
         end if
       else
         icalv=ical
         itgrv=itgr
         itgrb=1
         if(itgr.lt.0) then
           itgrb=-itgr
           itgrv=1
         end if
       end if
       iframe1=iframe+itgrv*(isum-1)
       iframe2=iframe1+itgrv-1

cc
       if(iifr.gt.1 .and. ical.le.0) then
          do i=0,336
             icald(i)=icald0(i)
          enddo
       end if
c
***** initialize rlph *****
       if (iifr.eq.1) rlph(1)=1.
c
c$$$       call snapuv17(iunit,filehd,itgrb,
c$$$     -      iframe1,iframe2,
c$$$     -      nantst,icalv,icald,rlph,
c$$$     -      ndut,msut,ndjst,msjst,msjsts,msjste,natt,nfalt,nfreq,
c$$$     -      flxr,flxl,mapra,maprb,mapla,maplb,nstat)
       call snapuv17(iunit,filehd,iframe1,iframe2,
     -     nantst,icalv,icald,rlph,1.,
     -     ndjst,msjst,msjsts,msjste,natt,nfalt,nfreq,nstat,
     -     flxr,flxl,mapra,maprb,mapla,maplb,nfrmst)

      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 (nstat.ne.3) then
       ndchar=13
      else
       ndchar=17
      endif

c$$$       write(6,*) 'end snapuv17',msut,ndut
c$$$       iyyyy=ndut/10000
c$$$       if(iyyyy.gt.69) then
c$$$          iyyyy=iyyyy+1900
c$$$       else
c$$$          iyyyy=iyyyy+2000
c$$$       endif
c$$$       immdd=mod(ndut,10000)
c$$$       write(fdate,'(i4.4,i4.4,a1,i6.6)') iyyyy,immdd,'_',
c$$$     -      int(msut/1000.)
c$$$       write(6,*)'fdate= ',fdate

       if(iifr.gt.1 .and. ical.lt.0) then
         ffact=1.
         if (natt.eq.10) ffact=9.8
         flxr1=flxr(1)*ffact
         flxl1=flxl(1)*ffact
       end if
c
c     array for 8Nd or 16Nd sampling data
c
       call cft512(mapra,maprb,-1)
       call cft512(mapla,maplb,-1)
       do j=1,512
       do i=1,512
          mapra1(i,j)=mapra(i,j)
          mapla1(i,j)=mapla(i,j)
       enddo
       enddo
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)
c       call caleph(dec1,dec2,dec3,ha1,ha2,ha3,
c     -     solr1,solr2,solr3,solp1,solp2,solp3,solb1,solb2,solb3,
c     -     msjst,dec,ha,solr,solp,solb,az,alt,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)') solp,solb,alt,az
ccc    write(6,'('' matrix ='',4f10.4)') (pmat(i),i=1,4)

       if(iifr.eq.1 .and. ical.le.0) then
         do i=1,336
            icald0(i)=icald(i)
         enddo
         alt0=alt
         solr0=solr
         flxr0=flxr(1)
         flxl0=flxl(1)
       end if

      if(iifr.gt.1 .and. ical.le.0) then
        bgrr=bgrr0
        dskr=dskr0
        bgrl=bgrl0
        dskl=dskl0
        dskbr=dskbr0
        pxew=pxew0
        pxns=pxns0
        go to 399
      end if
c
***** preclean process for brightness normalization *****
c
      nmaxa=2000
      nmaxs=2000
c
      if(iifr.eq.1 .or. ical.gt.0) then
        do j=1,512
           do i=1,512
              maprb(i,j)=mapra(i,j)
              maplb(i,j)=mapla(i,j)
           enddo
        enddo

        call clean17_512(3,maprb,pmat,solr*1.0125,
     -       nantst,gbeam,0.,0.,0.,
     -       nmaxa,bgrr,dskr,pir,plr,pwr,par,cmapr)
c

        call clean17_512(3,maplb,pmat,solr*1.0125,
     -       nantst,gbeam,0.,0.,0.,
     -       nmaxs,bgrl,dskl,pil,pll,pwl,pal,cmapl)
        if(iifr.eq.1) then
          bgrr0=bgrr
          dskr0=dskr
          bgrl0=bgrl
          dskl0=dskl
        end if
      end if
      
***** add and subtract rcp/lcp *****
c
      do j=1,512
      do i=1,512
         maprb(i,j)=mapra(i,j)-bgrr-(mapla(i,j)-bgrl)/dskl*dskr
         mapra(i,j)=(mapra(i,j)+mapla(i,j)/dskl*dskr)/2.
         mapla(i,j)=maprb(i,j)/2.
      enddo
      enddo
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 (nstat.eq.0) go to 100
c

      call clean17_512(2,mapra,pmat,solr*1.0125,
     -     nantst,gbeam,cradd,
     -     pxew,pxns,nmaxa,bgr,dskbr,pir,plr,pwr,par,cmapr)

      if(ical.le.0 ) then
        pxew0=pxew
        pxns0=pxns
        dskbr0=dskbr
      end if
c
***** first frame data, FFI images are made

      if (iifr.gt.0.or.isum.gt.1) go to 399
c
      call clean17_512(-1,mapla,pmat,solr*1.0125,
     -     nantst,gbeam,crsub,
     -     0.,0.,nmaxs,bgr,dskbr,pir,plr,pwr,par,cmapl)
      if (mbcyn.eq.1) dskbr=-dskbr

      call proj17(pmat,solp,dskbr,pxew,pxns,0.,0.,cmapr)
      call proj17(pmat,solp,dskbr,pxew,pxns,0.,0.,cmapl)
c
      write(fnum,'(i5.5)') iframe
      if(iifr.eq.1) then
        outfft=outf(1:loutf)//'/ifa'//datetime(1:ndchar)
c        outfft=outf(1:loutf)//'/'//'ipa'//fdate//'_reff'
      else
        outfft=outf(1:loutf)//'/ifa'//datetime(1:ndchar)
c        outfft=outf(1:loutf)//'/'//'ipa'//fdate
      end if
      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 puthdr(iounit,nbit,naxis,naxes,ndut,msut,
c     -     ndjst,msjst,msjsts,msjste,iframe,iframe+itgr-1,
c     -     0,natt,nfalt,nfreq,dattyp)
c
      if (cradd.ge.0.) then
       cr=cradd
      else
       cr=abs(cradd*dskbr)
      end if
      call putclp(iounit,nfreq,cr,nmaxa,pxew,pxns,ioffp,loffp,mbcyn)
      call putoep(iounit,solr,solp,solb,dec,ha,az,alt,za,pmat)
c      call putoep(iounit,solr,solp,solb,dec,ha,az,alt,pmat)
c
      call ftpdef(iounit,nbit,naxis,naxes,0,0,istat)
      call ftppre(iounit,0,1,naxes(1)*naxes(2),cmapr,istat)
      call ftclos(iounit,istat)
      if (istat.ne.0) go to 900
c     write(6,'('' file '',a40,'' created'')') outfft
      write(6,'('' out file '',a60)') outfft
c
      if(iifr.eq.1) then
        outfft=outf(1:loutf)//'/ifs'//datetime(1:ndchar)
c        outfft=outf(1:loutf)//'/'//'ips'//fdate//'_reff'
      else
        outfft=outf(1:loutf)//'/ifs'//datetime(1:ndchar)
c        outfft=outf(1:loutf)//'/'//'ips'//fdate
      end if
      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 puthdr(iounit,nbit,naxis,naxes,ndut,msut,
c     -     ndjst,msjst,msjsts,msjste,iframe,iframe+itgr-1,
c     -     2,natt,nfalt,nfreq,dattyp)
c
      if (crsub.ge.0.) then
       cr=crsub
      else
       cr=abs(crsub*dskbr)
      end if
      call putclp(iounit,nfreq,cr,nmaxs,pxew,pxns,ioffp,loffp,mbcyn)
      call putoep(iounit,solr,solp,solb,dec,ha,az,alt,za,pmat)
c      call putoep(iounit,solr,solp,solb,dec,ha,az,alt,pmat)
c
      call ftpdef(iounit,nbit,naxis,naxes,0,0,istat)
      call ftppre(iounit,0,1,naxes(1)*naxes(2),cmapl,istat)
      call ftclos(iounit,istat)
      if (istat.ne.0) go to 900
c     write(6,'('' file '',a40,'' created'')') outfft
      write(6,'('' out file '',a60)') outfft
c
  399 continue
c-------------------------------------------------------------
c-------------------------------------------------------------
c     16d imaging
c-------------------------------------------------------------
c-------------------------------------------------------------
c      if (mbcyn.eq.1) dskbr=-dskbr
      if (dskbr.lt.0) dskbr=-dskbr

***** Search max point of acp map
      if (iifr.eq.1) then
        if(i_alias.eq.1) then
          call make_mask(-1,mask_pfi,mapra,0.,0.)
        end if

        if(cenx.ne.999. .and. ceny.ne.999.) then
          xmc0=cenx
          ymc0=ceny
        else
          call mean(1,cmapr,amr,amxr,imxc,lmxc)
          xmc0=imxc
          ymc0=lmxc
        end if
        write(comment,*)'center coordinate(65,65) = ',xmc0,ymc0
      end if
***** coordinate of final map
      xmc=xmc0
      ymc=ymc0
c
***** max of CLEAN components
      nmax16r=5000
      nmax16l=5000
***** Caliculate max point of dirty map
ccc   write(6,*) solp,pmat
ccc   write(6,*)'initial xmc,ymc(512x512) ',xmc,ymc
      call cmap2dmap(170,pmat,solp,xmc,ymc)
ccc   write(6,*)'projected xmc,ymc(1024x1024) ',xmc,ymc

c
***** model sun without offset
       xs=xmc/2.-257.
       ys=ymc/2.-257.
c
c      maprb:dirty disk
c      maplb:flat disk
c
       call msun17(pmat,solr*1.0125,
     -     nantst,0.,0.,1.,maprb,maplb,mstype)
c
c      outfft2='mdisk0.fits'
c      call wt_fits(512,maprb,outfft2)
       call cft512(maprb,maplb,1)
       call jitter16(maprb,maplb,xs,ys)
       call cft512(maprb,maplb,-1)
c      outfft2='mdisk.fits'
c      call wt_fits(512,maprb,outfft2)
c
***** jitter correction

c      outfft2='pre_jitter.fits'
c      call wt_fits(512,mapra1,outfft2)
       call cft512(mapra1,dmy512,1)
       call jitter16(mapra1,dmy512,(pxew+xs),(pxns+ys))
       call cft512(mapra1,dmy512,-1)
c      outfft2='post_jitter.fits'
c      call wt_fits(512,mapra1,outfft2)

       call cft512(mapla1,dmy512,1)
       call jitter16(mapla1,dmy512,(pxew+xs),(pxns+ys))
       call cft512(mapla1,dmy512,-1)
c
c
c      maprb:dirty disk
***** R+L R-L
       if(ical.ge.0) then
         do j=1,512
         do i=1,512
            dmy512(i,j)=mapra1(i,j)-bgrr-(mapla1(i,j)-bgrl)/dskl*dskr
            mapra1(i,j)=((mapra1(i,j)+mapla1(i,j)/dskl*dskr)/2.-bgr)
     -           *10000./dskbr
            mapla1(i,j)=dmy512(i,j)/2.0*10000./dskbr
         enddo
         enddo
         if(iifr.ne.-1) then 
            do j=1,512
            do i=1,512
               mapra1(i,j)=mapra1(i,j)-maprb(i,j)/100.*10000.
            enddo
            enddo
         end if
       else

         fnorm=flxl1/flxl0*flxr0/flxr1*dskr0/dskl0
         dskbr=dskr*sin(alt0/rd)/sin(alt/rd)/(flxr1/flxr0)
         do j=1,512
         do i=1,512
            mapra1(i,j)=(mapra1(i,j)-flxr0/flxr1*bgrr)
            mapla1(i,j)=(mapla1(i,j)-flxl0/flxl1*bgrl)

            mapra1(i,j)=mapra1(i,j)/(flxr0/flxr1*dskr)*10000.
            mapla1(i,j)=mapla1(i,j)/(flxl0/flxl1*dskl)*10000.
            dmy512(i,j)=(mapra1(i,j)-mapla1(i,j))/2.
            mapra1(i,j)=mapra1(i,j)-maprb(i,j)/100.*10000.
            mapla1(i,j)=mapla1(i,j)-maprb(i,j)/100.*10000.
            mapra1(i,j)=(mapra1(i,j)+mapla1(i,j))/2.
            mapla1(i,j)=dmy512(i,j)
         enddo
         enddo
       end if
c

c
***** 
c
***** set pre CLEAN mode
      if(image_type.eq.1) then
        if(negative.eq.1) then
          menur=-1
          menul=-1
        else
          menur=1
          menul=1
        end if
      else
        if(negative.eq.1) then
          menur=-1
          menul=-1
        else
          menur=1
          menul=-1
        end if
      end if
c      outfft2='pre_pcleana.fits'
c      call wt_fits(512,mapra1,outfft2)
c      outfft2='pre_pcleans.fits'
c      call wt_fits(512,mapla1,outfft2)

      if(pcrit_r.gt.0) then   
        nmaxp=500000
        call max512(-1,mapra1,amx0,imx,lmx)
        amx0=amx0*c_ratio
        amxpr=amax1(amx0,pcrit_r)
c        amx0=maxval(mapra1)*c_ratio
c        amx1=abs(minval(mapra1))*c_ratio
c        amx=amax1(amx0,pcrit_r)
c        amxpr=amax1(amx,amx1)
          call pre_clean(-1,mapra1,maprb,257.,257.,amxpr,nmaxp,
     -         mask_all,mask_pfi,pmat,solp,iframe,outf,gainp,gbeam,
     -         taper,xd,istr_r)      
      else
         amxpr=pcrit_r
      endif
      if(pcrit_l.gt.0) then   
        nmaxp=500000
        call max512(-1,mapra1,amx0,imx,lmx)
        amx0=amx0*c_ratio
        amxpl=amax1(amx0,pcrit_r)
c        amx0=maxval(mapla1)*c_ratio
c        amx1=abs(minval(mapla1))*c_ratio
c        amx=amax1(amx0,pcrit_l)
c        amxpl=amax1(amx,amx1)
          call pre_clean(-1,mapla1,maplb,257.,257.,amxpl,nmaxp,
     -         mask_all,mask_pfi,pmat,solp,iframe,outf,gainp,gbeam,
     -         taper,xd,istr_l)     
      else
         amxpl=pcrit_l
      endif
c      outfft2='post_pcleana.fits'
c      call wt_fits(512,mapra1,outfft2)
c      outfft2='post_pcleans.fits'
c      call wt_fits(512,mapla1,outfft2)



      call cft512(mapra1,maprb,1)
      call cft512(mapla1,maplb,1)
      ixmc=int(xmc+0.5)
      iymc=int(ymc+0.5)
      ixmc=513
      iymc=513
      do j=1,1024
      do i=1,1024
         mapra3(i,j)=0.
         dummy(i,j)=0.
      enddo
      enddo
      do 151 j=1,512
      do 151 i=1,512
        mapra3(i+256,j+256)=mapra1(i,j)*mask_super(i,j)
        dummy(i+256,j+256)=maprb(i,j)*mask_super(i,j)
  151 continue
      call cft1024(mapra3,dummy,-1)
c     
      do j=-128,127
      do i=-128,127
         mapra4(i+129,j+129)=mapra3(ixmc+i,iymc+j)
      enddo
      enddo
c
      do j=1,1024
      do i=1,1024
         mapra3(i,j)=0.0
         dummy(i,j)=0.0
      enddo
      enddo

      do 152 j=1,512
      do 152 i=1,512
        mapra3(i+256,j+256)=mapla1(i,j)*mask_super(i,j)
        dummy(i+256,j+256)=maplb(i,j)*mask_super(i,j)
  152 continue
      call cft1024(mapra3,dummy,-1)
      do j=-128,127
      do i=-128,127
         mapla4(i+129,j+129)=mapra3(ixmc+i,iymc+j)
      enddo
      enddo

      write(fnum,'(i5.5)') iframe
c      outfft2='pre_256cleana.fits'
c      call wt_fits(256,mapra4,outfft2)
c      outfft2='pre_256cleans.fits'
c      call wt_fits(256,mapla4,outfft2)

c      write(6,*)'Start 256 CLEAN'
      if(cradd_r.gt.0) then
        call max256(-1,mapra4,amx,jmx,jmy)
        amxr1=amax1(amx*c_ratio,cradd_r)
ccc     write(6,*) '256x256 for ACP criterion = ',amxr1
        call fclean256(-1,mapra4,maprb4,gbeam2,cbeam2,5000,amxr1,
     -       icell,gain)
      else
        do j=1,256
        do i=1,256
           maprb4(i,j)=mapra4(i,j)
        enddo
        enddo
      end if
      if(cradd_l.gt.0) then
        call max256(-1,mapla4,amx,jmx,jmy)
        amxl1=amax1(amx*c_ratio,cradd_l)
ccc     write(6,*) '256x256 for SCP criterion = ',amxl1
        call fclean256(-1,mapla4,maplb4,gbeam2,cbeam2,5000,amxl1,
     -       icell,gain)
      else
        do j=1,256
        do i=1,256
           maplb4(i,j)=mapla4(i,j)
        enddo
        enddo
      end if
c      write(6,*)'End of 256 CLEAN'

c      outfft2='post_256cleana.fits'
c      call wt_fits(256,mapra4,outfft2)
c      outfft2='post_256cleans.fits'
c      call wt_fits(256,mapla4,outfft2)

      call proj_super(17,pmat,solp,maprb4)
      call proj_super(17,pmat,solp,maplb4)

      do j=-64,63
      do i=-64,63
         mapra2(i+65,j+65)=maprb4(129+i,129+j)
         mapla2(i+65,j+65)=maplb4(129+i,129+j)
      enddo
      enddo
c$$$      do j=-64,63
c$$$      do i=-64,63
c$$$         suma(i+65,j+65)=suma(i+65,j+65)+maprb4(129+i,129+j)
c$$$         sumb(i+65,j+65)=sumb(i+65,j+65)+maplb4(129+i,129+j)
c$$$      enddo
c$$$      enddo
c      mapra2=maprb4(129-64:129+63,129-64:129+63)
c      mapla2=maplb4(129-64:129+63,129-64:129+63)

  99  continue
c$$$      do j=1,128
c$$$      do i=1,128
c$$$         mapra2(i,j)=suma(i,j)/real(itgrav)
c$$$         mapla2(i,j)=sumb(i,j)/real(itgrav)
c$$$      enddo
c$$$      enddo


c     Set filename

      naxes1(1)=128
      naxes1(2)=128
      iounit2=3

      if (iifr.eq.1) then
c       outfft=outf(1:loutf)//'/ipa'//datetime(1:ndchar)//'_ref'
        goto 100
      else 
        outfft=outf(1:loutf)//'/ipa'//datetime(1:ndchar)
      end if
ccc   write(6,*) 'filename= ',outfft
c     Set image type
      if(image_type.eq.1) then
        itype=1
      else
        itype=3
      end if
      isam=8
      if(cradd_r.gt.0) then
        write(dattyp,*) 'cleaned_map'
      else
        write(dattyp,*) 'dirty_map'
      end if
ccc   write(6,*) 'ftinit'
      call ftinit(iounit2,outfft,2880,istat)
      if (istat.ne.0) go to 900
ccc   write(6,*) 'puthdr_pfi'
      call puthdr_pfi(iounit2,nbit,naxis,naxes1,ndut,msut,
     -     ndjst,msjst,msjsts,msjste,iframe,iframe+itgr*itgra-1,
     -     0,natt,nfalt,nfreq,dattyp,itype)
c      call     puthdr(iounit,nbit,naxis,naxes,
c     -     ndjst,msjst,msjsts,msjste,iframe,iframe+itgr-1,
c     -     0,natt,nfalt,nfreq,nstat,dattyp)
c
ccc   write(6,*) 'putclp_pfi'
      call putclp_pfi(iounit2,nfreq,cradd_r,nmaxa,pxew,pxns,ioffp,loffp,
     -     isam,xmc0,ymc0)
c      call putclp_cal(iounit2,icald)

ccc   write(6,*) 'putoep'
      call putoep(iounit2,solr,solp,solb,dec,ha,az,alt,za,pmat)
c      call putoep(iounit2,solr,solp,solb,dec,ha,az,alt,pmat)
c      call putclp_cal(iounit2,icald)
      nfrcal=iimg(iifr)
      call put_imparam(iounit2,nfrcal,icalv,itgrv,
     -     amxpr,amxr1,1.0125,1,
     -     icell,pxew,pxns,abs(dskbr))
c
      call ftpdef(iounit2,nbit,naxis,naxes1,0,0,istat)
      call ftppre(iounit2,0,1,naxes1(1)*naxes1(2),mapra2,istat)
c      call ftpcom(iounit2,pfipos,istat)
      call ftclos(iounit2,istat)

      if (istat.ne.0) go to 900
c     write(6,'('' file '',a40,'' created'')') outfft
      write(6,'('' out file '',a60)') outfft
c
      if (iifr.eq.1) then
        outfft=outf(1:loutf)//'/ips'//datetime(1:ndchar)//'_ref'
      else 
        outfft=outf(1:loutf)//'/ips'//datetime(1:ndchar)
      end if

      if(image_type.eq.1) then
        itype=2
      else
        itype=4
      end if

      if(cradd_l.gt.0) then
        write(dattyp,*) 'cleaned_mao'
      else
        write(dattyp,*) 'dirty_map'
      end if
      call ftinit(iounit2,outfft,2880,istat) 
      if (istat.ne.0) go to 900
      call puthdr_pfi(iounit2,nbit,naxis,naxes1,ndut,msut,
     -     ndjst,msjst,msjsts,msjste,iframe,iframe+itgr*itgra-1,
     -     0,natt,nfalt,nfreq,dattyp,itype)
      call putclp_pfi(iounit2,nfreq,cradd_l,nmaxa,pxew,pxns,ioffp,loffp,
     -     isam,xmc0,ymc0)
c      call putclp_cal(iounit2,icald)
      call putoep(iounit2,solr,solp,solb,dec,ha,az,alt,za,pmat)
c      call putoep(iounit2,solr,solp,solb,dec,ha,az,alt,pmat)
c      call putclp_cal(iounit2,icald)
      nfrcal=iimg(iifr)

      call put_imparam(iounit2,nfrcal,icalv,itgrv,
     -     amxpl,amxl1,1.0125,1,
     -     icell,pxew,pxns,abs(dskbr))
c
      call ftpdef(iounit2,nbit,naxis,naxes1,0,0,istat)
      call ftppre(iounit2,0,1,naxes1(1)*naxes1(2),mapla2,istat)
c      call ftpcom(iounit2,pfipos,istat)
      call ftclos(iounit2,istat)
      if (istat.ne.0) go to 900
c     write(6,'('' file '',a40,'' created'')') outfft
      write(6,'('' out file '',a60)') outfft
c
c
c-------------------------------------------------------------
c-------------------------------------------------------------
c     end of 16d imaging
c-------------------------------------------------------------
c-------------------------------------------------------------
c
c
  100 continue
c
  900 write(6,'('' istat ='',i8)') istat
ccc   write(6,*)'End.'
c
      end
c
*********************************************************************
*     wt_fits: read fits file(n x n array data only)
*      cnum: input - number of column and raw
*      data: input - array data, real data(cnum,cnum)
*      infil: input - FITS filename
*********************************************************************
      subroutine wt_fits(cnum,data,outfft)
c
      implicit NONE
      character*80 outfft,infil
      integer naxes(2),naxis,iounit,istat,nbit,cnum
      real data(cnum,cnum)
c
      iounit=2
      naxes(1)=cnum
      naxes(2)=cnum
      naxis=2
      nbit=-32
c
      infil=outfft
c
      call ftinit(iounit,infil,2880,istat)
      call ftpprh(iounit,1,nbit,naxis,naxes,0,0,0,istat)
      call ftpdef(iounit,nbit,naxis,naxes,0,0,istat)
      call ftppre(iounit,0,1,naxes(1)*naxes(2),data,istat)
      call ftclos(iounit,istat)
c
      write(6,'(''Write FITS file. istat= '',i3,'' file= '',a40)')
     -     istat,infil
c
      outfft=' '
      istat=0
      return 
      end