;+ ; NAME: ; slfcal ; PURPOSE: ; imaging of OVSA dat_uv2xy ; ovsa_fft_xy2uv ; slfcal_solve ; slfcal_vj2gn ; slfcal_residu ; ; OUTPUTS: ; map, misc ; COMMENTS: ; On the selfcal cycle: ; 1. get clean visibility using clean cmpts only (i.e. w/o backgr noise) ; 2. improve vis0_ij using the cleaned visibility cvis0_ij as a model ; 3. redo mapping using the renewed visibility (vis1_ij). ; 4. Clean the renewed map. ; ; Note that many things are redone, and repeated. For instance, ; cmap0 should be the same as cmap (input) & the same dbeam is ; calculated three times. The reason is that we want to keep ; using the basic routines developed for clean.-easy to understand ; how things are done in selfcal in comparison with CLEAN. ; However if efficiency is desired, had better replace ; ovsa_fft_uv2xy & ovsa_clean_only with one in which we do not compute dbeam. ; At least we do not call CLEAN itself because it includes determination ; of cbeam and clbox. ; ; Difference between pmap & cmap ; cmap = pmap convolved with cbeam + background noise (just clean) ; pmap = collection of point sources -decomposed (cleaner) ; ; SIDE EFFECTS: ; RESTRICTIONS: ; Works only when there are valid data at all 7 t least 1 time ; MODIFICATION HISTORY: ; Written 05-JAN-2001 by JL ; ; 06-Jun-2008 SDT ; Modified for use of the CLEAN ROI selection in imaging progs ; 22-Sep-2008 SDT ; Modified to use the L1 self-calibration scheme. The program still ; attempts 3 cycles of selfcal, the first to are phase only. In the ; third the user may choose to run a phase only or an amplitude and ; phase selfcal run. If the user chooses to take on an amplitude ; self-cal run the keyword in the call to selfcal_solveL1 must be ; changed from 'phase' to 'phamp' and all the lines relatinf to the ; modulus of the visibilities and its use as a factor must be uncommented ; so that an amplitude drift is corrected for. ; ; 14-Mar-2014 GN ; Replaced call to fft_uv2xy to the renamed ovsa_fft_uv2xy and call to fft_xy2uv to ; the renamed ovsa_fft_xy2uv to avoid potential conflict with the copies of the original ; fft_* routines located in the ssw\hessi tree ; 14-Mar-2014 DG ; Renamed some other routines with prefix ovsa_ for the same reason ; 09-May-2014 DG ; Found that Bmaj and Bmin were for some reason swapped, probably due to error in ; find_axes, which I fixed today and hence uncovered this problem. This is now fixed. ;- pro slfcal, cln, uv_ij, vis_ij, cmap3, misc,noprint=noprint ;----------------- ; set parameters ;----------------- ;m = n_elements(cmap[*,0]) ;not used. Can prob do away with cln if we use clean.pro m=256 ;m=512. niter = 2000 t = cln.t clgain = cln.gain clstop = cln.stop f_GHz = cln.f_GHz ;clbox = cln.box inroi = cln.inroi nroi = cln.nroi ;aa_pxl = cln.Bmin / cln.xyint ;bb_pxl = cln.Bmaj / cln.xyint aa_pxl = cln.Bmaj / cln.xyint bb_pxl = cln.Bmin / cln.xyint PA_deg = cln.PA nsorce = n_elements(clbox) xyint = cln.xyint ; Cellsiz [arcsec] TBmax = cln.TBmax ; Maximum TB [M K] f_GHz = cln.f_GHz ; Centr freq [GHz] tot_uv_siz= 1.0/(xyint/206265.) ;this is the original code uvint = tot_uv_siz/ float(m*2) iu = float(uv_ij)/uvint + m/2 iv = imaginary(uv_ij)/uvint + m/2 ;tfcon=3.14224 vis0_ij = vis_ij & vis_ij=0 ; just for clear notation nt = n_elements(uv_ij[*,0]) nb = n_elements(uv_ij[0,*]) na = 5 cl ={a:aa_pxl, b:bb_pxl, p:PA_deg, inroi:inroi,nroi:nroi, $ f_GHz:f_GHz,niter:niter, gain:clgain, stop:clstop, xyint:xyint} ;------------------------------- ; Exclude time interval where invalid points exists ; This part is done in imagr_vsel & also needs to be improved, so... ;------------------------------- ;------------------------- ; Initial CLEAN map to use to get model visibilities ;------------------------- ovsa_fft_uv2xy, f_GHz, m, vis0_ij, uv_ij, xyint, dbeam, dmap,wsum,/skipbls,noprint=noprint ;,is,js ovsa_clean_only, dmap, dbeam, cl, cmap0, cbeam, pmap0, maxitr0,Cfac,noprint=noprint ;------------------------- ; 1st cycle of selfcal, phase only ;------------------------- TB2flux =wsum/Cfac cvis0_ij = ovsa_fft_xy2uv(pmap0,TB2flux,iu,iv,noprint=noprint) ; see cmnt # 1, SDT vis1_ij = slfcal_solveL1(vis0_ij,cvis0_ij,'phase',noprint=noprint) ; see cmnt # 2 good=where(finite(vis1_ij),complement=bad) if (bad NE [-1]) then vis1_ij[bad]=complex(0.,0.) ovsa_fft_uv2xy, f_GHz, m, vis1_ij, uv_ij, xyint, dbeam, dmap,wsum,/skipbls,noprint=noprint ; cmnt # 3 ovsa_clean_only, dmap, dbeam, cl, cmap1, cbeam, pmap1, maxitr1,clnbm,noprint=noprint ; cmnt # 4 ;-------------------------l ; 2nd cycle of selfcal, phase only ;------------------------- TB2flux =wsum/Cfac cvis1_ij = ovsa_fft_xy2uv(pmap1,TB2flux,iu,iv,noprint=noprint) ; see cmnt # 1 vis2_ij = slfcal_solveL1(vis0_ij,cvis1_ij,'phase',noprint=noprint) ; see cmnt # 2 good=where(finite(vis2_ij),complement=bad) if (bad NE [-1]) then vis2_ij[bad]=complex(0.,0.) ovsa_fft_uv2xy, f_GHz, m, vis2_ij, uv_ij, xyint, dbeam, dmap,wsum,/skipbls,noprint=noprint ; cmnt # 3 ovsa_clean_only, dmap, dbeam, cl, cmap2, cbeam, pmap2, maxitr1,clnbm,noprint=noprint ; cmnt # 4 ;-------------------------l ; 3rd cycle of selfcal, amplitude ;------------------------- TB2flux =wsum/Cfac cvis2_ij = ovsa_fft_xy2uv(pmap2,TB2flux,iu,iv) ; see cmnt # 1, SDT ;modv0=sqrt(float(vis0_ij)^2+imaginary(vis0_ij)^2) ;modulus of all visibilities ;modvc=sqrt(float(cvis2_ij)^2+imaginary(cvis2_ij)^2) ;facc=avg(modv0)/avg(modvc) ;cvis2_ij=cvis2_ij*facc vis3_ij = slfcal_solveL1(vis0_ij,cvis2_ij,'phase',noprint=noprint) ; see cmnt # 2 good=where(finite(vis3_ij),complement=bad) if (bad NE [-1]) then vis3_ij[bad]=complex(0.,0.) ;modv0=sqrt(float(vis0_ij)^2+imaginary(vis0_ij)^2) ;modulus of all visibilities ;modv3=sqrt(float(vis3_ij)^2+imaginary(vis3_ij)^2) ;to correct for amplitude ;fac=avg(modv0)/avg(modv3) ;drift when doing amp slfcal ;vis3_ij=vis3_ij*fac ;Here I'll put an optional section to save out the sefl-caled uv data ;so that it can be imaged independently of the regular programs...temp. ;uvsav=dialog_pickfile(filter='*.sav',default_extension='.sav',/write, $ ; file='SCuvData.sav',/overwrite_prompt) ;save,filename=uvsav,vis3_ij,uv_ij,t,xyint,cl,f_GHz,m ovsa_fft_uv2xy, f_GHz, m, vis3_ij, uv_ij, xyint, dbeam, dmap,wsum,/skipbls,noprint=noprint ; cmnt # 3 ovsa_clean_only, dmap, dbeam, cl, cmap3, cbeam, pmap3, maxitr1,clnbm,noprint=noprint ; cmnt # 4 ;plot_vis2vis,t,uv_ij,vis1_ij,vis0_ij ;ovsa_plot_vis2vis_simpl,vis0_ij,vis1_ij ;---------- Show the results ------- big_tit ='I M A G R / S E L F C A L' ;tit3=['CLEAN MAP','1st SC-MAP','2nd SC-MAP'] ;map3=fltarr(m,m,3) ;map3[*,*,0]=cmap0 ;map3[*,*,1]=cmap1 ;map3[*,*,2]=cmap2 ;xyint = cln.xyint ;;;;winsiz = max(abs(clbox*xyint))*2.5 ; [arcsec] ;winsiz = max(abs(inroi*xyint))*2.5 ; [arcsec] ;SDT, ??? ;ovsa_plot_3maps, map3, tit3, big_tit, xyint, winsiz ;tit3=['CLEAN MAP','1st SC-MAP','2nd SC-MAP'] ;tit2=['CLEAN MAP','1st SC-MAP','1st SC-MAP'] ;map2=fltarr(m,m,3) ;map2[*,*,0]=cmap0 ;map2[*,*,1]=cmap1 ;map2[*,*,2]=cmap1 ;xyint = cln.xyint ;winsiz = max(abs(clbox*xyint))*2.5 ; [arcsec] ;winsiz = max(abs(inroi*xyint))*2.5 ; [arcsec] ;SDT, ??? ;ovsa_plot_3maps, map2, tit2, big_tit, xyint, winsiz ; * * * in lower panels compare differences in visibility ;window,1,xs=800,ys=400 ;!p.multi=[0,3,1] ;pha0=(atan(imaginary(vis0_ij),float(vis0_ij))/!dtor + 360.) mod 360. ;pha1=(atan(imaginary(vis1_ij),float(vis1_ij))/!dtor + 360.) mod 360. ;pha2=(atan(imaginary(vis2_ij),float(vis2_ij))/!dtor + 360.) mod 360. ;amp0=abs(vis0_ij) ;amp1=abs(vis1_ij) ;amp2=abs(vis2_ij) ;y1=(pha1-pha0+360.) mod 360. ;y2=(pha2-pha1+360.) mod 360. ;up=where(y1 gt 180.) & y1(up)=y1(up)-360. ;up=where(y2 gt 180.) & y2(up)=y2(up)-360. ;; before ploting, sort them according to amplitude ;sor = reverse(sort(amp0)) ;plot,y1(sor),psym=3,ytit='Phase Diff. [deg]', $ ; tit = 'Vis 1 - Vis 0' ;plot,y2(sor),psym=3,ytit='Phase Diff. [deg]', $ ; tit = 'Vis 2 - Vis 1' ;plot_io,amp1(sor)/amp0(sor),psym=1,tit='Amplitude Comparison' ; oplot,amp2(sor)/amp1(sor),psym=3 ;window,2,xs=800,ys=400 ;!p.multi=[0,2,1] ;plot,imaginary(vis0_ij(sor)),psym=1 ;oplot,imaginary(vis1_ij(sor)),psym=2 ;oplot,imaginary(vis2_ij(sor)),psym=3 ;plot,float(vis0_ij(sor)),psym=1 ;oplot,float(vis1_ij(sor)),psym=2 ;oplot,float(vis2_ij(sor)),psym=3 ;------------ organize outputs misc=['Remark : CLEAN + Self Calibration', $ 'Remark : SELFCAL, 1st cycle, iterat ='+string(maxitr1), $ 'Remark : SELFCAL, 1st cycle, TB,max ='+string(max(cmap1))] ;'Remark : SELFCAL, 1st cycle, TB,max ='+string(max(cmap1)), $ ;'Remark : SELFCAL, 2nd cycle, iterat ='+string(maxitr2), $ ;'Remark : SELFCAL, 2nd cycle, TB,max ='+string(max(cmap2))] end