program snap2d *********************************************************************** * create cleaned maps (1024x1024 etc.) * * ver 6.2 April 19, 1999 Y. Hanaoka * *********************************************************************** c integer naxes(2),nantst(84),iimg(1000), - pi(2000),pl(2000),pw(2000) 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(2000) 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/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 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 : ''$)') read(5,*) cradd 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(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 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 do 180 il=1,512,64 c 180 write(6,'('' '',8f9.2)') (mapa(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)') solp,solb,alt,az ccc write(6,'('' matrix ='',4f10.4)') (pmat(i),i=1,4) 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 clean34(menu,mapa,pmat,solr*1.0125, - nantst,gbeam,cradd,crfac,cra,pxew,pxns,corf, - nmaxa,bgr,dskbr,pi,pl,pw,pa,cmap) c do 350 il=1,512,64 c 350 write(6,'('' '',8f9.2)') (cmap(i,il),i=1,512,64) c if (mbcyn.eq.1) then pdskbr=-dskbr else pdskbr=dskbr endif call proj34(pmat,solp,pdskbr,pxew,pxns,xoffp,yoffp,cmap) c do 370 il=1,512,64 c 370 write(6,'('' '',8f9.2)') (cmap(i,il),i=1,1024,128) 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. 1024) 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)=cmap((1024-npix)/2+i,(1024-npix)/2+j) outfft=outf(1:loutf)//'/i'//fform//'z'//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,az,pmat) call putcor(iounit,xoffa,yoffa,real(npix)/2.+0.5, - real(npix)/2.+0.5,2.45553,2.45553,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 100 continue c 900 write(6,'('' istat ='',i8)') istat c end c c *********************************************************************** * putpdp - write promgram dependent parameters to fits file * * clean and projection parameters * * 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 clean criterion input * * 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 * * criter : 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,criter,ncmp) c character*44 cmnt c cmnt=' ' c call ftpkys(iunit,'progname','snap2d34 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,'dskbr',dskbr,2, - 'brightness of the dirty disk',istat) call ftpkyf(iunit,'ddcorr',corf,4, - 'correlation between dirty disk and model',istat) call ftpkyf(iunit,'criter',criter,2,'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