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