*********************************************************************** * helioglib50 - heliog library v5.0 * * 17/34 GHz Dual Frequency Observation * *********************************************************************** *********************************************************************** * msrc512 - dirty beam of gaussian sources * * nantst : input,integer*(84) antenna status * * gbeam : output,real(512,512,0:3) dirty beam of * * gaussian sources * * width : 0(point source),1,...,10 * * July 13,1992 Y. Hanaoka * * Feb. 27,1999 K. Fujiki * *********************************************************************** c subroutine msrc512(nantst,gbeam) c integer nantst(1) real gbeam(512,512,0:3),a(512,512),b(512,512) c do 100 iws=0,3 c ***** clean source ***** c do 200 il=1,512 do 200 i=1,512 a(i,il)=0. a(i,il)=0. 200 continue c if (iws.eq.0) then a(257,257)=100. else wsr=real(iws) sum=0. do 410 l=1,512 do 410 i=1,512 a(i,l)=a(i,l)+100.*exp(-((257.-real(i))*(257.-real(i)) - +(257.-real(l))*(257.-real(l)))/wsr/wsr) sum=sum+a(i,l) 410 continue do 420 l=1,512 do 420 i=1,512 a(i,l)=a(i,l)*100./sum 420 continue end if c c write(6,'('' '')') c do 400 il=254,260 c 400 write(6,'('' '',7f9.2)') (a(i,il),i=254,260) c ***** dirty map ***** c call cft512(a,b,1) call uvobs(nantst,1.,a,b) call cft512(a,b,-1) c do 300 j=1,512 do 300 i=1,512 300 gbeam(i,j,iws)=a(i,j) c 100 continue c return end *********************************************************************** * putclp - write clean and projection parameters to fits file * * iunit : input,integer fortran unit number for fits file * * nfreq : input,integer observation frequency * * criter : input,real clean criterion * * ncmp : input,integer number of clean components * * pxew : input,real x-offset of the center of dirty disk * * pxns : input,real y-offset of the center of dirty disk * * ioffp : input,integer x-offset of the center of clean disk * * loffp : input,integer y-offset of the center of clean disk * * mbcyn : input,integer main beam correction 1=yes, 2=no * * August 18, 1994 Y. Hanaoka * *********************************************************************** c subroutine putclp(iunit,nfreq,criter,ncmp,pxew,pxns, - ioffp,loffp,mbcyn) c character*44 cmnt c cmnt=' ' c call ftpkys(iunit,'bunit','k','disk = 10000 K',istat) call ftpkyf(iunit,'criter',criter,2,'clean criterion',istat) call ftpkyj(iunit,'ncompo',ncmp, - 'number of clean components',istat) 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) 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 c call ftpkyf(iunit,'crval1',0.,2,'disk center',istat) call ftpkyf(iunit,'crval2',0.,2,'disk center',istat) if (nfreq.eq.17) then call ftpkyf(iunit,'crpix1',257.-ioffp,2,cmnt,istat) call ftpkyf(iunit,'crpix2',257.-loffp,2,cmnt,istat) call ftpkyf(iunit,'cdelt1',4.91104,5,'arcsec',istat) call ftpkyf(iunit,'cdelt2',4.91104,5,'arcsec',istat) else call ftpkyf(iunit,'crpix1',513.-ioffp,2,cmnt,istat) call ftpkyf(iunit,'crpix2',513.-loffp,2,cmnt,istat) call ftpkyf(iunit,'cdelt1',2.45552,5,'arcsec',istat) call ftpkyf(iunit,'cdelt2',2.45552,5,'arcsec',istat) endif call ftpkys(iunit,'ctype1','solar-west',cmnt,istat) call ftpkys(iunit,'ctype2','solar-north',cmnt,istat) c if (istat.ne.0) write(6,'('' putclp : istat ='',i6)') istat return end c c *********************************************************************** * clean17_512 - clean process for 17 GHz images * * model source : gaussian source and flat disk * * menu : input,integer * * 2 : disk/pointsource(+) clean with restoration * * -2 : disk/pointsource(+/-) clean with restoration * * 22 : disk/pointsource(+) clean with point source * * restoration, without disk restoration * * -22 : disk/pointsource(+/-) clean with point source * * restoration, without disk restoration * * -1 : pointsource(+/-) clean with restore * * 1 : pointsource(+) clean with restore * * 3 : preclean only, criter=dskbr*5.+bgr * * map : input,real(512,512) dirty map * * output,real(512,512) resudual * * pmat : input,real(4) projection matrix * * solr : input,real radius of the model sun (arcsec) * * nantst : input,integer*(84) antenna status * * gbeam : input,real(512,512,0:10) dirty beam of gaussian * * sources * * criter : input,real clean criterion * * if < 0, cr=-criter*dskbr * * pxew : input/output,real x-offset of disk center (pixel) * * pxns : input/output,real y-offset of disk center (pixel) * * only valid for disk clean menu * * else when pxew=999. or pxns=999., x-offset and * * y-offset is calculated and output pxew and pxns * * is set to calculated offset. * * nmax : input,integer maximum number of clean components * * : output,integer number of clean components * * bgr : output,real background brightness * * dskbr : input,real disk brightness (valid for criter < 0 * * and menu=+/-1,+/-11 * * : output,real disk brightness * * pi : output,integer(nmax) x-coordinate of clean points * * pl : output,integer(nmax) y-coordinate of clean points * * pw : output,integer(nmax) width of clean points * * pa : output,real(nmax) intensity of clean points * * cmap : output,real(512,512) clean map * * November 25,1995 Y. Hanaoka * * Feb. 27 1999 K. Fujiki * *********************************************************************** c subroutine clean17_512(menu,map,pmat,solr, - nantst,gbeam,criter, - pxew,pxns,nmax,bgr,dskbr,pi,pl,pw,pa,cmap) c integer nantst(84),pi(nmax),pl(nmax),pw(nmax) real map(512,512),cmap(512,512),ddisk(512,512),cdisk(512,512), - cbeam(-30:30,-30:30),gbeam(512,512,0:3), - pmat(4),pa(nmax) c n=0 c if (mod(iabs(menu),10) .ge. 2) then c ********************************************* * preclean process for disk clean * ********************************************* c ***** model sun ***** call msun17(pmat,solr, - nantst,0.,0.,1.,ddisk,cdisk,mstype) c call clhis(map,mstype,300,bgr,dbr) dskbr=dbr-bgr ccc write(6,'('' background/disk(initial) :'',2f10.2)') ccc - bgr,dskbr if (criter.lt.0.) then cr=amax1(dskbr*5.+bgr,dskbr*(-criter)) c write(6,'('' preclean criterion :'',f10.1)') cr else cr=amax1(dskbr*5.+bgr,criter) endif c do 100 while(.true.) c ***** search maximum ***** call mean(menu,map,am,amx,imx,lmx) if (amx.gt.cr .and. n.lt.nmax) then c ***** subtract point source component ***** n=n+1 pi(n)=imx pl(n)=lmx call cleanp_512(2,map,gbeam,amx,imx,lmx,pw(n),pa(n)) c write(6,'('' '',i6,'': am,amx,w,a,imx,lmx ='', c - f10.2,f14.2,i3,f10.2,2i5)') c - n,am,amx,pw(n),pa(n),pi(n),pl(n) c else go to 110 end if c 100 continue 110 if (n.ne.0) then c write(6,'('' n ='',i6,'': am,amx,w,a,imx,lmx ='', c - f10.2,f14.2,i3,f10.2,2i5)') c - n,am,amx,pw(n),pa(n),pi(n),pl(n) else c write(6,'('' n ='',i6)') n end if c ********************************************* * disk clean process * ********************************************* c call clhis(map,mstype,300,bgr,dbr) dskbr=dbr-bgr ccc write(6,'('' background/disk :'',2f10.2)') ccc - bgr,dskbr c if (mod(iabs(menu),10).eq.3) return c ***** subtract disk component ***** c do 140 l=1,512 do 140 i=1,512 140 map(i,l)=map(i,l)-bgr c ***** calculate disk center position ***** if (pxew.eq.999. .or. pxns.eq.999.) then call possh17(map,dskbr,ddisk,pxew,pxns,corf) c write(6,'('' pxew,pxns ='',2f8.3)') pxew,pxns end if c ***** model sun with offset ***** call msun17(pmat,solr, - nantst,pxew,pxns,1.,ddisk,cdisk,mstype) c ***** subtract disk component ***** do 150 l=1,512 do 150 i=1,512 map(i,l)=map(i,l)-dskbr*ddisk(i,l)/100. 150 continue c end if c c ********************************************* * clean process * ********************************************* c ***** clean criterion ***** c if (criter.lt.0.) then cr=dskbr*(-criter) c write(6,'('' clean criterion :'',f10.1)') cr else cr=criter endif c do 160 while(.true.) c ***** search maximum ***** call mean(menu,map,am,amx,imx,lmx) if (abs(amx).gt.cr .and. n.lt.nmax) then c ***** subtract point source component ***** n=n+1 pi(n)=imx pl(n)=lmx call cleanp_512(2,map,gbeam,amx,imx,lmx,pw(n),pa(n)) c write(6,'('' '',i6,'': am,amx,w,a,imx,lmx ='', c - f10.2,f14.2,i3,f10.2,2i5)') c - n,am,amx,pw(n),pa(n),pi(n),pl(n) c else go to 170 end if c 160 continue 170 if (n.ne.0) then c write(6,'('' n ='',i6,'': am,amx,w,a,imx,lmx ='', c - f10.2,f14.2,i3,f10.2,2i5)') c - n,am,amx,pw(n),pa(n),pi(n),pl(n) else c write(6,'('' n ='',i6)') n end if c nmax=n if (abs(menu)/10.eq.1) return c c ********************************************* * restore process * ********************************************* c do 240 l=1,512 do 240 i=1,512 240 ddisk(i,l)=0. c do 500 iw=0,3 c ***** create clean beam ***** c if (iw.eq.0) then do 210 l=-30,30 do 210 i=-30,30 210 cbeam(i,l)=0. cbeam(0,0)=100. else sum=0. do 220 l=-30,30 do 220 i=-30,30 cbeam(i,l)=100.*exp(-(real(i)*real(i)+real(l)*real(l)) - /real(iw)/real(iw)) sum=sum+cbeam(i,l) 220 continue do 230 l=-30,30 do 230 i=-30,30 230 cbeam(i,l)=100./sum*cbeam(i,l) end if c ***** restore point sources ***** c do 250 ip=1,n if (pw(ip).eq.iw) then do 260 l=-30,30 do 260 i=-30,30 c 260 ddisk(pi(ip)+i,pl(ip)+l)=ddisk(pi(ip)+i,pl(ip)+l)+pa(ip) c - *cbeam(i,l)/gbeam(257,257,iw) ii=mod(pi(ip)+i+511,512)+1 ll=mod(pl(ip)+l+511,512)+1 ddisk(ii,ll)=ddisk(ii,ll) - +pa(ip)*cbeam(i,l)/gbeam(257,257,iw) 260 continue c write(6,'('' '',i4,5f10.3)') iw,(cbeam(ii,0),ii=-10,10,5) end if 250 continue c 500 continue c ***** restore disk component ***** c cfact=0. if (iabs(menu).eq.2) cfact=dskbr do 300 l=1,512 do 300 i=1,512 cdisk(i,l)=ddisk(i,l)+cfact*cdisk(i,l)/100. 300 continue c ***** convolution by gaussian beam ***** c do 340 ll=-10,10 do 340 ii=-10,10 340 cbeam(ii,ll)=0.07461*exp(-0.234403 - *(real(ll)*real(ll)+real(ii)*real(ii))) c do 310 l=1,512 do 310 i=1,512 310 cmap(i,l)=0. do 320 l=1,512 do 320 i=1,512 do 330 ll=-10,10 do 330 ii=-10,10 ib=mod(i+ii+511,512)+1 lb=mod(l+ll+511,512)+1 cmap(i,l)=cmap(i,l)+cdisk(ib,lb)*cbeam(ii,ll) 330 continue 320 continue c ***** add residual ***** do 350 l=1,512 do 350 i=1,512 cmap(i,l)=cmap(i,l)+map(i,l) 350 continue c return end c *********************************************************************** * cleanp - subtract a dirty beam component * * map : input,real(512,512) dirty map * * : output,real(512,512) dirty map, point source * * subtracted * * gbeam : input,real(512,512,0:10) dirty beam of gaussian * * sources * * amx : output,real maximum of map * * imx : output,integer x-coordinate of the maximum point * * lmx : output,integer y-coordinate of the maximum point * * iwid : output,integer width of subtracted point source * * ap : output,real intensity of subtracted point source * * July 17,1992 Y. Hanaoka * * Feb 27, 1999 K. Fujiki * *********************************************************************** c subroutine cleanp_512(id,map,gbeam,amx,imx,lmx,iwid,ap) c real map(512,512),gbeam(512,512,0:id) c do 200 iw=1,id do 210 j=-5,5 do 210 i=-5,5 210 if (abs(map(imx+i,lmx+j)) .lt. - gbeam(257+i,257+j,iw)/gbeam(257,257,iw) - *abs(map(imx,lmx))) goto 290 200 continue 290 iwid=iw-1 c ***** loop gain = 0.02 ***** ap=amx*0.02 c do 300 l=1,512 do 300 i=1,512 ib=mod(i-imx+768,512)+1 lb=mod(l-lmx+768,512)+1 300 map(i,l)=map(i,l)-ap*gbeam(ib,lb,iwid)/gbeam(257,257,iwid) c do 400 il=1,512,64 c 400 write(6,'('' '',8f9.2)') (map(i,il),i=1,512,64) c return end c *********************************************************************** * clean34_512 - clean process for 34 GHz images * * model source : gaussian source and flat disk * * menu : input,integer * * 2 : disk/pointsource(+) clean with restoration * * -2 : disk/pointsource(+/-) clean with restoration * * 22 : disk/pointsource(+) clean with point source * * restoration, without disk restoration * * -22 : disk/pointsource(+/-) clean with point source * * restoration, without disk restoration * * -1 : pointsource(+/-) clean with restore * * 1 : pointsource(+) clean with restore * * 3 : preclean only, criter=dskbr*5.+bgr * * map : input,real(512,512) dirty map * * output,real(512,512) resudual * * pmat : input,real(4) projection matrix * * solr : input,real radius of the model sun (arcsec) * * nantst : input,integer*(84) antenna status * * gbeam : input,real(512,512,0:10) dirty beam of gaussian * * sources * * criter : input,real clean criterion * * if < 0, cr=-criter*dskbr * * pxew : input/output,real x-offset of disk center (pixel) * * pxns : input/output,real y-offset of disk center (pixel) * * only valid for disk clean menu * * else when pxew=999. or pxns=999., x-offset and * * y-offset is calculated and output pxew and pxns * * is set to calculated offset. * * nmax : input,integer maximum number of clean components * * : output,integer number of clean components * * bgr : output,real background brightness * * dskbr : input,real disk brightness (valid for criter < 0 * * and menu=+/-1,+/-11 * * : output,real disk brightness * * pi : output,integer(nmax) x-coordinate of clean points * * pl : output,integer(nmax) y-coordinate of clean points * * pw : output,integer(nmax) width of clean points * * pa : output,real(nmax) intensity of clean points * * cmap : output,real(1024,1024) clean map * * November 25,1995 Y. Hanaoka * * Feb 27, 1999 K. Fujiki * *********************************************************************** c subroutine clean34_512(menu,map,pmat,solr, - nantst,gbeam,criter, - pxew,pxns,nmax,bgr,dskbr,pi,pl,pw,pa,cmap) c integer nantst(84),pi(nmax),pl(nmax),pw(nmax) real map(512,512),cmap(1024,1024), - ddisk(512,512),cdisk(1024,1024), - cbeam(-30:30,-30:30),gbeam(512,512,0:3), - pmat(4),pa(nmax) c n=0 c if (mod(iabs(menu),10) .ge. 2) then c ********************************************* * preclean process for disk clean * ********************************************* c ***** model sun ***** call msun34(pmat,solr, - nantst,0.,0.,1.,ddisk,cdisk,mstype) c call clhis(map,mstype,40,bgr,dbr) dskbr=dbr-bgr ccc write(6,'('' background/disk(initial) :'',2f10.2)') ccc - bgr,dskbr if (criter.lt.0.) then cr=amax1(dskbr*5.+bgr,dskbr*(-criter)) c write(6,'('' preclean criterion :'',f10.1)') cr else cr=amax1(dskbr*5.+bgr,criter) endif c do 100 while(.true.) c ***** search maximum ***** call mean(menu,map,am,amx,imx,lmx) if (amx.gt.cr .and. n.lt.nmax) then c ***** subtract point source component ***** n=n+1 pi(n)=imx pl(n)=lmx call cleanp_512(2,map,gbeam,amx,imx,lmx,pw(n),pa(n)) c write(6,'('' '',i6,'': am,amx,w,a,imx,lmx ='', c - f10.2,f14.2,i3,f10.2,2i5)') c - n,am,amx,pw(n),pa(n),pi(n),pl(n) c else go to 110 end if c 100 continue 110 if (n.ne.0) then c write(6,'('' n ='',i6,'': am,amx,w,a,imx,lmx ='', c - f10.2,f14.2,i3,f10.2,2i5)') c - n,am,amx,pw(n),pa(n),pi(n),pl(n) else c write(6,'('' n ='',i6)') n end if c ********************************************* * disk clean process * ********************************************* c call clhis(map,mstype,40,bgr,dbr) dskbr=dbr-bgr ccc write(6,'('' background/disk :'',2f10.2)') ccc - bgr,dskbr c if (mod(iabs(menu),10).eq.3) return c ***** subtract disk component ***** c do 140 l=1,512 do 140 i=1,512 140 map(i,l)=map(i,l)-bgr c ***** calculate disk center position ***** if (pxew.eq.999. .or. pxns.eq.999.) then call possh34(map,dskbr,ddisk,pxew,pxns,corf) c write(6,'('' pxew,pxns ='',2f8.3)') pxew,pxns end if c ***** model sun with offset ***** write(6,*)'OFFSET BY FUJIKI' write(6,*)'OFFSET BY FUJIKI' write(6,*)'OFFSET BY FUJIKI' write(6,*)pxew,pxns call msun34(pmat,solr, - nantst,pxew,pxns,1.,ddisk,cdisk,mstype) c ***** subtract disk component ***** do 150 l=1,512 do 150 i=1,512 map(i,l)=map(i,l)-dskbr*ddisk(i,l)/100. 150 continue c end if c c ********************************************* * clean process * ********************************************* c if (criter.lt.0.) then cr=dskbr*(-criter) c write(6,'('' clean criterion :'',f10.1)') cr else cr=criter endif c do 160 while(.true.) c ***** search maximum ***** call mean(menu,map,am,amx,imx,lmx) if (abs(amx).gt.cr .and. n.lt.nmax) then c ***** subtract point source component ***** n=n+1 pi(n)=imx pl(n)=lmx call cleanp_512(2,map,gbeam,amx,imx,lmx,pw(n),pa(n)) c write(6,'('' '',i6,'': am,amx,w,a,imx,lmx ='', c - f10.2,f14.2,i3,f10.2,2i5)') c - n,am,amx,pw(n),pa(n),pi(n),pl(n) c else go to 170 end if c 160 continue 170 if (n.ne.0) then c write(6,'('' n ='',i6,'': am,amx,w,a,imx,lmx ='', c - f10.2,f14.2,i3,f10.2,2i5)') c - n,am,amx,pw(n),pa(n),pi(n),pl(n) else c write(6,'('' n ='',i6)') n end if c nmax=n if (abs(menu)/10.eq.1) return c c ********************************************* * restore process * ********************************************* c do 240 l=1,512 do 240 i=1,512 240 ddisk(i,l)=0. c do 500 iw=0,3 c ***** create clean beam ***** c if (iw.eq.0) then do 210 l=-30,30 do 210 i=-30,30 210 cbeam(i,l)=0. cbeam(0,0)=100. else sum=0. do 220 l=-30,30 do 220 i=-30,30 cbeam(i,l)=100.*exp(-(real(i)*real(i)+real(l)*real(l)) - /real(iw)/real(iw)) sum=sum+cbeam(i,l) 220 continue do 230 l=-30,30 do 230 i=-30,30 230 cbeam(i,l)=100./sum*cbeam(i,l) end if c ***** restore point sources ***** c do 250 ip=1,n if (pw(ip).eq.iw) then do 260 l=-30,30 do 260 i=-30,30 c 260 ddisk(pi(ip)+i,pl(ip)+l)=ddisk(pi(ip)+i,pl(ip)+l)+pa(ip) c - *cbeam(i,l)/gbeam(257,257,iw) ii=mod(pi(ip)+i+511,512)+1 ll=mod(pl(ip)+l+511,512)+1 ddisk(ii,ll)=ddisk(ii,ll) - +pa(ip)*cbeam(i,l)/gbeam(257,257,iw) 260 continue c write(6,'('' '',i4,5f10.3)') iw,(cbeam(ii,0),ii=-10,10,5) end if 250 continue c 500 continue c ***** restore disk component ***** c cfact=0. if (iabs(menu).eq.2) cfact=dskbr do 300 l=1,1024 do 300 i=1,1024 cdisk(i,l)=ddisk(mod(i+511,512)+1,mod(l+511,512)+1) - +cfact*cdisk(i,l)/100. 300 continue c ***** convolution by gaussian beam ***** c do 340 ll=-10,10 do 340 ii=-10,10 340 cbeam(ii,ll)=0.07461*exp(-0.234403 - *(real(ll)*real(ll)+real(ii)*real(ii))) c do 310 l=1,1024 do 310 i=1,1024 310 cmap(i,l)=0. do 320 l=1,1024 do 320 i=1,1024 do 330 ll=-10,10 do 330 ii=-10,10 ib=mod(i+ii+1023,1024)+1 lb=mod(l+ll+1023,1024)+1 cmap(i,l)=cmap(i,l)+cdisk(ib,lb)*cbeam(ii,ll) 330 continue 320 continue c ***** add residual ***** do 350 l=1,1024 do 350 i=1,1024 cmap(i,l)=cmap(i,l)+map(mod(i+511,512)+1,mod(l+511,512)+1) 350 continue c return end *********************************************************************** * proj16_1 - projection correction for a snapshot map * * solar north pole is to the top of the output image * * intensity is normalized to disk brightness = 10000. * * intensity information from nearby 3x3 pixels * * pixel size of the output map = 4.91104 x 4.91104 arcsec^2 * * (identical to the SXT half resolution images) * * pmat : input,real(4) projection matrix * * solp : input,real polar angle of the sun (degree) * * dskbr : input,real abs(dskbr):disk brightness of * * the snapshot map * * if dskbr<0., main beam correction are done * * ioff : input,integer x-offset of the disk center * * on the original map * * loff : input,integer l-offset of the disk center * * on the original map * * ioffp : input,integer x-offset of the center * * of the projected map from the disk center * * loffp : input,integer l-offset of the center * * of the projected map from the disk center * * map : input,real(512,512) snapshot map * * output,real(512,512) projected map * * April 29,1993 Y. Hanaoka * *********************************************************************** c subroutine proj16_1(ifreq,pmat,solp,dskbr,x0,y0,xx0,yy0,map) c real map(384,384),mapp(384,384),pmat(4) real x0,y0,xx0,yy0 integer ifreq c rd=57.29577951 c cc if(ifreq.eq.17) then do 500 l=1,384 do 500 i=1,384 x=(cos(solp/rd)*real(i-193) - -sin(solp/rd)*real(l-193))*4.91104/4.64947 y=(sin(solp/rd)*real(i-193) - +cos(solp/rd)*real(l-193))*4.91104/4.64947 xx=pmat(1)*x+pmat(3)*y yy=1.000112494*(pmat(2)*x+pmat(4)*y) ix=int(xx+192.5)+1 il=int(yy+192.5)+1 if (ix.gt.0 .and. ix.le.384 .and. il.gt.0 .and. - il.le.384) then ix=int(xx+192.5)+1 il=int(yy+192.5)+1 sum=0. swt=0. do 510 ll=il-1,il+1 do 510 ii=ix-1,ix+1 wt=exp(-((xx-real(ii-193))*(xx-real(ii-193)) - +(yy-real(ll-193))*(yy-real(ll-193))) - /0.5/0.5) swt=swt+wt sum=sum+wt*map(mod(ii+383,384)+1, - mod(ll+383,384)+1) 510 continue mapp(i,l)=sum/swt else mapp(i,l)=0. end if 500 continue end if if(ifreq.eq.34) then afact=(4.91104/2.)/(4.64947/33.8*17.) do 530 l=1,384 do 530 i=1,384 x=(cos(solp/rd)*real(i-193) - -sin(solp/rd)*real(l-193))*afact y=(sin(solp/rd)*real(i-193) - +cos(solp/rd)*real(l-193))*afact xx=pmat(1)*x+pmat(3)*y yy=1.000112494*(pmat(2)*x+pmat(4)*y) ix=int(xx+192.5)+1 il=int(yy+192.5)+1 if (ix.gt.0 .and. ix.le.384 .and. il.gt.0 .and. - il.le.384) then c rx=amod(xx+pxew*4.+193.+384.,384.) c ry=amod(yy+pxns*4.+193.+384.,384.) rx=amod(xx+193.+384.,384.) ry=amod(yy+193.+384.,384.) ix=int(rx+0.5) il=int(ry+0.5) sum=0. swt=0. do 540 ll=il-1,il+1 do 540 ii=ix-1,ix+1 c wt=exp(-((rx-real(ii-193))*(rx-real(ii-193)) c - +(ry-real(ll-193))*(ry-real(ll-193)))/0.5/0.5/16.) wt=exp(-((rx-real(ii))*(rx-real(ii)) - +(ry-real(ll))*(ry-real(ll)))/0.5/0.5) swt=swt+wt sum=sum+wt*map(mod(ii+383,384)+1, - mod(ll+383,384)+1) c write(6,*) i,l,rx,ry,xx,yy,ii,ll,wt, c - (rx-real(ii)),(ry-real(ll)) 540 continue ccc if (swt.lt.1.e-5) write(6,*) ix,il,swt,sum mapp(i,l)=sum/swt else mapp(i,l)=0. end if 530 continue end if c do 520 l=1,384 do 520 i=1,384 520 map(i,l)=mapp(i,l) c return end c c ******************************************************************** *********************************************************************** * cmap2dmap - projection correction for a snapshot map at 17 GHz * * solar north pole is to the top of the output image * * intensity is normalized to disk brightness = 10000. * * intensity information from nearby 3x3 pixels * * pixel size of the output map = 4.91104 x 4.91104 arcsec^2 * * (identical to the SXT half resolution images) * * pmat : input,real(4) projection matrix * * solp : input,real polar angle of the sun (degree) * * pxew : input,real x-offset of the disk center * * on the original map * * pxns : input,real l-offset of the disk center * * on the original map * *********************************************************************** c subroutine cmap2dmap(ifreq,pmat,solp,xd,yd) c real pmat(4) c,pmat1(4),fac integer ifreq c rd=57.29577951 pi=3.141592653 if (ifreq.eq.17) then x=(cos(solp/rd)*real(xd-257) - -sin(solp/rd)*real(yd-257))*4.91104/4.64947 y=(sin(solp/rd)*real(xd-257) - +cos(solp/rd)*real(yd-257))*4.91104/4.64947 xd=pmat(1)*x+pmat(3)*y+257. yd=1.000112494*(pmat(2)*x+pmat(4)*y)+257. end if if(ifreq.eq.34) then afact=(4.91104/2.)/(4.64947/33.8*17.) x=(cos(solp/rd)*real(xd-513) - -sin(solp/rd)*real(yd-513))*afact y=(sin(solp/rd)*real(xd-513) - +cos(solp/rd)*real(yd-513))*afact xd=pmat(1)*x+pmat(3)*y+513. yd=1.000112494*(pmat(2)*x+pmat(4)*y)+513. end if ***** super uniform weighting if (ifreq.eq.170) then xd=xd*2. yd=yd*2. x=(cos(solp/rd)*real(xd-513) - -sin(solp/rd)*real(yd-513))*4.91104/4.64947 y=(sin(solp/rd)*real(xd-513) - +cos(solp/rd)*real(yd-513))*4.91104/4.64947 xd=pmat(1)*x+pmat(3)*y+513. yd=1.000112494*(pmat(2)*x+pmat(4)*y)+513. end if c return end c *********************************************************************** * full_uv - create UV mask for use antennas * * * * nantst : input,integer(84) antenna status * * mask : output,real(512,512) UV mask * * March 5,1996 K. Fujiki * *********************************************************************** c subroutine full_uv(nantst,mask) integer nantst(84) real mask(512,512) integer ipos(-28:28) data ipos/-160,-144,-128,-112,-96,-80,-64, - -56,-48,-40,-32, - -28,-24,-20,-16,-14,-12, - -11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1, - 0,1,2,3,4,5,6,7,8,9,10,11, - 12,14,16,20,24,28, - 32,40,48,56, - 64,80,96,112,128,144,160/ c do 10 l=1,512 do 10 i=1,512 10 mask(i,l)=0.0 c do 20 l=-27,27 do 20 i=-28,28 20 mask(ipos(i)+257,ipos(l)+257)=1.0 c ***** 0-set for correlations of N-kei unuse antennas ***** do 300 iex=1,28 if (nantst(iex).ne.0) then do 400 i=1,512 mask(i,ipos(iex-1)+257)=0. 400 continue do 410 i=1,512 mask(i,-ipos(iex-1)+257)=0. 410 continue end if 300 continue c ***** 0-set for correlations of E-kei unuse antennas ***** do 310 iex=29,56 if (nantst(iex).ne.0) then do 420 l=258,512 mask(ipos(57-iex)+257,l)=0. 420 continue do 430 l=1,256 mask(-ipos(57-iex)+257,l)=0. 430 continue end if 310 continue c ***** 0-set for correlations of W-kei unuse antennas ***** do 320 iex=57,84 if (nantst(iex).ne.0) then do 440 l=257,512 mask(-ipos(iex-56)+257,l)=0. 440 continue do 450 l=1,257 mask(ipos(iex-56)+257,l)=0. 450 continue end if 320 continue c return end ************************************************************************* * jitter - jittercorrection * * nfreq : input,integer observation frequency * * pmat : input,real(4) projection matrix * * solr : input,real radius of the model sun (arcsec) * * nantst : input,integer*(84) antenna status * * pxew : input real x-offset of the center * * pxns : input,real y-offset of the center * * mapd : output,real(512,512) dirty map * * mapc : output,real(512,512) clean map * * mshis : output,integer(4) histogram of model sun * * October 18,1995 Y. Hanaoka * ************************************************************************* c subroutine fshift(mapd,mapa,pxew,pxns,dim) c integer dim real mapd(dim,dim),mapa(dim,dim) real pxew,pxns c rd=57.29577951 pi=3.141592654 ccc write(6,*) 'in jitter16 ',pxew,pxns c c ix0=dim/2+1 ***** position shift of model sun ***** c do 600 l=1,dim do 610 i=1,dim am=sqrt(mapd(i,l)*mapd(i,l)+mapa(i,l)*mapa(i,l)) if (am.lt.1.e-4) then ph=0. else ph=acos(max(min(mapd(i,l)/am,1.),-1.)) if (mapa(i,l).lt.0.) ph=2*pi-ph end if ph1=pxew*real(i-ix0)/real(dim)*360./rd ph2=pxns*real(l-ix0)/real(dim)*360./rd ph=ph+pxew*real(i-ix0)/real(dim)*360./rd - +pxns*real(l-ix0)/real(dim)*360./rd mapd(i,l)=am*cos(ph) mapa(i,l)=am*sin(ph) 610 continue c write(6,*) 'i,l,phx,phy= ',i,l,ph1,ph2 600 continue c return end c ************************************************************************* * jitter - jittercorrection * * nfreq : input,integer observation frequency * * pmat : input,real(4) projection matrix * * solr : input,real radius of the model sun (arcsec) * * nantst : input,integer*(84) antenna status * * pxew : input real x-offset of the center * * pxns : input,real y-offset of the center * * mapd : output,real(512,512) dirty map * * mapc : output,real(512,512) clean map * * mshis : output,integer(4) histogram of model sun * * October 18,1995 Y. Hanaoka * ************************************************************************* c subroutine jitter16(mapd,mapa,pxew,pxns) c real mapd(512,512),mapa(512,512) c rd=57.29577951 pi=3.141592654 ccc write(6,*) 'in jitter16 ',pxew,pxns c c ***** position shift of model sun ***** c do 600 l=1,512 do 610 i=1,512 am=sqrt(mapd(i,l)*mapd(i,l)+mapa(i,l)*mapa(i,l)) if (am.lt.1.e-4) then ph=0. else ph=acos(max(min(mapd(i,l)/am,1.),-1.)) if (mapa(i,l).lt.0.) ph=2*pi-ph end if ph1=pxew*real(i-257)/512.*360./rd ph2=pxns*real(l-257)/512.*360./rd ph=ph+pxew*real(i-257)/512.*360./rd - +pxns*real(l-257)/512.*360./rd mapd(i,l)=am*cos(ph) mapa(i,l)=am*sin(ph) 610 continue c write(6,*) 'i,l,phx,phy= ',i,l,ph1,ph2 600 continue c return end c c c *********************************************************************** * proj34 - projection correction for a snapshot map * * solar north pole is to the top of the output image * * intensity is normalized to disk brightness = 10000. * * intensity information from nearby 3x3 pixels * * pixel size of the output map = 4.91104 x 4.91104 arcsec^2 * * (identical to the SXT half resolution images) * * pmat : input,real(4) projection matrix * * solp : input,real polar angle of the sun (degree) * * dskbr : input,real abs(dskbr):disk brightness of * * the snapshot map * * if dskbr<0., main beam correction are done * * pxew : input,real x-offset of the disk center * * on the original map * * pxns : input,real l-offset of the disk center * * on the original map * * ioffp : input,integer x-offset of the center * * of the projected map from the disk center * * loffp : input,integer l-offset of the center * * of the projected map from the disk center * * map : input,real(512,512) snapshot map * * mapkk : output,real(1024,1024) projected map * * November 22,1995 Y. Hanaoka * *********************************************************************** c subroutine proj34_1(pmat,solp,dskbr,pxew,pxns,ioffp,loffp,map) c real map(1024,1024),mapp(1024,1024),pmat(4) c rd=57.29577951 pi=3.141592653 c ***** c afact=(4.91104/2.)/(4.64947/33.8*17.) x=(cos(solp/rd)*real(ioffp) - -sin(solp/rd)*real(loffp))*afact y=(sin(solp/rd)*real(ioffp) - +cos(solp/rd)*real(loffp))*afact pxx=pmat(1)*x+pmat(3)*y pyy=1.000112494*(pmat(2)*x+pmat(4)*y) do 500 l=1,1024 do 500 i=1,1024 x=(cos(solp/rd)*real(i-513) - -sin(solp/rd)*real(l-513))*afact y=(sin(solp/rd)*real(i-513) - +cos(solp/rd)*real(l-513))*afact xx=pmat(1)*x+pmat(3)*y yy=1.000112494*(pmat(2)*x+pmat(4)*y) ix=int(xx+512.5)+1 il=int(yy+512.5)+1 if (ix.gt.0 .and. ix.le.1024 .and. il.gt.0 .and. - il.le.1024) then rx=amod(xx+pxx+pxew+513.+1024.,1024.) ry=amod(yy+pyy+pxns+513.+1024.,1024.) ix=int(rx+0.5) il=int(ry+0.5) sum=0. swt=0. do 510 ll=il-1,il+1 do 510 ii=ix-1,ix+1 wt=exp(-((rx-real(ii))*(rx-real(ii))+(ry-real(ll)) - *(ry-real(ll)))/0.5/0.5) swt=swt+wt sum=sum+wt*map(mod(ii+1023,1024)+1, - mod(ll+1023,1024)+1) 510 continue if (swt.lt.1.e-5) write(6,'('' ix,il ='',2i5)') ix,il mapp(i,l)=sum/swt else mapp(i,l)=0. end if 500 continue c do 520 l=1,1024 do 520 i=1,1024 520 map(i,l)=mapp(i,l)/dskbr*10000. c return end ***************************************************************** * shift1 - position shift * * * * ioffm : input,real * * loffm : input,real * * map : input,real(512,512) * * shift_map: output,real(512,512) * * June 13,1995 koshix * ***************************************************************** subroutine shift1(ioffm,loffm,map,shift_map) integer u,v real ioffm,loffm real map,shift_map real map_real,map_imaginary real pai real cosine,sine real uv_cos,uv_sin,shift_uv_cos,shift_uv_sin dimension map(128,128) dimension shift_map(128,128) dimension map_real(128,128) dimension map_imaginary(128,128) dimension uv_cos(128,128) dimension uv_sin(128,128) dimension shift_uv_cos(128,128) dimension shift_uv_sin(128,128) parameter (pai=3.1415927) do 200 i=1,128 do 100 j=1,128 map_real(i,j)=map(i,j) map_imaginary(i,j)=0. 100 continue 200 continue call cft128(map_real,map_imaginary,1) do 400 i=1,128 do 300 j=1,128 uv_cos(i,j)=map_real(i,j) uv_sin(i,j)=map_imaginary(i,j) 300 continue 400 continue do 600 u=-64,63 do 500 v=-64,63 cosine=cos(2.*pai*(ioffm*real(u)+loffm*real(v))/128.) sine=-1.*sin(2.*pai*(ioffm*real(u)+loffm*real(v))/128.) shift_uv_cos(u+65,v+65)= & uv_cos(u+65,v+65)*cosine- & uv_sin(u+65,v+65)*sine shift_uv_sin(u+65,v+65)= & uv_cos(u+65,v+65)*sine+ & uv_sin(u+65,v+65)*cosine 500 continue 600 continue call cft128(shift_uv_cos,shift_uv_sin,-1) do 800 i=1,128 do 700 j=1,128 shift_map(i,j)=shift_uv_cos(i,j) 700 continue 800 continue return end *********************************************************************** * make_dtmp- create dirty map (512x512) * * * * map: input real array- original map * * rmap: oputput real array- dirty map * * mask :real input- masking array * * March 5,1996 K. Fujiki * *********************************************************************** c subroutine make_dtmp(map,rmap,mask) c real map(512,512),rmap(512,512),mask(512,512) real a(512,512),b(512,512) c do 100 l=1,512 do 100 i=1,512 a(i,l)=map(i,l) 100 continue c call cft512(a,b,1) c do 200 l=1,512 do 200 i=1,512 a(i,l)=a(i,l)*mask(i,l) b(i,l)=b(i,l)*mask(i,l) 200 continue c call cft512(a,b,-1) c do 300 l=1,512 do 300 i=1,512 rmap(i,l)=a(i,l) 300 continue c return end subroutine gauss_source(map,peak,width,x,y) real map(512,512),peak,width,x,y do 110 l=1,512 do 100 i=1,512 map(i,l)=0.0 100 continue 110 continue if(width.ne.0) then do 210 l=1,512 do 200 i=1,512 fin=real((i-x)*(i-x)+(l-y)*(l-y))/width/width if(fin.lt.50) then map(i,l)=peak*exp(-fin) else map(i,l)=0.0 end if 200 continue 210 continue else map(int(x),int(y))=peak end if return end *********************************************************************** * writefits- create dirty map (512x512) * * * * map: input real array- original map * * rmap: oputput real array- dirty map * * mask :real input- masking array * * March 5,1996 K. Fujiki * *********************************************************************** subroutine writefits(istat,ndut,msut,ndjst,msjst, - msjsts,msjste,iframe,itgr,natt,nfalt,nfreq,dattyp,cradd, - nmaxa,pxew,pxns,ioffp,loffp,mbcyn, - solr,solp,solb,dec,ha,az,alt,pmat, - cmapr,idim,outfft) real cmapr(idim,idim),pmat(4) integer idim,naxes(2) character*80 outfft character*9 dattyp c c nbit=-32 naxis=2 naxes(1)=idim naxes(2)=idim iounit=2 istat=0 write(6,*)'dim= ',idim write(6,*) outfft c c call ftinit(iounit,outfft,2880,istat) if (istat.ne.0) then write(6,*) 'File initialization failed!! istat=',istat go to 900 end if c c dattyp='CLEAN_MAP' call puthdr2(iounit,nbit,naxis,naxes,ndut,msut, - ndjst,msjst,msjsts,msjste,iframe,iframe+itgr-1, - 0,natt,nfalt,nfreq,dattyp) c c if (cradd.ge.0.) then cr=cradd else cr=abs(cradd*dskbr) endif c c c call putclp(iounit,nfreq,cr,nmaxa,pxew,pxns,ioffp,loffp,mbcyn) call putoep(iounit,solr,solp,solb,dec,ha,az,alt,pmat) c 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) then write(6,*)'Fail to write FITS file!! istat=',istat go to 900 end if write(6,'('' file '',a40,'' created'')') outfft c c 900 continue c c return end c *********************************************************************** * puthdr - 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 * * ndut : input,integer observation date in ut (yymmdd) * * msut : input,integer center of observation time * * in ut (hhmmsssss) * * ndjst : input,integer observation date in jst (yymmdd) * * msjst : input,integer center of observation time * * in jst (hhmmsssss) * * msjsts : input,integer start of observation time * * in jst (hhmmsssss) * * msjste : input,integer end of observation time * * in jst (hhmmsssss) * * 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,integer attenuator 0 or 10 (dB) * * nfalt : input,integer frequency alternation * * 0 : fix, 1 : alternate * * nfreq : input,integer observation frequency * * dattyp : input,character*18 data type * * October 27, 1995 Y. Hanaoka * *********************************************************************** c subroutine puthdr2(iunit,nbit,naxis,naxes,ndut,msut, - ndjst,msjst,msjsts,msjste,nfrst,nfrend, - npol,natt,nfalt,nfreq,dattyp) c integer naxes(2) character*18 dtobs,tmobs,jstdt,jsttm,jsts,jste, - pol,att,alt,freq,dattyp character*44 cmnt c cmnt=' ' c ***** write primary header information ***** call ftpprh(iunit,1,nbit,naxis,naxes,0,0,0,istat) c write(dtobs,'(i2.2,''/'',i2.2,''/'',i2.2)') mod(ndut,100), - mod(ndut,10000)/100,ndut/10000 c write(tmobs,'(i2.2,'':'',i2.2,'':'',f6.3)') msut/10000000, c - mod(msut,10000000)/100000,real(mod(msut,100000))/1000. write(tmobs,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msut/10000000,mod(msut,10000000)/100000, - mod(msut,100000)/1000,mod(msut,1000) write(jstdt,'(i2.2,''/'',i2.2,''/'',i2.2)') mod(ndjst,100), - mod(ndjst,10000)/100,ndjst/10000 c write(jsttm,'(i2.2,'':'',i2.2,'':'',f6.3)') msjst/10000000, c - mod(msjst,10000000)/100000,real(mod(msjst,100000))/1000. write(jsttm,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msjst/10000000,mod(msjst,10000000)/100000, - mod(msjst,100000)/1000,mod(msjst,1000) c write(jsts,'(i2.2,'':'',i2.2,'':'',f6.3)') msjsts/10000000, c - mod(msjsts,10000000)/100000,real(mod(msjsts,100000))/1000. write(jsts,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msjsts/10000000,mod(msjsts,10000000)/100000, - mod(msjsts,100000)/1000,mod(msjsts,1000) c write(jste,'(i2.2,'':'',i2.2,'':'',f6.3)') msjste/10000000, c - mod(msjste,10000000)/100000,real(mod(msjste,100000))/1000. write(jste,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msjste/10000000,mod(msjste,10000000)/100000, - mod(msjste,100000)/1000,mod(msjste,1000) 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 write(att,'(i2.2,''dB'')') natt if (nfalt.eq.0) then write(alt,'(''fix'')') else write(alt,'(''alt'')') endif write(freq,'(i2.2,''GHz'')') nfreq c call ftpkys(iunit,'date-obs',dtobs, cmnt,istat) call ftpkys(iunit,'time-obs',tmobs, 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) call ftpkys(iunit,'polariz', pol, cmnt,istat) call ftpkys(iunit,'att-10db',att, cmnt,istat) call ftpkys(iunit,'obs-mode',alt, cmnt,istat) call ftpkys(iunit,'obs-freq',freq, 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) cmnt='coded by K. Fujiki' call ftpkys(iunit,'pversion','6.0',cmnt,istat) c if (istat.ne.0) write(6,'('' puthdr : istat ='',i6)') istat return end c *********************************************************************** * puthdr - 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 * * ndut : input,integer observation date in ut (yymmdd) * * msut : input,integer center of observation time * * in ut (hhmmsssss) * * ndjst : input,integer observation date in jst (yymmdd) * * msjst : input,integer center of observation time * * in jst (hhmmsssss) * * msjsts : input,integer start of observation time * * in jst (hhmmsssss) * * msjste : input,integer end of observation time * * in jst (hhmmsssss) * * 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,integer attenuator 0 or 10 (dB) * * nfalt : input,integer frequency alternation * * 0 : fix, 1 : alternate * * nfreq : input,integer observation frequency * * dattyp : input,character*18 data type * * October 27, 1995 Y. Hanaoka * *********************************************************************** c subroutine puthdr_pfi(iunit,nbit,naxis,naxes,ndut,msut, - ndjst,msjst,msjsts,msjste,nfrst,nfrend, - npol,natt,nfalt,nfreq,dattyp,image_type) c integer naxes(2) character*18 dtobs,tmobs,jstdt,jsttm,jsts,jste, - pol,att,alt,freq,dattyp,aliasing character*44 cmnt c cmnt=' ' c ***** write primary header information ***** call ftpprh(iunit,1,nbit,naxis,naxes,0,0,0,istat) c c ndmy=ndut/10000 c write(6,*) 'ndmy= ',ndmy,ndut c if(ndmy.lt.69) then c ndmy=ndmy+2000 c else c ndmy=ndmy+1900 c endif c write(dtobs,'(i4.4,''-'',i2.2,''-'',i2.2)') ndmy, write(dtobs,'(i4.4,''-'',i2.2,''-'',i2.2)') ndut/10000, - mod(ndut,10000)/100,mod(ndut,100) c write(tmobs,'(i2.2,'':'',i2.2,'':'',f6.3)') msut/10000000, c - mod(msut,10000000)/100000,real(mod(msut,100000))/1000. write(tmobs,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msut/10000000,mod(msut,10000000)/100000, - mod(msut,100000)/1000,mod(msut,1000) c ndmy=ndjst/10000 c write(6,*) 'ndjst= ',ndmy,ndut c if(ndmy.lt.69) then c ndmy=ndmy+2000 c else c ndmy=ndmy+1900 c endif c write(jstdt,'(i4.4,''-'',i2.2,''-'',i2.2)') ndmy, write(jstdt,'(i4.4,''-'',i2.2,''-'',i2.2)') ndjst/10000, - mod(ndjst,10000)/100,mod(ndjst,100) c write(jsttm,'(i2.2,'':'',i2.2,'':'',f6.3)') msjst/10000000, c - mod(msjst,10000000)/100000,real(mod(msjst,100000))/1000. write(jsttm,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msjst/10000000,mod(msjst,10000000)/100000, - mod(msjst,100000)/1000,mod(msjst,1000) c write(jsts,'(i2.2,'':'',i2.2,'':'',f6.3)') msjsts/10000000, c - mod(msjsts,10000000)/100000,real(mod(msjsts,100000))/1000. write(jsts,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msjsts/10000000,mod(msjsts,10000000)/100000, - mod(msjsts,100000)/1000,mod(msjsts,1000) c write(jste,'(i2.2,'':'',i2.2,'':'',f6.3)') msjste/10000000, c - mod(msjste,10000000)/100000,real(mod(msjste,100000))/1000. write(jste,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msjste/10000000,mod(msjste,10000000)/100000, - mod(msjste,100000)/1000,mod(msjste,1000) c if (image_type.eq.1) then write(pol,'(''rcp'')') else if (image_type.eq.2) then write(pol,'(''lcp'')') else if (image_type.eq.3) then write(pol,'(''r+l'')') else if (image_type.eq.4) then write(pol,'(''r-l'')') else write(pol,'(''unknown'')') end if ccc write(6,*) 'alias ',i_alias if(i_alias.eq.1) then write(aliasing,'(''on'')') else write(aliasing,'(''off'')') end if write(att,'(i2.2,''dB'')') natt if (nfalt.eq.0) then write(alt,'(''fix'')') else write(alt,'(''alt'')') endif write(freq,'(i2.2,''GHz'')') nfreq c call ftpkys(iunit,'date-obs',dtobs, cmnt,istat) call ftpkys(iunit,'time-obs',tmobs, 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',att, cmnt,istat) call ftpkys(iunit,'obs-mode',alt, cmnt,istat) call ftpkys(iunit,'obs-freq',freq, 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) call ftpkys(iunit,'polariz',pol,cmnt,istat) call ftpkys(iunit,'hdrident','HeliogFITS 2.0',cmnt,istat) cmnt=' ' c if (istat.ne.0) write(6,'('' puthdr : istat ='',i6)') istat return end c *********************************************************************** * put_imparam - write FITS header for image information * * iunit * nfrcal * ical * itgr * pcritr * criter * mbeamc * icell * ddoff1 * ddoff2 * dskbr * *********************************************************************** subroutine put_imparam(iunit,nfrcal,ical,itgr,pcriter,criter, - solfac, - mbeamc,icell,ddoff1,ddoff2,dskbr) c real pcriter,criter,solfac,ddoff1,ddoff2,dskbr integer iunit,ical,nfrcal character*44 cmnt cmnt='calibration frame' call ftpkyj(iunit,'nfrcal',nfrcal,cmnt,istat) cmnt='integration for calibration' call ftpkyj(iunit,'ical',ical,cmnt,istat) cmnt='integration for restoration' call ftpkyj(iunit,'itgr',itgr,cmnt,istat) c cmnt='integration for after restoration' c call ftpkyj(iunit,'icln',icln,cmnt,istat) cmnt='CLEAN criterion for pre-CLEAN' call ftpkyf(iunit,'pcriter',pcriter,1,cmnt,istat) cmnt='CLEAN criterion for CLEAN' call ftpkyf(iunit,'criter',criter,1,cmnt,istat) cmnt='radius correction factor' call ftpkyf(iunit,'solr-fac',solfac,5,cmnt,istat) cmnt='x-offset of the dirty disk' call ftpkyf(iunit,'ddoff1',ddoff1,5,cmnt,istat) cmnt='y-offset of the dirty disk' call ftpkyf(iunit,'ddoff2',ddoff2,5,cmnt,istat) cmnt='dirty disk brightness' call ftpkyf(iunit,'dskbr',dskbr,1,cmnt,istat) cmnt='disk = 10000K' call ftpkys(iunit,'bunit','K',cmnt,istat) cmnt='cell size for weighting' call ftpkyj(iunit,'cellsize',icell,cmnt,istat) cmnt='main beam correction' call ftpkys(iunit,'mbeamc','yes',cmnt,istat) cmnt='1999-02-28 version coded by K. Fujiki' call ftpkys(iunit,'progname','snap2d17suw',cmnt,istat) end *********************************************************************** * puthdr - 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 * * ndut : input,integer observation date in ut (yymmdd) * * msut : input,integer center of observation time * * in ut (hhmmsssss) * * ndjst : input,integer observation date in jst (yymmdd) * * msjst : input,integer center of observation time * * in jst (hhmmsssss) * * msjsts : input,integer start of observation time * * in jst (hhmmsssss) * * msjste : input,integer end of observation time * * in jst (hhmmsssss) * * 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,integer attenuator 0 or 10 (dB) * * nfalt : input,integer frequency alternation * * 0 : fix, 1 : alternate * * nfreq : input,integer observation frequency * * dattyp : input,character*18 data type * * October 27, 1995 Y. Hanaoka * *********************************************************************** c subroutine puthdr_pfi_old(iunit,nbit,naxis,naxes,ndut,msut, - ndjst,msjst,msjsts,msjste,nfrst,nfrend, - npol,natt,nfalt,nfreq,dattyp,icell,image_type,i_alias) c integer naxes(2) character*18 dtobs,tmobs,jstdt,jsttm,jsts,jste, - pol,att,alt,freq,dattyp,aliasing,cell character*44 cmnt c cmnt=' ' c ***** write primary header information ***** call ftpprh(iunit,1,nbit,naxis,naxes,0,0,0,istat) c write(dtobs,'(i2.2,''/'',i2.2,''/'',i2.2)') mod(ndut,100), - mod(ndut,10000)/100,ndut/10000 c write(tmobs,'(i2.2,'':'',i2.2,'':'',f6.3)') msut/10000000, c - mod(msut,10000000)/100000,real(mod(msut,100000))/1000. write(tmobs,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msut/10000000,mod(msut,10000000)/100000, - mod(msut,100000)/1000,mod(msut,1000) write(jstdt,'(i2.2,''/'',i2.2,''/'',i2.2)') mod(ndjst,100), - mod(ndjst,10000)/100,ndjst/10000 c write(jsttm,'(i2.2,'':'',i2.2,'':'',f6.3)') msjst/10000000, c - mod(msjst,10000000)/100000,real(mod(msjst,100000))/1000. write(jsttm,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msjst/10000000,mod(msjst,10000000)/100000, - mod(msjst,100000)/1000,mod(msjst,1000) c write(jsts,'(i2.2,'':'',i2.2,'':'',f6.3)') msjsts/10000000, c - mod(msjsts,10000000)/100000,real(mod(msjsts,100000))/1000. write(jsts,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msjsts/10000000,mod(msjsts,10000000)/100000, - mod(msjsts,100000)/1000,mod(msjsts,1000) c write(jste,'(i2.2,'':'',i2.2,'':'',f6.3)') msjste/10000000, c - mod(msjste,10000000)/100000,real(mod(msjste,100000))/1000. write(jste,'(i2.2,'':'',i2.2,'':'',i2.2,''.'',i3.3)') - msjste/10000000,mod(msjste,10000000)/100000, - mod(msjste,100000)/1000,mod(msjste,1000) c if (image_type.eq.1) then write(pol,'(''rcp'')') else if (image_type.eq.2) then write(pol,'(''lcp'')') else if (image_type.eq.3) then write(pol,'(''acp'')') else if (image_type.eq.4) then write(pol,'(''scp'')') else write(pol,'(''unknown'')') end if ccc write(6,*) 'alias ',i_alias if(i_alias.eq.1) then write(aliasing,'(''on'')') else write(aliasing,'(''off'')') end if write(att,'(i2.2,''dB'')') natt if (nfalt.eq.0) then write(alt,'(''fix'')') else write(alt,'(''alt'')') endif write(freq,'(i2.2,''GHz'')') nfreq c call ftpkys(iunit,'date-obs',dtobs, cmnt,istat) call ftpkys(iunit,'time-obs',tmobs, 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',att, cmnt,istat) call ftpkys(iunit,'obs-mode',alt, cmnt,istat) call ftpkys(iunit,'obs-freq',freq, 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) call ftpkys(iunit,'polariz',pol,cmnt,istat) c cmnt='PFI center X' c call ftpkyf(iunit,'crpix1',xcen,2,cmnt,istat) c cmnt='PFI center Y' c call ftpkyf(iunit,'crpix1',ycen,2,cmnt,istat) cmnt='coded by K. Fujiki' call ftpkys(iunit,'pversion','6.0',cmnt,istat) write(cell,'(i3)') icell cmnt='cell size for super uniform' call ftpkys(iunit,'cellsize',cell,cmnt,istat) cmnt='aliasing reduction' call ftpkys(iunit,'alising',aliasing,cmnt,istat) cmnt=' ' c if (istat.ne.0) write(6,'('' puthdr : istat ='',i6)') istat return end c *********************************************************************** * putclp - write clean and projection parameters to fits file * * iunit : input,integer fortran unit number for fits file * * nfreq : input,integer observation frequency * * criter : input,real clean criterion * * ncmp : input,integer number of clean components * * pxew : input,real x-offset of the center of dirty disk * * pxns : input,real y-offset of the center of dirty disk * * ioffp : input,integer x-offset of the center of clean disk * * loffp : input,integer y-offset of the center of clean disk * * mbcyn : input,integer main beam correction 1=yes, 2=no * * August 18, 1994 Y. Hanaoka * *********************************************************************** c subroutine putclp_pfi(iunit,nfreq,criter,ncmp,pxew,pxns, - ioffp,loffp,isam,xcen,ycen) c character*44 cmnt c cmnt=' ' c c call ftpkys(iunit,'bunit','k','disk = 10000 K',istat) c call ftpkyf(iunit,'criter',criter,2,'clean criterion',istat) c call ftpkyj(iunit,'ncompo',ncmp, c - 'number of clean components',istat) c call ftpkyf(iunit,'ddoff1',pxew,2, c - 'x-offset of the dirty disk',istat) c call ftpkyf(iunit,'ddoff2',pxns,2, c - 'y-offset of the dirty disk',istat) c c call ftpkyf(iunit,'crval1',0.,2,'disk center',istat) c call ftpkyf(iunit,'crval2',0.,2,'disk center',istat) c cmnt='from sun center (deg.)' cmnt='arcsec' call ftpkyf(iunit,'crval1',(xcen-257)*4.911,2,cmnt,istat) c cmnt='from sun center (deg.)' cmnt='arcsec' call ftpkyf(iunit,'crval2',(ycen-257)*4.911,2,cmnt,istat) cmnt=' ' if (nfreq.eq.17) then cmnt=' ' call ftpkyf(iunit,'crpix1',64.5,2,cmnt,istat) cmnt=' ' call ftpkyf(iunit,'crpix2',64.5,2,cmnt,istat) cmnt=' ' call ftpkyf(iunit,'cdelt1',2.45552,5,'arcsec',istat) call ftpkyf(iunit,'cdelt2',2.45552,5,'arcsec',istat) else cmnt=' ' call ftpkyf(iunit,'crpix1',64.5,2,cmnt,istat) cmnt=' ' call ftpkyf(iunit,'crpix2',64.5,2,cmnt,istat) cmnt=' ' call ftpkyf(iunit,'cdelt1',1.22776,5,'arcsec',istat) call ftpkyf(iunit,'cdelt2',1.22776,5,'arcsec',istat) endif call ftpkys(iunit,'ctype1','solar-west',cmnt,istat) call ftpkys(iunit,'ctype2','solar-north',cmnt,istat) c if (istat.ne.0) write(6,'('' putclp : istat ='',i6)') istat return end c *********************************************************************** * putclp - write clean and projection parameters to fits file * * iunit : input,integer fortran unit number for fits file * * nfreq : input,integer observation frequency * * criter : input,real clean criterion * * ncmp : input,integer number of clean components * * pxew : input,real x-offset of the center of dirty disk * * pxns : input,real y-offset of the center of dirty disk * * ioffp : input,integer x-offset of the center of clean disk * * loffp : input,integer y-offset of the center of clean disk * * mbcyn : input,integer main beam correction 1=yes, 2=no * * August 18, 1994 Y. Hanaoka * *********************************************************************** c subroutine putclp_pfi_old(iunit,nfreq,criter,ncmp,pxew,pxns, - ioffp,loffp,isam,xcen,ycen) c character*44 cmnt c cmnt=' ' c c call ftpkys(iunit,'bunit','k','disk = 10000 K',istat) call ftpkyf(iunit,'criter',criter,2,'clean criterion',istat) c call ftpkyj(iunit,'ncompo',ncmp, c - 'number of clean components',istat) c call ftpkyf(iunit,'ddoff1',pxew,2, c - 'x-offset of the dirty disk',istat) c call ftpkyf(iunit,'ddoff2',pxns,2, c - 'y-offset of the dirty disk',istat) c c call ftpkyf(iunit,'crval1',0.,2,'disk center',istat) c call ftpkyf(iunit,'crval2',0.,2,'disk center',istat) cmnt='PFI center X' call ftpkyf(iunit,'pficenx',xcen,2,cmnt,istat) cmnt='PFI center Y' call ftpkyf(iunit,'pficeny',ycen,2,cmnt,istat) cmnt=' ' if (nfreq.eq.17) then cmnt='FFI center X' call ftpkyf(iunit,'crpix1',257.-ioffp,2,cmnt,istat) cmnt='FFI center Y' call ftpkyf(iunit,'crpix2',257.-loffp,2,cmnt,istat) cmnt=' ' call ftpkyf(iunit,'cdelt1',2.45552,5,'arcsec',istat) call ftpkyf(iunit,'cdelt2',2.45552,5,'arcsec',istat) else cmnt='FFI center X' call ftpkyf(iunit,'crpix1',513.-ioffp,2,cmnt,istat) cmnt='FFI center Y' call ftpkyf(iunit,'crpix2',513.-loffp,2,cmnt,istat) cmnt=' ' call ftpkyf(iunit,'cdelt1',1.22776,5,'arcsec',istat) call ftpkyf(iunit,'cdelt2',1.22776,5,'arcsec',istat) endif call ftpkys(iunit,'ctype1','solar-west',cmnt,istat) call ftpkys(iunit,'ctype2','solar-north',cmnt,istat) c if (istat.ne.0) write(6,'('' putclp : istat ='',i6)') istat return end c *********************************************************************** * putclp - write clean and projection parameters to fits file * * iunit : input,integer fortran unit number for fits file * * nfreq : input,integer observation frequency * * criter : input,real clean criterion * * ncmp : input,integer number of clean components * * pxew : input,real x-offset of the center of dirty disk * * pxns : input,real y-offset of the center of dirty disk * * ioffp : input,integer x-offset of the center of clean disk * * loffp : input,integer y-offset of the center of clean disk * * mbcyn : input,integer main beam correction 1=yes, 2=no * * August 18, 1994 Y. Hanaoka * *********************************************************************** c subroutine putclp_cal(iunit,icald) c integer*2 icald(336) character*44 cmnt character comm0 character*2 comm1 character*6 comm3 do 10 i=1,84 if(i.lt.29) then comm0='N' write(comm1,'(i2.2)') i else if(i.lt.57) then write(comm1,'(i2.2)') (i-28) comm0='W' else write(comm1,'(i2.2)') (i-56) comm0='E' end if cmnt='phase of '//comm0//comm1 comm3='ph-'//comm0//comm1 call ftpkyf(iunit,comm3,icald(i)/100.,2,cmnt,istat) 10 continue do 20 i=1,84 if(i.lt.29) then comm0='N' write(comm1,'(i2.2)') i else if(i.lt.57) then write(comm1,'(i2.2)') (i-28) comm0='W' else write(comm1,'(i2.2)') (i-56) comm0='E' end if cmnt='amplitude*100 of '//comm0//comm1 comm3='ap-'//comm0//comm1 call ftpkyf(iunit,comm3,icald(i+84)/100.,2,cmnt,istat) 20 continue cmnt=' ' c if (istat.ne.0) write(6,'('' putclp_cal : istat ='',i6)') istat return end c c *********************************************************************** * mask_super - make mask with super-uniform weighting * * mask_super: output - mask with super-uniform weighting * * mask: input - mask with natural weighting * * icell: input - cell size to calculate the weight * * Feb 1,1997 K. Fujiki * *********************************************************************** subroutine make_super(mask_super,mask,icell) c real mask_super(512,512),mask(512,512) integer icell integer ipos(-28:28) data ipos/-160,-144,-128,-112,-96,-80,-64, - -56,-48,-40,-32, - -28,-24,-20,-16,-14,-12, - -11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1, - 0,1,2,3,4,5,6,7,8,9,10,11, - 12,14,16,20,24,28, - 32,40,48,56, - 64,80,96,112,128,144,160/ c c do j=1,512 do i=1,512 mask_super(i,j)=mask(i,j) enddo enddo c is=icell/2 do 100 j=-27,27 do 100 i=-28,28 ii=ipos(i)+257 jj=ipos(j)+257 n=0 do 200 l=jj-is,jj+is do 200 k=ii-is,ii+is if(mask(k,l).eq.1) n=n+1 200 continue if(n.eq.0) then write(6,*) ipos(i),ipos(j),' is unused.' else mask_super(ii,jj)=mask_super(ii,jj)/real(n) end if 100 continue do j=1,512 do i=1,512 mask_super(i,j)=mask_super(i,j)*icell*icell enddo enddo c end *********************************************************************** * msrc_super - dirty beam of gaussian sources for SUW * * nantst : input,integer*(84) antenna status * * gbeam : output,real(512,512,0:10) dirty beam of * * gaussian sources * * width : 0(point source),1,...,10 * * July 13,1992 Y. Hanaoka * * Modified Jan 31,1997 K. Fujiki * *********************************************************************** c subroutine msrc_super(nantst,gbeam,cbeam,mask) c integer nantst(1) real gbeam(512,512,0:5),a(1024,1024),b(1024,1024), - cbeam(256,256,0:5),mask(512,512),c(1024,1024) c do 100 iws=0,5 c ***** clean source ***** c do j=1,1024 do i=1,1024 a(i,j)=0. b(i,j)=0. c(i,j)=0. enddo enddo do j=-256,255 do i=-256,255 c(513+i,513+j)=mask(i+257,j+257) enddo enddo c if (iws.eq.0) then a(513,513)=100. else wsr=real(iws) sum=0. do 410 l=1,1024 do 410 i=1,1024 a(i,l)=a(i,l)+100.*exp(-((513.-real(i))*(513.-real(i)) - +(513.-real(l))*(513.-real(l)))/wsr/wsr) sum=sum+a(i,l) 410 continue c do j=1,1024 do i=1,1024 a(i,j)=a(i,j)*100./sum enddo enddo c end if c do j=1,256 do i=1,256 cbeam(i,j,iws)=a(384+i,384+j) enddo enddo c ***** dirty map ***** c call cft1024(a,b,1) do j=1,1024 do i=1,1024 a(i,j)=a(i,j)*c(i,j)/mask(257,257) b(i,j)=b(i,j)*c(i,j)/mask(257,257) enddo enddo call cft1024(a,b,-1) c do 300 j=1,512 do 300 i=1,512 300 gbeam(i,j,iws)=a(i+256,j+256) c 100 continue c return end c c c c *********************************************************************** * fclean256 - CLEAN process for extended source * * model source : gaussian source and flat disk * * menu : input,integer * * 2 : disk/pointsource(+) CLEAN with restore * * 1 : pointsource(+) CLEAN with restore * * -1 : pointsource(+/-) CLEAN with restore * * -2 : disk/pointsource(+/-) CLEAN with restore * * map : input,real(512,512) dirty map * * ddisk : input,real(512,512) dirty map of model disk * * cdisk : input,real(512,512) CLEAN map of model disk * * gaine : input,real CLEAN loop gain * * criter : input,real CLEAN criterion * * nantst : input,integer(84) antenna status * * nmax : input,integer maximum number of CLEAN components * * bg : output,real background brightness * * dskbr : output,real disk brightness * * extend : output,real(512,512) CLEAN points * * for extended source * * cmap : output,real(512,512) CLEAN map * * April 24, 1996 arranged by fujiki * *********************************************************************** subroutine fclean256(menu,map,cmap,gbeam,cbeam1,nmax,crite,icell, - gain) c integer nmax,ne real scale,crite real map(256,256),cmap(256,256) real cbeam(-40:40,-40:40) real emap(256,256),remap(256,256),map_status(256,256) real point(256,256),extend(256,256),gbeam(512,512,0:5) real cbeam1(256,256,0:5) ne=0 npix=0 do j=1,256 do i=1,256 remap(i,j)=0. extend(i,j)=0. enddo enddo call max256(menu,map,amx,imx,lmx) c cr=abs(amx)*0.01 cr=crite if(cr .lt. abs(amx)*c_ratio) then cr=abs(amx)*c_ratio end if nnn=0 do 1100 while(.true.) ***** search maximum ***** nnn=nnn+1 call max256(menu,map,amx,imx,lmx) if(abs(amx).gt.cr .and. ne.lt.nmax) then call src_width256(map,imx,lmx,gbeam,iwid) c iwid=0 call cleanpf256(map,gbeam,amx,imx,lmx,iwid,emap,scale, - cbeam1,gain) npix=1 ne=ne+npix do j=1,256 do i=1,256 remap(i,j)=remap(i,j)+emap(i,j)*scale extend(i,j)=extend(i,j)+map_status(i,j) enddo enddo if(mod(nnn,20).eq.0 .or. nnn.eq.1) then c write(6,'('' '',i6,'': amx,imx,lmx,npix,scale ='', c - f14.2,2i5,i6,f7.2,i2)')ne,amx,imx,lmx,npix, c - scale,iwid end if else go to 1200 endif 1100 continue 1200 if(ne.ne.0) then ccc write(6,'('' ne ='',i6,'' : amx ='',f11.2)')ne,amx else ccc write(6,'('' ne ='',i6)') ne endif c nmax=ne c********************************************* c* restore process * c********************************************* do j=1,256 do i=1,256 cmap(i,j)=0. point(i,j)=0. enddo enddo c***** restore extended source ***** do j=1,256 do i=1,256 point(i,j)=point(i,j)+remap(i,j) enddo enddo ***** convolution by gaussian beam ***** c c cbeam(:,:)=0. c cbeam(0,0)=1.0 c sum=1. c goto 2500 sum=0.0 do 340 ll=-40,40 do 339 ii=-40,40 c xa=2.8494093 c ya=3.1584145 if(icell.eq.3) then xa=2.5 ya=2.5 else if(icell.eq.5) then xa=2.3 ya=2.3 else if(icell.eq.7) then xa=2. ya=2. else if(icell.eq.9) then xa=1.8 ya=1.8 end if fin=0.5*(real(ii)*real(ii)/xa/xa+real(ll)*real(ll)/ya/ya) if (fin.le.50) then cbeam(ii,ll)=exp(-fin) sum=sum+cbeam(ii,ll) else cbeam(ii,ll)=0.0 end if 339 continue 340 continue 2500 continue do 346 ll=-40,40 do 345 ii=-40,40 cbeam(ii,ll)=cbeam(ii,ll)/sum 345 continue 346 continue c do 2400 l=1,256 do 2300 i=1,256 do 2200 ll=-40,40 do 2100 ii=-40,40 ib=mod(i+ii+255,256)+1 lb=mod(l+ll+255,256)+1 cmap(i,l)=cmap(i,l)+point(ib,lb)*cbeam(ii,ll) c cmap(i,l)=cmap(i,l)+point(i,l) 2100 continue 2200 continue 2300 continue 2400 continue c***** add residual ***** do j=1,256 do i=1,256 cmap(i,j)=cmap(i,j)+map(i,j) enddo enddo return end *********************************************************************** * cleanpf256 - subtract a dirty beam component * * map : input,real(512,512) dirty map * * : output,real(512,512) dirty map, point source * * subtracted * * gbeam : input,real(512,512,0:10) dirty beam of gaussian * * sources * * amx : output,real maximum of map * * imx : output,integer x-coordinate of the maximum point * * lmx : output,integer y-coordinate of the maximum point * * iwid : output,integer width of subtracted point source * * ap : output,real intensity of subtracted point source * * July 17,1992 Y. Hanaoka * *********************************************************************** subroutine cleanpf256(map,gbeam,amx,imx,lmx,iwid,emap,scale, - cbeam,gain) c integer ib,lb real map(256,256),gbeam(512,512,0:5),emap(256,256),beam(256,256), - diff(256,256),cbeam(256,256,0:5),add(256,256) c c***** loop gain = 0.02 ***** ap=amx*gain do j=1,256 do i=1,256 emap(i,j)=0. enddo enddo c ii=257-imx jj=257-lmx do j=1,256 do i=1,256 ib=mod(i+129-imx+255,256)+1 lb=mod(j+129-lmx+255,256)+1 map(i,j)=map(i,j)-ap*gbeam(ii+i,jj+j,iwid)/gbeam(257,257,iwid) emap(i,j)=ap*cbeam(ib,lb,iwid)/gbeam(257,257,iwid) enddo enddo c$$$ do j=1,256 c$$$ do i=1,256 c$$$ ib=mod(i-imx+384,256)+1 c$$$ lb=mod(l-lmx+384,256)+1 c$$$ map(i,j)=map(i,j)-ap*gbeam(258-imx,513-lmx,iwid) c$$$ - /gbeam(257,257,iwid) c$$$ emap(i,j)=emap(i,j)+ap*cbeam(ib,lb,iwid)/gbeam(257,257,iwid) c$$$ enddo c$$$ enddo scale=1.0 return end c$$$*********************************************************************** c$$$* cleanp - subtract a dirty beam component * c$$$* map : input,real(512,512) dirty map * c$$$* : output,real(512,512) dirty map, point source * c$$$* subtracted * c$$$* gbeam : input,real(512,512,0:10) dirty beam of gaussian * c$$$* sources * c$$$* amx : output,real maximum of map * c$$$* imx : output,integer x-coordinate of the maximum point * c$$$* lmx : output,integer y-coordinate of the maximum point * c$$$* iwid : output,integer width of subtracted point source * c$$$* ap : output,real intensity of subtracted point source * c$$$* July 17,1992 Y. Hanaoka * c$$$*********************************************************************** c$$$ subroutine cleanpf256(map,gbeam,amx,imx,lmx,iwid,emap,scale, c$$$ - cbeam,gain) c$$$c c$$$ real map(256,256),gbeam(512,512,0:5),emap(256,256), c$$$ - diff(256,256),cbeam(256,256,0:5),add(256,256) c$$$c c$$$c***** loop gain = 0.02 ***** c$$$ ap=amx*gain c$$$ emap=0. c$$$c c$$$ ii=257-imx+1 c$$$ jj=257-lmx+1 c$$$ diff=ap*gbeam(ii:ii+255,jj:jj+255,iwid)/gbeam(257,257,iwid) c$$$ add=ap*cbeam(:,:,iwid)/gbeam(257,257,iwid) c$$$ add=cshift(add,129-imx,1) c$$$ add=cshift(add,129-lmx,2) c$$$ c$$$ map=map-diff c$$$ c$$$ emap=emap+add c$$$ scale=1.0 c$$$ c$$$ return c$$$ end ***************************************************************** * src_width256 : estimate source width for 256x256 image * map- input,real dirty map * imx- input,integer source position * lmx- input,integer source position * gbeam- input,real dirty beam * iwid- output,integer source width * Feb. 28, 1999 K. Fujiki ***************************************************************** subroutine src_width256(map,imx,lmx,gbeam,iwid) real map(256,256),gbeam(512,512,0:5) do 20 iw=1,5 do 20 j=-5,5 do 20 i=-5,5 c imx1=mod(imx+i*2+255,256)+1 lmx1=mod(lmx+j*2+255,256)+1 if (abs(map(imx1,lmx1)) .lt. - gbeam(257+i*2,257+j*2,iw)/gbeam(257,257,iw) - *abs(map(imx,lmx))) then goto 29 endif 20 continue 29 iwid=iw-1 return end *********************************************************************** * proj16_1 - projection correction for a snapshot map * * solar north pole is to the top of the output image * * intensity is normalized to disk brightness = 10000. * * intensity information from nearby 3x3 pixels * * pixel size of the output map = 4.91104 x 4.91104 arcsec^2 * * (identical to the SXT half resolution images) * * pmat : input,real(4) projection matrix * * solp : input,real polar angle of the sun (degree) * * dskbr : input,real abs(dskbr):disk brightness of * * the snapshot map * * if dskbr<0., main beam correction are done * * ioff : input,integer x-offset of the disk center * * on the original map * * loff : input,integer l-offset of the disk center * * on the original map * * ioffp : input,integer x-offset of the center * * of the projected map from the disk center * * loffp : input,integer l-offset of the center * * of the projected map from the disk center * * map : input,real(512,512) snapshot map * * output,real(512,512) projected map * * April 29,1993 Y. Hanaoka * *********************************************************************** c subroutine proj_super(ifreq,pmat,solp,map) c real map(256,256),mapp(256,256),pmat(4) integer ifreq c rd=57.29577951 c write(6,*) 'start projection',solp,pmat c cc do j=1,256 do i=1,256 mapp(i,j)=0.0 enddo enddo do 500 l=1,256 do 500 i=1,256 x=(cos(solp/rd)*real(i-129.) - -sin(solp/rd)*real(l-129.))*4.91104/4.64947 y=(sin(solp/rd)*real(i-129.) - +cos(solp/rd)*real(l-129.))*4.91104/4.64947 xx=pmat(1)*x+pmat(3)*y yy=1.000112494*(pmat(2)*x+pmat(4)*y) ix=int(xx+128.5)+1 il=int(yy+128.5)+1 if (ix.gt.0 .and. ix.le.256 .and. il.gt.0 .and. - il.le.256) then ix=int(xx+128.5)+1 il=int(yy+128.5)+1 sum=0. swt=0. do 510 ll=il-1,il+1 do 510 ii=ix-1,ix+1 wt=exp(-((xx-real(ii-129.))*(xx-real(ii-129.)) - +(yy-real(ll-129.))*(yy-real(ll-129.)))/0.5/0.5) swt=swt+wt sum=sum+wt*map(mod(ii+255,256)+1, - mod(ll+255,256)+1) 510 continue mapp(i,l)=sum/swt else mapp(i,l)=0. end if 500 continue c write(6,*) 'exit projection' do 520 l=1,256 do 520 i=1,256 520 map(i,l)=mapp(i,l) c return end c c *********************************************************************** * pre_clean - First step CLEAN * menu * mapra input,real dirty map * output,real CLEAN map * maprb output,real residual * xmc0 input,real center of PFI * ymc0 input,real center of PFI * cradd input,real CLEAN criterion * nmax input,real maximum CLEAN component * mask_all input,real sampling function of UV components * mask input,real mask for image plane * pmat input,real projection matrix * solp input,real position angle * iframe input,integer frame number * outf input,character file name header * gain input,real gain factor * gbeam input,real dirty beam * taper input,real taper at xd (see below) * xd input,real xd harmonics * istr * Feb. 28, 1999 K. Fujiki *********************************************************************** subroutine pre_clean(menu,mapra,maprb,xmc0,ymc0,cradd, - nmax,mask_all,mask,pmat,solp,iframe,outf,gain,gbeam, - taper,xd,istr) real mapra(512,512),mask_all(512,512),maprb(512,512), - maskc(512,512),pmat(4),aa(512,512),gbeam(512,512,0:10), - bb(512,512),dd(512,512),taper,xd, - mask(512,512),solp character*80 outf,outfft character*5 fnum ccc write(6,*)'Start Pre-CLEAN...' c write(6,*)'Target ',xmc0,ymc0 ccc write(6,*)'CLEAN Criterion and Limit ',cradd,nmax ccc write(6,*)'TAPER ',taper,xd write(fnum,'(i5.5)') iframe loutf=index(outf,' ')-1 do j=1,512 do i=1,512 maprb(i,j)=mapra(i,j) enddo enddo ***** calculate PFI image center in dirty map xmc1=xmc0 ymc1=ymc0 call cmap2dmap(17,pmat,solp,xmc1,ymc1) c write(6,*)'image center= ',xmc1,ymc1 call make_mask(1,mask,maskc,xmc1,ymc1) call max512(-1,maprb,amx0,imx,lmx) do j=1,512 do i=1,512 bb(i,j)=maprb(i,j)*maskc(i,j) enddo enddo call max512(-1,bb,amx1,imx,lmx) ccc write(6,*)'peak brightness of full dirty map =',amx0,amx0/200. ccc write(6,*)'peak brightness of out of target =',amx1 if(abs(amx1).lt.abs(amx0/200.)) then write(6,*)'skip pre CLEAN process' go to 100 end if ***** CLEAN by Steer argolithm call fclean512(menu,maprb,gain,cradd,mask_all,nmax,aa,gbeam, - istr) ***** maprb : residual do j=1,512 do i=1,512 bb(i,j)=maprb(i,j) enddo enddo if(taper.eq.0.0 .and. xd.eq.0.0) then do j=1,512 do i=1,512 bb(i,j)=0. dd(i,j)=0. enddo enddo else if(taper.lt.1.) then ***** bb : residual with taper (taper at xd d) ***** dd : residual without real structure call cft512(bb,dd,1) call filter4(taper,xd,bb) call filter4(taper,xd,dd) call cft512(bb,dd,-1) ***** mesh pattern do j=1,512 do i=1,512 dd(i,j)=maprb(i,j)-bb(i,j) enddo enddo else do j=1,512 do i=1,512 dd(i,j)=0. enddo enddo end if end if ***** Target image aa:CLEAN map, bb:residual c c outfft='in_pcrn_cmap.fits' c call wt_fits(512,aa,outfft) do j=1,512 do i=1,512 maprb(i,j) = aa(i,j) * abs(maskc(i,j)-1.) + bb(i,j) enddo enddo c c outfft='in_pcrn_cmap_mask.fits' c call wt_fits(512,maprb,outfft) c outfft='in_mask.fits' c call wt_fits(512,maskc,outfft) call cft512(maprb,aa,1) do j=1,512 do i=1,512 maprb(i,j)=maprb(i,j)*mask_all(i,j) aa(i,j)=aa(i,j)*mask_all(i,j) enddo enddo call cft512(maprb,aa,-1) c outfft='in_pcrn_cmap_mask_beam.fits' c call wt_fits(512,maprb,outfft) c do j=1,512 do i=1,512 mapra(i,j)=maprb(i,j) maprb(i,j)=bb(i,j) enddo enddo c c mapra=maprb c maprb=bb c outfft='in_pcrn_ra.fits' c call wt_fits(512,mapra,outfft) c outfft='in_pcrn_rb.fits' c call wt_fits(512,maprb,outfft) 100 continue return end ************************************************************************* * make_mask - make mask to pick up CLEAN component of target region * isam - input,integer * mask - input,real input mask * rmask - output,real output mask * xx - real,input center of PFI * yy - real,input center of PFI * Feb. 28 1999 K. Fujiki ************************************************************************* subroutine make_mask(isam,mask,rmask,xx,yy) integer ix,iy,isam real rmask(512,512),aa(512,512), - mask(512,512),cbeam(512,512),xx,yy do j=1,512 do i=1,512 rmask(i,j)=0.0 enddo enddo if (isam.lt.0) then ***** mask size of super uniform weighting ii=64 do 100 j=1,512 do 100 i=1,512 if(i.ge.257-ii .and. i.le.257+ii-1 .and. - j.ge.257-ii .and. j.le.257+ii-1 ) then mask(i,j)=0.0 else mask(i,j)=1.0 end if aa(i,j)=0. 100 continue sum=0 call filter3(0.1,250.,cbeam) call cft512(mask,rmask,1) do j=1,512 do i=1,512 mask(i,j)=mask(i,j)*cbeam(i,j) rmask(i,j)=0. enddo enddo call cft512(mask,rmask,-1) do j=1,512 do i=1,512 rmask(i,j)=mask(i,j) enddo enddo else ix=int(257-xx+0.5) iy=int(257-yy+0.5) c aa=cshift(mask,ix,1) c rmask=cshift(aa,iy,2) do j=1,512 do i=1,512 ib=mod(i-ix+511,512)+1 lb=mod(j-iy+511,512)+1 rmask(i,j)=mask(ib,lb) enddo enddo end if return end *********************************************************************** * clean17f - clean process for 17 GHz images * * model source : gaussian source and flat disk * * menu : input,integer * * -1 : pointsource(+/-) clean with restore * * 1 : pointsource(+) clean with restore * * map : input,real(512,512) dirty map * * output,real(512,512) resudual * * gbeam : input,real(512,512,0:10) dirty beam of gaussian * * sources * * criter : input,real clean criterion * * if < 0, cr=-criter*dskbr * * nmax : input,integer maximum number of clean components * * cmap : output,real(512,512) clean map * * March 19,1996 K. Fujiki * *********************************************************************** c subroutine clean17f(menu,map,criter,nmax,dskbr,gbeam,cmap) c integer pi(10000),pl(10000),pw(10000) real map(512,512),cmap(512,512),ddisk(512,512), - cbeam(-30:30,-30:30),gbeam(512,512,0:10), - pa(10000),sum c n=0 c ********************************************* * clean process * ********************************************* c ***** clean criterion ***** c if (criter.lt.0.) then cr=dskbr*(-criter) c write(6,'('' clean criterion :'',f10.1)') cr else cr=criter endif c do 160 while(.true.) c ***** search maximum ***** call mean(menu,map,am,amx,imx,lmx) if (abs(amx).gt.cr .and. n.lt.nmax) then c ***** subtract point source component ***** n=n+1 pi(n)=imx pl(n)=lmx call cleanpf(map,gbeam,amx,imx,lmx,pw(n),pa(n),emap) if (mod(n,10).eq.0) then c write(6,'('' n ='',i6,'': am,amx,w,a,imx,lmx ='', c - f10.2,f14.2,i3,f10.2,2i5)') c - n,am,amx,pw(n),pa(n),pi(n),pl(n) end if c else go to 170 end if c 160 continue 170 if (n.ne.0) then c write(6,'('' n ='',i6,'': am,amx,w,a,imx,lmx ='', c - f10.2,f14.2,i3,f10.2,2i5)') c - n,am,amx,pw(n),pa(n),pi(n),pl(n) else c write(6,'('' n ='',i6)') n end if c nmax=n if (abs(menu)/10.eq.1) return c c ********************************************* * restore process * ********************************************* c do 240 l=1,512 do 240 i=1,512 240 ddisk(i,l)=0. c do 500 iw0=0,1 iw=iw0 c if(iw.gt.10) go to 500 c ***** create clean beam ***** c if (iw.eq.0) then do 210 l=-30,30 do 210 i=-30,30 210 cbeam(i,l)=0. cbeam(0,0)=100. sum=100. else sum=0. do 220 l=-30,30 do 220 i=-30,30 fin=(real(i)*real(i)+real(l)*real(l)) - /real(iw)/real(iw) if(fin.le.40 .and. fin.gt.0) then cbeam(i,l)=100.*exp(-fin) else cbeam(i,l)=0.0 end if sum=sum+cbeam(i,l) 220 continue c do 230 l=-30,30 do 230 i=-30,30 cbeam(i,l)=100.*cbeam(i,l)/sum 230 continue end if cc c***** restore point sources ***** cc do 250 ip=1,n c if(pw(ip).ne.0) write(6,*)'pw= ',pw(ip),ip,n if(pw(ip).eq.iw) then do 260 l=-30,30 do 260 i=-30,30 ii=mod(pi(ip)+i+511,512)+1 ll=mod(pl(ip)+l+511,512)+1 ddisk(ii,ll)=ddisk(ii,ll)+pa(ip) - *cbeam(i,l)/gbeam(257,257,iw) 260 continue cc write(6,'('' '',i4,5f10.3)') iw,(cbeam(ii,0),ii=-10,10,5) end if 250 continue 500 continue c ***** convolution by gaussian beam ***** c c$$$ go to 349 c$$$c c$$$ do 340 ll=-10,10 c$$$ do 340 ii=-10,10 c$$$ 340 cbeam(ii,ll)=0.07461*exp(-0.234403 c$$$ - *(real(ll)*real(ll)+real(ii)*real(ii))) c$$$c c$$$ do 310 l=1,512 c$$$ do 310 i=1,512 c$$$ 310 cmap(i,l)=0. c$$$ do 320 l=1,512 c$$$ do 320 i=1,512 c$$$ do 330 ll=-10,10 c$$$ do 330 ii=-10,10 c$$$ ib=mod(i+ii+511,512)+1 c$$$ lb=mod(l+ll+511,512)+1 c$$$ cmap(i,l)=cmap(i,l)+ddisk(ib,lb)*cbeam(ii,ll) c$$$ 330 continue c$$$ 320 continue c$$$c c$$$c ***** add residual ***** c$$$ 349 continue c do 350 l=1,512 do 350 i=1,512 cmap(i,l)=ddisk(i,l)+map(i,l) 350 continue c return end c *********************************************************************** * cleanpf - subtract a dirty beam component * * map : input,real(512,512) dirty map * * : output,real(512,512) dirty map, point source * * subtracted * * gbeam : input,real(512,512,0:10) dirty beam of gaussian * * sources * * amx : output,real maximum of map * * imx : output,integer x-coordinate of the maximum point * * lmx : output,integer y-coordinate of the maximum point * * iwid : output,integer width of subtracted point source * * emap : output,integer clean component map * * Feb. 28 1999 K. Fujiki * *********************************************************************** c subroutine cleanpf(map,gbeam,amx,imx,lmx,iwid,emap,scale) c real map(512,512),gbeam(512,512,0:10),emap(512,512) c ***** loop gain = 0.02 ***** ap=amx*0.02 do j=1,512 do i=1,512 emap(i,j)=0. enddo enddo c do 300 l=1,512 do 300 i=1,512 ib=mod(i-imx+768,512)+1 lb=mod(l-lmx+768,512)+1 300 map(i,l)=map(i,l)-ap*gbeam(ib,lb,iwid)/gbeam(257,257,iwid) c do 400 il=1,512,64 c 400 write(6,'('' '',8f9.2)') (map(i,il),i=1,512,64) c emap(imx,lmx)=emap(imx,lmx)+ap scale=100./gbeam(257,257,iwid) return end ***************************************************************** * filter3 - high cut filtering * * * * taper : input,real taper at 160d * * uv : input/output,real(512,512) * * Jun 13,1995 koshix * * Modified Sep 30,1996 Fujiki * ***************************************************************** subroutine filter3(taper,xd,uv) real taper,sigma,r,xd real uv(512,512) if(taper .ne. 1.) then sigma=xd/sqrt(alog(1./taper)) do 200 i=1,512 do 100 j=1,512 xx=real(i-257) yy=real(j-257) r=sqrt(xx*xx+yy*yy) fin=(r*r)/(sigma*sigma) if(fin.lt.80) then uv(i,j)=exp(-1.*fin) else uv(i,j)=0.0 end if 100 continue 200 continue else do j=1,512 do i=1,512 uv(i,j)=1.0 enddo enddo end if return end ***************************************************************** * filter4 - high cut filtering * * * * taper : input,real taper at 160d * * uv : input/output,real(512,512) * * Jun 13,1995 koshix * * Modified Sep 30,1996 Fujiki * ***************************************************************** subroutine filter4(taper,xd,uv) real taper,sigma,r,xd real uv(512,512) if(taper .ne. 1.) then sigma=xd/sqrt(alog(1./taper)) do 200 i=1,512 do 100 j=1,512 xx=real(i-257) yy=real(j-257) r=sqrt(xx*xx+yy*yy) fin=(r*r)/(sigma*sigma) if(fin.lt.80) then uv(i,j)=uv(i,j)*exp(-1.*fin) else uv(i,j)=0.0 end if 100 continue 200 continue end if return end *********************************************************************** * fclean512 - CLEAN process for extended source * * model source : gaussian source and flat disk * * menu : input,integer * * 2 : disk/pointsource(+) CLEAN with restore * * 1 : pointsource(+) CLEAN with restore * * -1 : pointsource(+/-) CLEAN with restore * * -2 : disk/pointsource(+/-) CLEAN with restore * * map : input,real(512,512) dirty map * * ddisk : input,real(512,512) dirty map of model disk * * cdisk : input,real(512,512) CLEAN map of model disk * * gaine : input,real CLEAN loop gain * * criter : input,real CLEAN criterion * * nantst : input,integer(84) antenna status * * nmax : input,integer maximum number of CLEAN components * * bg : output,real background brightness * * dskbr : output,real disk brightness * * extend : output,real(512,512) CLEAN points * * for extended source * * cmap : output,real(512,512) CLEAN map * * June 13,1995 arranged by koshix * *********************************************************************** c subroutine fclean512(menu,map,gaine,crite,mask,nmax,cmap,gbeam, - istr) integer nmax,ne,istr real scale,crite real map(512,512),cmap(512,512),mask(512,512) real gbeam(512,512,0:3) real emap(512,512),remap(512,512),map_status(512,512) real extend(512,512) ccc write(6,*)'Criterion= ',crite ccc write(6,*)'Max Components= ',nmax ne=0 do j=1,512 do i=1,512 remap(i,j)=0. extend(i,j)=0. cmap(i,j)=0. enddo enddo gaine=0.2 c call mean(menu,map,am,amx,imx,lmx) call max512(menu,map,amx,imx,lmx) cr=crite nnn=0 do 1100 while(.true.) ***** search maximum ***** nnn=nnn+1 call max512(menu,map,amx,imx,lmx) if(abs(amx).gt.cr .and. ne.lt.nmax) then call src_width(map,imx,lmx,gbeam,iwid) if((iwid.lt.2 .and. abs(amx).gt.30000.) - .or. istr.eq.0) then c if(.false.) then c if(.true.) then iwid=0 gaine=0.02 npix=1 call cleanpf(map,gbeam,amx,imx,lmx,iwid,emap,scale) ne=ne+1 do j=1,512 do i=1,512 remap(i,j)=remap(i,j)+emap(i,j)*scale enddo enddo else trim=0. if(iwid.eq.0) then gaine=0.08 else if(iwid.eq.2) then gaine=0.1 else if(iwid.eq.3) then gaine=0.15 else gaine=0.2 end if c gaine=0.2 call contour_trim512(map,gaine,trim,npix,amx, - map_status,emap) ne=ne+npix call cleane2_512(map,mask,amx,gaine,emap,scale) do j=1,512 do i=1,512 remap(i,j)=remap(i,j)+emap(i,j)*scale enddo enddo end if if(mod(nnn,200).eq.0.or.iwid.gt.1.or.nnn.eq.1) then ccc write(6,'('' '',2i7,'': amx,npix,scale ='', ccc - f14.2,i6,f7.2,f7.2,i2)')nnn,ne,amx, ccc - npix,scale,gaine,iwid end if else go to 1200 endif 1100 continue 1200 if(ne.ne.0) then ccc write(6,'('' ne ='',i6,'' : amx ='',f11.2)')ne,amx else ccc write(6,'('' ne ='',i6)') ne endif c nmax=ne ********************************************* * restore process * ********************************************* ***** restore extended source ***** do j=1,512 do i=1,512 cmap(i,j)=remap(i,j) enddo enddo return end *********************************************************************** * cleane2_512 - subtract extended source Ver.2 * * map : input,real(512,512) dirty map * * : output,real(512,512) dirty map, extended source * * subtracted * * nantst : input,integer(84) antenna status * * amx : input,real maximum of map * * gaine : input,real CLEAN loop gain for extended source * * emap : input,real(512,512) extended source * * scale : output,real scaling factor for extended source * * June 13,1995 koshix * *********************************************************************** subroutine cleane2_512(map,mask,amx,gaine,emap,scale) real gaine,amx,dmx,scale real map(512,512),emap(512,512),mask(512,512) real a(512,512),b(512,512) do j=1,512 do i=1,512 a(i,j)=emap(i,j) b(i,j)=0. enddo enddo call cft512(a,b,1) do j=1,512 do i=1,512 a(i,j)=a(i,j)*mask(i,j) b(i,j)=b(i,j)*mask(i,j) enddo enddo c a(257,257)=0.0 c b(257,257)=0.0 call cft512(a,b,-1) dmx=0. do 400 j=1,512 do 300 i=1,512 if(abs(a(i,j)) .gt. abs(dmx)) dmx=a(i,j) 300 continue 400 continue do j=1,512 do i=1,512 map(i,j)=map(i,j)-a(i,j)/abs(dmx)*abs(amx)*gaine enddo enddo c sum1=0. c sum2=0. c c do 800 i=1,512 c do 700 j=1,512 c c sum1=sum1+emap(i,j) c sum2=sum2+a(i,j)/abs(dmx)*abs(amx)*gaine c c700 continue c800 continue c c scale=abs(sum2/sum1) scale=abs(amx)/abs(dmx)*gaine return end ***************************************************************** * contour_trim - find contour trim * * * * map : input,real(128,128) dirty map * * gaine : input,real CLEAN loop gain * * for extended source * * npix : output,integer number of pixels(1.-gaine) * * 2 : extended source * * map_status : output,integer(128,128) * * CLEAN point for extended source * * e_map : output,real(128,128) extended component * * March 30 ,1995 Fujiki * ***************************************************************** subroutine contour_trim512(map,gaine,trim,npix,amx, - map_status,emap) integer i,j integer npix real gaine real map(512,512) real emap(512,512),map_status(512,512) ********** Contour Trim(1-gaine) of Extended Source ********** if(trim.eq.0) trim=1.-gaine c write(6,*)'trim level= ',trim do j=1,512 do i=1,512 emap(i,j)=0. map_status(i,j)=0. emap(i,j)=map(i,j)/amx-trim enddo enddo do 300 j=1,512 do 400 i=1,512 c emap(i,j)=map(i,j)/amx-trim map_status(i,j)=real((int(sign(1.,emap(i,j)))+1)/2) 400 continue 300 continue do j=1,512 do i=1,512 emap(i,j)=map(i,j)*map_status(i,j)*gaine enddo enddo npix=0 do 900 j=1,512 do 800 i=1,512 npix=npix+int(map_status(i,j)) 800 continue 900 continue return end c$$$***************************************************************** c$$$* contour_trim - find contour trim * c$$$* * c$$$* map : input,real(512,512) dirty map * c$$$* gaine : input,real CLEAN loop gain * c$$$* for extended source * c$$$* imx : input,integer * c$$$* x-coordinate of the maximum point * c$$$* lmx : input,integer * c$$$* y-coordinate of the maximum point * c$$$* npix : output,integer number of pixels(1.-gaine) * c$$$* 2 : extended source * c$$$* map_status : output,integer(512,512) * c$$$* CLEAN point for extended source * c$$$* e_map : output,real(512,512) extended component * c$$$* June 13,1995 koshix * c$$$***************************************************************** c$$$ c$$$ subroutine contour_trim5121(map,gaine,imx,lmx,npix, c$$$ & map_status,emap) c$$$ c$$$ integer key c$$$ integer i,j,k,l,imx,lmx,ib,lb c$$$ integer npix c$$$ integer n,buffer_num c$$$ integer ibuffer(262144),jbuffer(262144) c$$$ integer map_status(512,512) c$$$ real gaine c$$$ real map(512,512) c$$$ real emap(512,512) c$$$ c$$$ c$$$********** Contour Trim(1-gaine) of Extended Source ********** c$$$ c$$$ do 200 i=1,512 c$$$ do 100 j=1,512 c$$$ c$$$ map_status(i,j)=0 c$$$ emap(i,j)=0. c$$$ c$$$100 continue c$$$200 continue c$$$ c$$$ do 300 n=1,262144 c$$$ c$$$ ibuffer(n)=-1 c$$$ jbuffer(n)=-1 c$$$ c$$$300 continue c$$$ c$$$ map_status(imx,lmx)=1 c$$$ buffer_num=1 c$$$ ibuffer(buffer_num)=imx c$$$ jbuffer(buffer_num)=lmx c$$$ c$$$ key=1 c$$$ c$$$ do 700 while(key .eq. 1) c$$$ c$$$ key=0 c$$$ c$$$ do 600 n=1,buffer_num c$$$ c$$$ do 500 k=-1,1 c$$$ do 400 l=-1,1 c$$$ c$$$ ib=mod(ibuffer(n)+k+511,512)+1 c$$$ lb=mod(jbuffer(n)+l+511,512)+1 c$$$ c$$$ if(abs(map(ib,lb)) .ge. c$$$ & abs(map(imx,lmx))*(1.-gaine) .and. c$$$ & map_status(ib,lb) .ne. 1) then c$$$ c$$$ map_status(ib,lb)=1 c$$$ buffer_num=buffer_num+1 c$$$ ibuffer(buffer_num)=ib c$$$ jbuffer(buffer_num)=lb c$$$ c$$$ key=1 c$$$ c$$$ endif c$$$ c$$$400 continue c$$$500 continue c$$$ c$$$600 continue c$$$ c$$$700 continue c$$$ c$$$ npix=buffer_num c$$$ c$$$ do 900 i=1,512 c$$$ do 800 j=1,512 c$$$ c$$$ if(map_status(i,j) .eq. 1) then c$$$ c$$$ emap(i,j)=map(i,j) c$$$ c$$$ endif c$$$ c$$$800 continue c$$$900 continue c$$$ c$$$ return c$$$ end ********************************************************* subroutine src_width(map,imx,lmx,gbeam,iwid) real map(512,512),gbeam(512,512,0:3) do 200 iw=1,3 do 210 j=-5,5 do 210 i=-5,5 210 if (abs(map(imx+i,lmx+j)) .lt. - gbeam(257+i,257+j,iw)/gbeam(257,257,iw) - *abs(map(imx,lmx))) goto 290 200 continue 290 iwid=iw-1 return end c*********************************************************************** c* cleane2 - subtract extended source Ver.2 * c* map : input,real(512,512) dirty map * c* : output,real(512,512) dirty map, extended source * c* subtracted * c* nantst : input,integer(84) antenna status * c* amx : input,real maximum of map * c* gaine : input,real CLEAN loop gain for extended source * c* emap : input,real(512,512) extended source * c* scale : output,real scaling factor for extended source * c* June 13,1995 koshix * c*********************************************************************** subroutine cleane2(map,mask128,amx,gaine,emap,scale) real gaine,amx,dmx,scale real map(128,128),emap(128,128),mask128(128,128) real a(128,128),b(128,128) do j=1,128 do i=1,128 a(i,j)=emap(i,j) b(i,j)=0. enddo enddo call cft128(a,b,1) do j=1,128 do i=1,128 a(i,j)=a(i,j)*mask128(i,j) b(i,j)=b(i,j)*mask128(i,j) enddo enddo c a(65,65)=0.0 c b(65,65)=0.0 call cft128(a,b,-1) dmx=0. do 400 j=1,128 do 300 i=1,128 if(abs(a(i,j)) .gt. abs(dmx)) dmx=a(i,j) 300 continue 400 continue do j=1,128 do i=1,128 map(i,j)=map(i,j)-a(i,j)/abs(dmx)*abs(amx)*gaine enddo enddo scale=abs(amx)/abs(dmx)*gaine return end c$$$c***************************************************************** c$$$c* contour_trim - find contour trim * c$$$c* * c$$$c* map : input,real(512,512) dirty map * c$$$c* gaine : input,real CLEAN loop gain * c$$$c* for extended source * c$$$c* imx : input,integer * c$$$c* x-coordinate of the maximum point * c$$$c* lmx : input,integer * c$$$c* y-coordinate of the maximum point * c$$$c* npix : output,integer number of pixels(1.-gaine) * c$$$c* 2 : extended source * c$$$c* map_status : output,integer(512,512) * c$$$c* CLEAN point for extended source * c$$$c* e_map : output,real(512,512) extended component * c$$$c* June 13,1995 koshix * c$$$c***************************************************************** c$$$ c$$$ subroutine contour_trim(map,gaine,imx,lmx,npix, c$$$ & map_status,emap) c$$$ c$$$ integer key c$$$ integer i,j,k,l,imx,lmx,ib,lb c$$$ integer npix c$$$ integer n,buffer_num c$$$ integer ibuffer(16384),jbuffer(16384) c$$$ integer map_status(128,128) c$$$ real gaine c$$$ real map(128,128) c$$$ real emap(128,128) c$$$ c$$$ c$$$c********** Contour Trim(1-gaine) of Extended Source ********** c$$$ c$$$ do 200 i=1,128 c$$$ do 100 j=1,128 c$$$ c$$$ map_status(i,j)=0 c$$$ emap(i,j)=0. c$$$ c$$$100 continue c$$$200 continue c$$$ c$$$ do 300 n=1,16384 c$$$ c$$$ ibuffer(n)=-1 c$$$ jbuffer(n)=-1 c$$$ c$$$300 continue c$$$ c$$$ map_status(imx,lmx)=1 c$$$ buffer_num=1 c$$$ ibuffer(buffer_num)=imx c$$$ jbuffer(buffer_num)=lmx c$$$ c$$$ key=1 c$$$ c$$$ do 700 while(key .eq. 1) c$$$ c$$$ key=0 c$$$ c$$$ do 600 n=1,buffer_num c$$$ do 500 k=-1,1 c$$$ do 400 l=-1,1 c$$$ c$$$ ib=mod(ibuffer(n)+k+127,128)+1 c$$$ lb=mod(jbuffer(n)+l+127,128)+1 c$$$ c$$$ if(abs(map(ib,lb)) .ge. c$$$ & abs(map(imx,lmx))*(1.-gaine) .and. c$$$ & map_status(ib,lb) .ne. 1) then c$$$ c$$$ map_status(ib,lb)=1 c$$$ buffer_num=buffer_num+1 c$$$ ibuffer(buffer_num)=ib c$$$ jbuffer(buffer_num)=lb c$$$ c$$$ key=1 c$$$ c$$$ endif c$$$ c$$$400 continue c$$$500 continue c$$$ c$$$600 continue c$$$ c$$$700 continue c$$$ c$$$ npix=buffer_num c$$$ c$$$ do 900 i=1,128 c$$$ do 800 j=1,128 c$$$ c$$$ if(map_status(i,j) .eq. 1) then c$$$ c$$$ emap(i,j)=map(i,j) c$$$ c$$$ endif c$$$ c$$$800 continue c$$$900 continue c$$$ c$$$ return c$$$ end c*********************************************************************** c* mean128_1 - calculate mean and search maximum * c* menu : input,integer clean menu * c* map : input,real(128,128) dirty map * c* am : output,real mean of map * c* amx : output,real maximum of map * c* imx : output,integer x-coordinate of the maximum point * c* lmx : output,integer y-coordinate of the maximum point * c* March 18, 1995 K. Fujiki * c*********************************************************************** subroutine mean128_1(menu,map,am,amx,imx,lmx) c real map(128,128) c am=0. amx=0. if (menu.gt.0) then do 100 l=1,128 do 100 i=1,128 am=am+map(i,l) if (amx.lt.map(i,l)) then amx=map(i,l) imx=i lmx=l end if 100 continue else do 110 l=1,128 do 110 i=1,128 am=am+map(i,l) if (abs(amx).lt.abs(map(i,l))) then amx=map(i,l) imx=i lmx=l end if 110 continue end if am=am/128./128. c return end c c$$$c***************************************************************** c$$$c* contour_trim - find contour trim * c$$$c* * c$$$c* map : input,real(128,128) dirty map * c$$$c* gaine : input,real CLEAN loop gain * c$$$c* for extended source * c$$$c* npix : output,integer number of pixels(1.-gaine) * c$$$c* 2 : extended source * c$$$c* map_status : output,integer(128,128) * c$$$c* CLEAN point for extended source * c$$$c* e_map : output,real(128,128) extended component * c$$$c* March 30 ,1995 Fujiki * c$$$c***************************************************************** c$$$ subroutine contour_trim2(map,gaine,npix,amx,map_status,emap) c$$$ c$$$ integer i,j c$$$ integer npix c$$$ real gaine c$$$ real map(128,128) c$$$ real emap(128,128),map_status(128,128) c$$$ c$$$c********** Contour Trim(1-gaine) of Extended Source ********** c$$$ c$$$ do j=1,128 c$$$ do i=1,128 c$$$ emap(i,j)=0. c$$$ map_status(i,j)=0. c$$$ enddo c$$$ enddo c$$$ c$$$ do 300 j=1,128 c$$$ do 400 i=1,128 c$$$ c$$$c map_status(i,j)=real(int(map(i,j)/amx + gaine)) c$$$ emap(i,j)=map(i,j)/amx-1.+gaine c$$$ map_status(i,j)=real((int(sign(1.,emap(i,j)))+1)/2) c$$$ 400 continue c$$$ 300 continue c$$$ c$$$ do j=1,128 c$$$ do i=1,128 c$$$ emap(i,j)=map(i,j)*map_status(i,j) c$$$ enddo c$$$ enddo c$$$ c$$$ npix=0 c$$$ c$$$ do 900 j=1,128 c$$$ do 800 i=1,128 c$$$ c$$$ npix=npix+int(map_status(i,j)) c$$$ c$$$ 800 continue c$$$ 900 continue c$$$ c$$$ return c$$$ end c***************************************************************** c* src_width128 - estimate the size of gaussian source * c* * c* map : input,real(128,128) dirty map * c* imx,lmx: input,integer source position * c* gbeam : input,real dirty beam * c* iwid : output,integer source size * c* March 30 ,1995 Fujiki * c***************************************************************** subroutine src_width128(map,imx,lmx,gbeam,iwid) real map(128,128),gbeam(128,128,0:40) do 20 iw=1,40 do 20 j=-5,5 do 20 i=-5,5 c imx1=mod(imx+i*4+127,128)+1 lmx1=mod(lmx+j*4+127,128)+1 if (abs(map(imx1,lmx1)) .lt. - gbeam(65+i*4,65+j*4,iw)/gbeam(65,65,iw) - *abs(map(imx,lmx))) then goto 29 endif 20 continue 29 iwid=iw-1 return end