diff --git a/c2jtimeline.pro b/c2jtimeline.pro new file mode 100644 index 0000000..a2e7c44 --- /dev/null +++ b/c2jtimeline.pro @@ -0,0 +1,419 @@ +pro c2jtimeline, pickpath=pickpath, range=range, linf=linf + + ; This procedure performs a linear fitting of the cnt2Jy values, considering the + ; lists produced by reducecals.pro, named 'stack_lin_final.txt' and 'stack_cub_final.txt' + ; and grouping data according to a time range defined by the user (expressed in hours). + ; By default, this gap is set to 24 hours. + ; + ; Users must specify if they want to extract linearly fitted cnt2Jy values (linfit value inside each gap), + ; by setting the /linf + ; + ; See comments to identify where the program still needs to be edited. + ; + ; Notice: the program automatically handles the measurements achieved on TP-like files resulting from the + ; conversion of SARDARA acquisitions, which have raw counts levels in the + ; orders of magnitude of 10E+06, 10E+07, producing properly formatted output tables. + ; + ; Authors: Marcello Giroletti, Simona Righini + ; Last edited: Mar 3, 2017 + ; + + sep=path_sep() + + if keyword_set(pickpath) then begin + workpath=dialog_pickfile(/DIRECTORY) + endif else begin + cd, current=curr + workpath=curr+sep + endelse + + + ; if not keyword_set(range) then range=0.08333333 ; 5 minutes + ; if not keyword_set(range) then range=0.16666666 ; 10 minutes + ; if not keyword_set(range) then range=0.25 ; 15 minutes + ; if not keyword_set(range) then range=0.50 ; 30 minutes + ; if not keyword_set(range) then range=1.00 ; 1 hour + if not keyword_set(range) then range=24.0 ; 24 hours + + + if keyword_set(linf) then mode='linf' else mode='aver' + + infiles=strarr(2) + outfiles=strarr(2) + outplot=strarr(2) + infiles[0]='cal_stack_lin_final.txt' + infiles[1]='cal_stack_cub_final.txt' + outfiles[0]='cnt2jy_exlin.txt' + outfiles[1]='cnt2jy_excub.txt' + outplot[0]='cnt2Jy_lin.ps' + outplot[1]='cnt2Jy_cub.ps' + + for f=0, 1 do begin + + checkinput=file_search(workpath+infiles[f],count=fin,/test_read) + if fin eq 0 then begin + print, ' ' + print, '***************************' + print, 'File ',infiles[f],' not found in folder ', workpath + print, '***************************' + break + endif + + ; read and fits the c2Jy from the calibrator scans as a constant or as a function of time + + readcol, workpath+infiles[f], format='(a12,a7,D12.5,d8.1,d8.4,d8.4,d8.4,d8.4,d8.4,d8.4,d8.4,d8.4,d8.4)', name, hms, time, freq, el_d, counts_0, counts_1, err_cnt_0, err_cnt_1, cnt2Jy_0_el, cnt2Jy_1_el, err_cnt2Jy_0_el, err_cnt2Jy_1_el, /SILENT, count=ngood + if ngood eq 0 then begin + readcol, workpath+infiles[f], format='(a12,a7,D12.5,d8.1,d8.4,E12.6,E12.6,E12.6,E12.6,E12.6,E12.6,E12.6,E12.6)', name, hms, time, freq, el_d, counts_0, counts_1, err_cnt_0, err_cnt_1, cnt2Jy_0_el, cnt2Jy_1_el, err_cnt2Jy_0_el, err_cnt2Jy_1_el, /SILENT, count=ngood + endif + + day=strcompress(string(long(time[0])),/remove_all) + t_h=24*(time-floor(time[0])) + timespan=(time[-1]-time[0])*24.0d + elapsed=(time-time[0])*24.0 + + num=n_elements(name) + intnum=ceil(timespan/range) ; number of time intervals with "gap" length + + if timespan eq 0.0 then intnum=1 ; to handle cases where only one calibration measurement is available + + + ; finding the number of unique calibrator names + nameu = name[UNIQ(name, SORT(name))] + numcals=n_elements(nameu) + + ; setup for the output text and graphic file + openw, Unit, workpath+outfiles[f], /GET_LUN + printf, Unit, FORMAT = '(a,a)','Selected interpolation mode = ',mode + printf, Unit, FORMAT = '(a,a,a)','Time range = ',strcompress(string(range,format='(D6.3)'),/remove_all),' hours' + printf, Unit, FORMAT = '(a,a)', 'Number of calibration solution intervals = ',strcompress(string(intnum,format='(I3)'),/remove_all) + + if counts_0[0] le 3000 then begin + ; printf, Unit, FORMAT = '(a2,a9,2(a11),9(a11))',"n","Freq","t_i","t_f","m_0","q_0","e_m_0","e_q_0","m_1","q_1","e_m_1","e_q_1" + printf, Unit, ' n Freq t_i t_f m_0 q_0 e_m_0 e_q_0 m_1 q_1 e_m_1 e_q_1' + endif else begin + printf, Unit, ' n Freq t_i t_f m_0 q_0 e_m_0 e_q_0 m_1 q_1 e_m_1 e_q_1' + endelse + loadct, 13, /silent, ncolors=numcals + + set_plot,'ps' + device, /COLOR, bits=8 + device, file=workpath+outplot[f], /COLOR, /landscape + !p.multi=[1,1,0] + ; overall plot of full dataset + plot, elapsed, cnt2Jy_0_el, yrange=[0,max(cnt2Jy_0_el)+1], xstyle=1, psym=2, ytitle="cnt2Jy", xtitle="Elapsed time [hh.h] since "+day + oplot, elapsed, cnt2Jy_1_el, psym=6 + ; preparing multi-plot page for details + !p.multi=[0,2,2] + + for j=0, intnum-1 do begin + + thisgap=where(elapsed ge elapsed[0]+range*j and elapsed le elapsed[0]+range*(j+1)) ; these are the measurements contained in the time interval at hand + + if thisgap[0] ne -1 then begin + + ti=time[thisgap[0]] + tf=time[thisgap[-1]] + + index0 = counts_0[thisgap] + index1 = counts_1[thisgap] + t_0=time[thisgap] + t_1=time[thisgap] + el_0=el_d[thisgap] + el_1=el_d[thisgap] + c2j_0=cnt2Jy_0_el[thisgap] + c2j_1=cnt2Jy_1_el[thisgap] + err_c2j_0=err_cnt2Jy_0_el[thisgap] + err_c2j_1=err_cnt2Jy_1_el[thisgap] + C_0=dblarr(2) ; for actual interpolation, to be written in output file + C_1=dblarr(2) ; for actual interpolation, to be written in output file + C_sig0=dblarr(2) ; for actual interpolation, to be written in output file + C_sig1=dblarr(2) ; for actual interpolation, to be written in output file + Ch_0=dblarr(2) ; for plot only, in elapsed time + Ch_1=dblarr(2) ; for plot only, in elapsed time + Ch_sig0=dblarr(2) ; for plot only, in elapsed time + Ch_sig1=dblarr(2) ; for plot only, in elapsed time + + + ; selecting only the meaningful values (discarding the "-99.00" cases) + + valid0=WHERE(index0 GT 0, valnum0) + valid1=WHERE(index1 GT 0, valnum1) + + if valnum0 ge 1 then begin + t_0val=dblarr(valnum0) + el_0val=dblarr(valnum0) + c2j_0va=dblarr(valnum0) + err_c2j_0val=dblarr(valnum0) + th_0val=dblarr(valnum0) + endif + + if valnum1 ge 1 then begin + t_1val=dblarr(valnum1) + el_1val=dblarr(valnum1) + c2j_1va=dblarr(valnum1) + err_c2j_1val=dblarr(valnum1) + th_1val=dblarr(valnum1) + endif + + e=findgen(1000) ; for elevation axis representation + e=e/1000.0*90 + + if valid0[0] ne -1 then begin + + ; there are usable measurements + t_0val=t_0[valid0] + el_0val=el_0[valid0] + c2j_0val=c2j_0[valid0] + err_c2j_0val=err_c2j_0[valid0] + th_0val=elapsed[valid0] + + if valnum0 gt 1 then begin + ; there are at least two usable measurements + w_0=1/(err_c2j_0val/c2j_0val)^2 + meanerr,c2j_0val,err_c2j_0val,x_0mean,sigma_0m,sigma_0d + C_0 = LINFIT( t_0val, c2J_0val, CHISQ=chisq, /DOUBLE, MEASURE_ERRORS=err_c2J_0val, PROB=lik_0, SIGMA=C_sig0 ) + Ch_0 = LINFIT(th_0val, c2J_0val, CHISQ=chisq_h, /DOUBLE, MEASURE_ERRORS=err_c2J_0val_h, PROB=lik_0_h, SIGMA=Ch_sig0) + el0plot=el_0val + c2J0plot=c2j_0val + errc2j0plot=err_c2j_0val + th0plot=th_0val + endif else begin + ;only one measurement is available: linear fitting cannot be performed! + x_0mean=c2j_0val[0] + sigma_0d=err_c2j_0val[0] + sigma_0m=err_c2j_0val[0] + C_0[0]=c2j_0val[0] + C_0[1]=0.0 + C_sig0[0]=err_c2j_0val[0] + C_sig0[1]=0.0 + Ch_0[0]=c2j_0val[0] + Ch_0[1]=0.0 + Ch_sig0[0]=err_c2j_0val[0] + Ch_sig0[1]=0.0 + el0plot=[el_0[valid0],el_0[valid0]] + c2J0plot=[c2j_0val,c2j_0val] + th0plot=[th_0val,th_0val] + errc2j0plot=[err_c2j_0val,err_c2j_0val] + endelse + + t=th_0val + y_0=Ch_0[0]+Ch_0[1]*t + y_0min=Ch_0[0]+(Ch_0[1])*th_0val+Ch_sig0[0] + y_0max=Ch_0[0]+(Ch_0[1])*th_0val-Ch_sig0[0] + + ; printouts for test usage only + ; print, '***********************' + ; print, 'Interval ',j+1 + ; print, 'thisgap ', thisgap + ; print, 'ti, tf ', ti, tf + ; print, 'c2J_0 ',c2j_0 + ; print, 'c2J_0val ',c2j_0val + ; print, 'valid0 ', valid0 + ; print, 'th_0val ',th_0val + ; print, '***********************' + ; print, ' ' + ; stop + + ystep=(max(c2J0plot)-min(c2J0plot))/10.0 + labstep=ystep + if ystep eq 0 then begin + ystep=0.1 + labstep=ystep/numcals + endif + plot, el0plot, c2J0plot, /NODATA, xrange=[0,90], yrange=[min(c2J0plot-ystep),max(c2J0plot+ystep)], ytitle="cnt2Jy_0", xtitle="Elevation (deg)" + for i=0,n_elements(nameu)-1 do begin + thiscal=where(name[valid0] eq nameu[i], namenum) + if namenum ge 1 then begin + oplot, el0plot[thiscal], c2J0plot[thiscal], color=i, psym=2 + oploterror, el0plot[thiscal], c2J0plot[thiscal], el_0val[thiscal]*0, errc2j0plot[thiscal], color=i, psym=2 + xyouts, 0, max(c2J0plot)-i*labstep, nameu[i], color=i, charsize=0.8, charthick=2 + endif + endfor + ; oplot, e, y_0, linestyle=0 + ; oplot, e, y_0+sigma_0d, linestyle=1 + ; oplot, e, y_0-sigma_0d, linestyle=1 + + + plot, th0plot, c2J0plot, /NODATA, xs=1, yrange=[min(c2J0plot-ystep),max(c2J0plot+ystep)], ytitle="cnt2Jy_0", xtitle="Elapsed hours since "+day, XTICKFORMAT='(f5.2)' + for i=0,n_elements(nameu)-1 do begin + thiscal=where(name[valid0] eq nameu[i], namenum) + if namenum ge 1 then begin + oplot, th0plot[thiscal], c2J0plot[thiscal], color=i, psym=2 + oploterror, th0plot[thiscal], c2J0plot[thiscal], th0plot[thiscal]*0, errc2j0plot[thiscal], color=i, psym=2 + endif + endfor + + oplot, th_0val, y_0, linestyle=0 + oplot, th_0val, y_0min, linestyle=1 + oplot, th_0val, y_0max, linestyle=1 + + endif else begin + ; there are no usable measurements for CH0! + print, infiles[f]+' Time range '+string(j, format='(I2)')+' - Invalid measurements for CH0' + endelse + + + if valid1[0] ne -1 then begin + ; there are usable measurements + t_1val=t_1[valid1] + el_1val=el_1[valid1] + c2j_1val=c2j_1[valid1] + err_c2j_1val=err_c2j_1[valid1] + th_1val=elapsed[valid1] + + if valnum1 gt 1 then begin + ; there are at least two usable measurements + w_1=1/(err_c2j_1val/c2j_1val)^2 + meanerr,c2j_1val,err_c2j_1val,x_1mean,sigma_1m,sigma_1d + C_1 = LINFIT(t_1val, c2J_1val, CHISQ=chisq, /DOUBLE, MEASURE_ERRORS=err_c2J_1val, PROB=lik_1, SIGMA=C_sig1 ) + Ch_1 = LINFIT(th_1val, c2J_1val, CHISQ=chisq_h, /DOUBLE, MEASURE_ERRORS=err_c2J_1val_h, PROB=lik_1_h, SIGMA=Ch_sig1) + el1plot=el_1val + c2J1plot=c2j_1val + errc2j1plot=err_c2j_1val + th1plot=th_1val + endif else begin + ;only one measurement is available: linear fitting cannot be performed! + x_1mean=c2j_1val[0] + sigma_1d=err_c2j_1val[0] + sigma_1m=err_c2j_1val[0] + C_1[0]=c2j_1val[0] + C_1[1]=0.0 + C_sig1[0]=err_c2j_1val[0] + C_sig1[1]=0.0 + Ch_1[0]=c2j_1val[0] + Ch_1[1]=0.0 + Ch_sig1[0]=err_c2j_1val[0] + Ch_sig1[1]=0.0 + el1plot=[el_1[valid1],el_1[valid1]] + c2J1plot=[c2j_1val,c2j_1val] + th1plot=[th_1val,th_1val] + errc2j1plot=[err_c2j_1val,err_c2j_1val] + endelse + + t=th_1val + y_1=Ch_1[0]+Ch_1[1]*th_1val + y_1min=Ch_1[0]+(Ch_1[1])*th_1val+Ch_sig1[0] + y_1max=Ch_1[0]+(Ch_1[1])*th_1val-Ch_sig1[0] + + ystep=(max(c2J1plot)-min(c2J1plot))/10.0 + labstep=ystep + if ystep eq 0 then begin + ystep=0.1 + labstep=ystep/numcals + endif + plot, el1plot, c2J1plot, /NODATA, xrange=[0,90], yrange=[min(c2J1plot-ystep),max(c2J1plot+ystep)], ytitle="cnt2Jy_1", xtitle="Elevation (deg)" + for i=0,n_elements(nameu)-1 do begin + thiscal=where(name[valid1] eq nameu[i], namenum) + if namenum ge 1 then begin + oplot, el1plot[thiscal], c2J1plot[thiscal], color=i, psym=6 + oploterror, el1plot[thiscal], c2J1plot[thiscal], el1plot[thiscal]*0, errc2J1plot[thiscal], color=i, psym=6 + xyouts, 0, max(c2J1plot)-i*labstep, nameu[i], color=i, charsize=0.8, charthick=2 + endif + endfor + + + plot, th1plot, c2J1plot, /NODATA, xs=1, yrange=[min(c2J1plot-ystep),max(c2J1plot+ystep)], ytitle="cnt2Jy_1", xtitle="Elapsed hours since "+day, XTICKFORMAT='(f5.2)' + for i=0,n_elements(nameu)-1 do begin + thiscal=where(name[valid1] eq nameu[i], namenum) + if namenum ge 1 then begin + oplot, th1plot[thiscal], c2J1plot[thiscal], color=i, psym=6 + oploterror, th1plot[thiscal], c2J1plot[thiscal], th1plot[thiscal]*0, errc2J1plot[thiscal], color=i, psym=6 + endif + endfor + oplot, th_1val, y_1, linestyle=0 + oplot, th_1val, y_1min, linestyle=1 + oplot, th_1val, y_1max, linestyle=1 + + endif else begin + ; there are no usable measurements for CH1! + print, infiles[f]+' Time range '+string(j, format='(I2)')+' - Invalid measurements for CH1' + endelse + + + case mode of + 'aver': begin + printf, Unit, FORMAT = '(i4," ",d8.1,2(d11.4),8(g13.7))',j+1,freq[0],ti,tf,0,x_0mean,0,sigma_0m,0,x_1mean,0,sigma_1m + end + 'linf': begin + printf, Unit, FORMAT = '(i4," ",d8.1,2(d11.4),8(g13.7))',j+1,freq[0],ti,tf,C_0[1],C_0[0],C_sig0[1],C_sig0[0],C_1[1],C_1[0],C_sig1[1],C_sig1[0] + end + endcase + + endif else begin + print, 'No measurements in time interval number ', strcompress(string(j+1),/remove_all) + printf, Unit, 'Interval number ',strcompress(string(j+1),/remove_all),' does not contain measurements' + print, 'Going tp the next step' + endelse + + endfor + + endfor + + device, /close + close, /ALL + + print, ' ******************************' + print, ' **** C2JTIMELINE is DONE ****' + print, ' *** Next step is RUNTARGET ***' + print, ' ******************************' + + return +end + + + +;pro fit,ascissa,somma,ch, + + +pro meanerr,x,sigmax,xmean,sigmam,sigmad + ;+ + ; NAME: + ;meanerr + ; PURPOSE: (one line) + ;Calculate the mean and estimated errors for a set of data points + ; DESCRIPTION: + ;This routine is adapted from Program 5-1, XFIT, from "Data Reduction + ;and Error Analysis for the Physical Sciences", p. 76, by Philip R. + ;Bevington, McGraw Hill. This routine computes the weighted mean using + ;Instrumental weights (w=1/sigma^2). + ; CATEGORY: + ;Statistics + ; CALLING SEQUENCE: + ;meanerr,x,sigmax,xmean,sigmam,sigmad + ; INPUTS: + ;x - Array of data points + ;sigmax - array of standard deviations for data points + ; OPTIONAL INPUT PARAMETERS: + ;None. + ; KEYWORD PARAMETERS: + ;None. + ; OUTPUTS: + ;xmean - weighted mean + ;sigmam - standard deviation of mean + ;sigmad - standard deviation of data + ; COMMON BLOCKS: + ;None. + ; SIDE EFFECTS: + ;None. + ; RESTRICTIONS: + ;None. + ; PROCEDURE: + ; MODIFICATION HISTORY: + ;Written by Marc W. Buie, Lowell Observatory, 1992 Feb 20 + ;- + + if n_elements(x) eq 1 then begin + xmean = x[0] + sigmam = sigmax[0] + sigmad = sigmax[0] + endif else begin + weight = 1.0/sigmax^2 + sum = total(weight) + if sum eq 0.0 then print,'MEANERR: sum is zero.' + sumx = total(weight*x) + xmean = sumx/sum + sigmam = sqrt(1.0/sum) + sigmad = sqrt(total((x-xmean)^2)/(n_elements(x)-1)) + endelse + +end \ No newline at end of file diff --git a/calflux.pro b/calflux.pro new file mode 100644 index 0000000..97efb1b --- /dev/null +++ b/calflux.pro @@ -0,0 +1,337 @@ +pro calflux, namecal, nu, intSnu, intSnu_err + + ; Last Edited: Oct 6, 2017 + + ; Adapted from SDI_CALFLUX_RELISE_V1: CREATION DATE 4/5/2016; VERSION 1 BY Elise Egron, Alberto Pellizzoni. + ; Addition of Ott et al flux density estimate: Simona Righini - such part of the code is now muted, though + ; P&B *interpolated* flux densities are provided. + + ; EXAMPLE: calflux, '3C295', 7.2 + ; EXAMPLE (with outputs): calflux, '3C295', 7.2, Snu, intSnu, Snu_err, intSnu_err, pbchoice + + ; INPUTS: + ; namecal = calibrator name (string) + ; nu = required frequency (GHz) + + ; OUTPUTS: + ; + ; intSnu: flux density (Jy) at nu (error in intSnu_err ) + ; ("INTERPOLATION" based on polynomial coeff. from Perley et al. 2013) + ; + ; AVAILABLE CALIBRATORS; 3C123, 3C286, 3C295, 3C48, 3C147, NGC7027 + + + ; ****************************** + + ; CALIBRATORS' DATA from Perley et al. 2013 + + ; frequency table in GHz (33 entries) + + f_tab=[0.3275,1.015,1.275,1.465,1.865,2.565,3.565,4.535,4.835,4.885,6.135,6.885,7.465,8.435,8.485,8.735,11.06,12.890,14.635,14.715,14.915,14.965,17.422,18.230,18.485,18.585,20.485,22.460,22.835,24.450,25.836,26.485,28.450] + + + ; CALIBRATOR DATA: fluxes (Jy), flux errors (Jy), polynomial coeff. (log Jy), polynomial coeff. errors (log Jy), RA/DEC coords. J2000 (deg) + + ; ****** 3C123 ****** + + S_3C123=[145.0,66.2,46.6,47.8,38.7,28.9,21.4,16.9,16.0,15.88,12.81,11.20,11.01,9.20,9.10,8.86,6.73,6.035,5.34,5.02,5.132,5.092,4.272,4.181,4.090,3.934,3.586,3.297,3.334,2.867,2.697,2.716,2.436] + + S_3C123_err=[4.3,4.3,3.2,0.5,0.6,0.3,0.8,0.2,0.2,0.1,0.15,0.14,0.2,0.04,0.15,0.05,0.15,0.15,0.05,0.05,0.025,0.028,0.07,0.07,0.055,0.055,0.055,0.022,0.06,0.03,0.06,0.05,0.06] + + A_3C123=[1.8077d,-0.8018d,-0.1157d,0.d] + A_3C123_err=[0.0036d,0.0081d,0.0047d,0.d] + + coord_3C123=[69.268230d,29.670505d] + + ; NOTES: + ; missing flux value at 12.890 GHz: estimated (6.73+5.34)/2.=6.035 + ; missing flux value at 18.230 GHz: estimated (4.272+4.090)/2.=4.181 + ; missing error value at 12.890 GHz: estimated 0.15 + ; missing error value at 18.230 GHz: estimated 0.07 + ; coords. 04:37:04.3752, +29:40:13.818 + + + ; ****** 3C286 ****** + + S_3C286=[26.1,18.4,13.8,15.0,13.2,10.9,9.5,7.68,7.33,7.297,6.49,5.75,5.70,5.059,5.045,4.930,4.053,3.662,3.509,3.375,3.399,3.387,2.980,2.860,2.925,2.880,2.731,2.505,2.562,2.387,2.181,2.247,2.079] + + S_3C286_err=[0.8,4.3,2.0,0.2,0.2,0.2,0.1,0.1,0.2,0.046,0.15,0.05,0.10,0.021,0.07,0.024,0.08,0.070,0.040,0.040,0.016,0.015,0.04,0.045,0.045,0.04,0.05,0.0160,0.05,0.03,0.06,0.05,0.05] + + A_3C286=[1.2515d,-0.4605d,-0.1715d,0.0336d] + A_3C286_err=[0.0048d,0.0163d,0.0208d,0.0082d] + + coord_3C286=[202.784533d,30.509155d] + + ; NOTES: + ; coords. 13:31:08.2879, +30:30:32.958 + + + ; ****** 3C295 ****** + + S_3C295=[60.8,30.8,21.5,22.2,17.9,12.8,9.62,6.96,6.45,6.37,4.99,4.21,4.13,3.319,3.295,3.173,2.204,1.904,1.694,1.630,1.626,1.617,1.311,1.222,1.256,1.221,1.089,0.952,0.967,0.861,0.770,0.779,0.689] + + S_3C295_err=[1.8,7.3,3.0,0.5,0.3,0.2,0.2,0.09,0.15,0.04,0.05,0.05,0.07,0.014,0.05,0.016,0.05,0.04,0.04,0.03,0.008,0.007,0.025,0.05,0.020,0.015,0.015,0.005,0.015,0.020,0.02,0.0200,0.020] + + A_3C295=[1.4866d,-0.7871d,-0.3440d,0.0749d] + A_3C295_err=[0.0036d,0.0110d,0.0160d,0.0070d] + + coord_3C295=[212.835279d, 52.202644d] + + ; NOTES: + ; coords. 14:11:20.467, +52:12:09.52 + + + ; ****** 3C48 ****** + + A_3C48=[1.3324d, -0.7690d, -0.1950d, 0.059d] + A_3C48_err=[0.005d,0.01d,0.02d,0.02d] ; NOT REPORTED IN PERLEY (Derived by ee, ap) + + coord_3C48=[24.422081d, 33.159759d] + + ; NOTES: + ; polynomial coeff. for 2012 (table 11 p.13 Perley et al. 2013); ERRORS NOT AVAILABLE ! + ; + ; coords. 01:37:41.2994, +33:09:35.132 + + + ; ****** 3C147 ****** + + A_3C147=[1.4616d, -0.7187d, -0.2424d, 0.079d] + A_3C147_err=[0.01d, 0.02d, 0.02d, 0.03d] ; NOT REPORTED IN PERLEY (Derived by ee, ap) + + coord_3C147=[85.650575d, 49.852009d] + + ; NOTES: + ; polynomial coeff. for 2012 (table 11 p.13 Perley et al. 2013); ERRORS NOT AVAILABLE ! + ; + ; coords. 05:42:36.1379, +49:51:07.233 + + + ; ****** NGC7027 ****** + + f_tab_NGC7027=[1.465, 4.885, 8.435, 14.965, 22.460, 43.340] ; required frequency table + + S_NGC7027=[1.561d, 5.392d, 5.851d, 5.702d, 5.480d, 5.091d] + S_NGC7027_err=[2., 5., 5., 7., 13., 42.]/1000d + + S_NGC7027_delta=[3.6, -3.0, -7.1, -7.4, -7.0, -5.0]/1000d + + S_NGC7027=S_NGC7027+S_NGC7027_delta*12. ; VALID FOR 2012 + + f00=5.5215-3.8*15./1000. + f11=5.9294-7.12*15./1000. + + coord_NGC7027=[316.756638d, 42.236163d] + + ; NOTES: + ; polynomial coeff. NOT AVAILABLE + ; Assuming interpolated values for approx. 2012 (table 12 p.17 Perley et al. 2013) ! + ; + ; coords. 21:07:01.5930, +42:14:10.186 + + + ;*************************** + + + S_name=['3C123','3C286','3C295', '3C48', '3C147', 'NGC7027'] + + nn=n_elements(S_name) + + nfreq_tab=n_elements(f_tab) + S_tot=fltarr(nn,2,nfreq_tab) + S_tot(0,0,*)=S_3C123 + S_tot(1,0,*)=S_3C286 + S_tot(2,0,*)=S_3C295 + S_tot(0,1,*)=S_3C123_err + S_tot(1,1,*)=S_3C286_err + S_tot(2,1,*)=S_3C295_err + + + if (strtrim(strupcase(namecal),2)) eq 'NGC7027' then begin + f_tab=f_tab_NGC7027 + nfreq_tab=n_elements(f_tab_NGC7027) + S_tot=fltarr(nn,2,nfreq_tab) + S_tot(5,0,*)=S_NGC7027 + S_tot(5,1,*)=S_NGC7027_err + endif + + + A_tot=dblarr(nn,2,4) + A_tot(0,0,*)=A_3C123 + A_tot(1,0,*)=A_3C286 + A_tot(2,0,*)=A_3C295 + A_tot(3,0,*)=A_3C48 + A_tot(4,0,*)=A_3C147 + A_tot(0,1,*)=A_3C123_err + A_tot(1,1,*)=A_3C286_err + A_tot(2,1,*)=A_3C295_err + A_tot(3,1,*)=A_3C48_err + A_tot(4,1,*)=A_3C147_err + + Coord_tot=dblarr(nn,2) + Coord_tot(0,*)=coord_3C123 + Coord_tot(1,*)=coord_3C286 + Coord_tot(2,*)=coord_3C295 + Coord_tot(3,*)=coord_3C48 + Coord_tot(4,*)=coord_3C147 + Coord_tot(5,*)=coord_NGC7027 + + wname=where(strtrim(strupcase(namecal),2) eq S_name,ncal) + + if ncal eq 0 then begin + print + print,'CALIBRATOR NOT FOUND!' + print + stop + endif + + S=S_tot(wname,0,*) + S_err=S_tot(wname,1,*) + + A=A_tot(wname,0,*) + A_err=A_tot(wname,1,*) + + ra=Coord_tot(wname,0) + dec=Coord_tot(wname,1) + + wmin=where(f_tab le nu) + numin=max(f_tab(wmin)) + wmin=where(f_tab eq numin) + + wmax=where(f_tab ge nu) + numax=min(f_tab(wmax)) + wmax=where(f_tab eq numax) + + Snu=(S(wmax)*(nu-numin)+S(wmin)*(numax-nu))/(numax-numin) + Snu_err=(S_err(wmax)*(nu-numin)+S_err(wmin)*(numax-nu))/(numax-numin) + Snu=Snu(0) + Snu_err=Snu_err(0) + + intSnu_=A(0)+A(1)*alog10(nu)+A(2)*(alog10(nu))^2.+A(3)*(alog10(nu))^3. + intSnu=(10.d0)^intSnu_ + intSnu_err=A_err(0) ; TBD more precise error propagation ****** !!!!!! ******* + intSnu=intSnu(0) + intSnu_err=intSnu_err(0) + + if not keyword_set(verb) then verb=1 ; for standalone usage of procedure, otherwise it fails + + if verb eq 1 then begin + +; print +; print,'**************************************************' +; print,'CALIBRATOR: ', strtrim(namecal,2) +; print,'**************************************************' +; print,'>>> Fluxes by Perley et al. 2013' + +; if intSnu_ ne 0 then begin +;; print, 'INTERPOLATED FLUX (Jy): ', strtrim(string(intSnu),2), ' +/- ', strtrim(string(intSnu_err),2) +; endif +; +; if Snu ne 0 then begin +;; print, 'EXTRAPOLATED FLUX (Jy): ',strtrim(string(Snu),2),' +/- ',strtrim(string(Snu_err),2) +; endif +; +; if Snu ne 0 and intSnu_ ne 0 then begin +;; print, 'DISCREPANCY (Jy): ', strtrim(string(abs(Snu-intSnu)),2) +;; print, 'DISCREPANCY %: ', strtrim(string(abs(Snu-intSnu)/Snu*100.),2) +; pbchoice='INT' +; endif +; +; if Snu eq 0 or Snu eq '-NaN' then begin +;; print, 'EXTRAPOLATED FLUX NOT AVAILABLE!' +; Snu=intSnu +; Snu_err=intSnu_err +; pbchoice='INT' +; endif +; + if A(0) eq 0 and A(1) eq 0 and A(2) eq 0 and A(3) eq 0 then begin + ; print, 'INTERPOLATED FLUX NOT AVAILABLE! Assigning extrapolated one' + intSnu=Snu + intSnu_err=Snu_err + pbchoice='EXT' + endif + + +; OTT ET AL. PART + +; print,'--------------------------------------------------' +; print,'>>> Fluxes by Ott et al. 1994' + + case strupcase(namecal) of + '3C48': begin + a=2.465 + b=-0.004 + c=-0.125 + low=1.408 + high=23.780 + end + + '3C123': begin + a=2.525 + b=0.246 + c=-0.1638 + low=1.408 + high=23.780 + end + + '3C147': begin + a=2.806 + b=-0.140 + c=-0.1031 + low=0.500 + high=23.780 + end + + '3C286': begin + a=0.956 + b=0.584 + c=-0.1644 + low=1.408 + high=43.200 + end + + '3C295': begin + a=1.470 + b=0.765 + c=-0.2545 + low=1.408 + high=32.000 + end + + 'NGC7027': begin + a=1.322 + b=-0.134 + c=0.0 + low=10.550 + high=43.200 + end + + endcase + + ott_flux=10^(a+b*alog10(nu*1000.0)+c*alog10(nu*1000.0)*alog10(nu*1000.0)) + + inrange='Y' + ; warning on improper use of polynomial, i.e. when observed frequency is outside the validity range + if (nu lt low) or (nu gt high) then begin +; print, 'WARNING: requested frequency lies +; print, ' outside Ott et al. range + inrange='N' + endif + + +; print, 'OTT INTERPOLATED FLUX (Jy): ', strtrim(string(ott_flux,format='(D7.4)'),2) +; print, 'DISCREPANCY WITH PERLEY (Jy): ', strtrim(string(abs(ott_flux-intSnu),format='(D7.4)'),2) +; print, 'DISCREPANCY WITH PERLEY %: ', strtrim(string(abs(ott_flux-intSnu)/ott_flux*100.,format='(D5.2)'),2) +; print,'**************************************************' +; print + + endif + + ;stop + + + +end + + + diff --git a/calibfit.pro b/calibfit.pro new file mode 100644 index 0000000..2339c9f --- /dev/null +++ b/calibfit.pro @@ -0,0 +1,179 @@ +pro calibresiduals, X,Y,datafit,Unit2,peak,err, x_peak, rasd, decsd, tipo, SNR, peak_cnt, err_cnt, Level, n_off, off, residual + + ;Minor section, to compute the rms-noise and the signal-to-noise ratio + + fit = LINFIT(X[0:20], Y[0:20], YFIT=datalinfit) + subtracted = Y[0:20]-datalinfit + mom = moment(subtracted) + noise = sqrt(mom[1]) + residual = Y-datafit + mom = moment(residual) + resnoise = sqrt(mom[1]) + PRINTF, Unit2, "noise: ", noise, " counts; noise residuals: ", resnoise, " counts" + ; ;;print, "noise: ", noise, " counts; noise residuals: ", resnoise, " counts" + SNR=peak/noise + PRINTF, Unit2, "SNR: ", SNR + peak_cnt = peak + err_cnt = err + Level = (Y[-1] + Y[0])/2. + + case tipo of + 0: begin + meanlat=x_peak/180d*!dpi ; radians + off=(x_peak-decsd) ; degrees + end + 1: begin + meanlat=decsd ; radians + off=(x_peak-rasd)*cos(meanlat) ; degrees + end + endcase + n_off=1 + + return +end + + + +pro getcnt2jy, flux, peak_cnt, err_cnt, tau0, datael, cnt2Jy, err_cnt2Jy, Unit2, Level + + ; Minor section, for the computation of conversion factor (counts --> Jy), + ; inclusive of opacity correction. + + peak_cnt=peak_cnt*exp(mean(tau0)/sin(datael)) + err_cnt=err_cnt*exp(mean(tau0)/sin(datael)) + cnt2Jy = flux/peak_cnt + err_cnt2Jy = cnt2Jy*err_cnt/peak_cnt + + if flux gt 1000 then begin + cnt2Jy = -99.000 + err_cnt2Jy = -99.000 + endif + + PRINTF, Unit2, " " + PRINTF, Unit2, "cnt2Jy = (", cnt2Jy, "+-", err_cnt2Jy, " [Jy/cnt]" + PRINTF, Unit2, " " + PRINTF, Unit2, "Level = ", Level, " [cnt]" + PRINTF, Unit2, " " + return +end + + + +pro calibfit, scanflag,stacflag,polyflag,section,tipo,allpath,namefile,Out3,flux,datael,tau0, xx,yy,ii,ff,x_mid,Nsamples,sd_sub,gpos,decsd,rasd, fwhm, n_off, off, p, e, c, d, SNR, plo, doplot + + +; Main procedure, devoted to the fitting operations, both in the linear and cubic cases, for targets. +; +; Authors: Marcello Giroletti, Simona Righini +; Last edited: Oct 5, 2017 +; + + ; scanflag = scan buono (=1) o cattivo (=0) + ; stacflag = single o stacked scan + ; polyflag = gaussiana più grado del polinomio (primo='linear', terzo='cubic') + ; section = Ch_0 o Ch_1 + ; tipo = flag di direzione - convertito qua sotto in direflag + if (tipo eq 1) then direflag = 'RA' else direflag = 'DEC' + + ; common logistics + sep=path_sep() + if scanflag eq 0 then begin + + n_off = 0 + off = 0. + p = -99. + e = -99. + c = -99. + d = -99. + SNR = -99. + if doplot eq 'y' then dataplot, 0, [0,1], [0,1], namefile, allpath, section, 0, 0, 0, 0, 0, direflag, polyflag, plo, stacflag + + endif else begin + + decs=decsd*!dpi/180.0 + media = (yy[ii]+yy[ff])/2. + ; define new coordinates in which the source is centred at x=0 and the baseline is about 0 at the source position + X=xx[ii:ff]-gpos + Y=yy[ii:ff]-media + + if (polyflag eq 'linear') then begin + + datafit = GAUSSFIT(X, Y, A_fit, ESTIMATES=[(Y[x_mid]),0.0,sd_sub,(X[0]*Y[-1]-X[-1]*Y[0])/(X[0]-X[-1]),(Y[-1]-Y[0])/(X[-1]-X[0])], NTERMS=5, SIGMA=A_sig) +; print, ' Ampl(cnt) Peak_pos(deg) sigma(deg) q(cnt) m(cnt/deg)' +; print, 'Estimates =',(Y[x_mid]),0.0,sd_sub,(X[0]*Y[-1]-X[-1]*Y[0])/(X[0]-X[-1]),(Y[-1]-Y[0])/(X[-1]-X[0]) +; print, 'Fit =',A_fit +; print, 'Sigma =',A_sig + if (tipo eq 0) then begin + fwhm=2*SQRT(2*ALOG(2))*A_fit[2]*60.0 ; in arcmin, if fit is done on coordinates in degrees ; *(speed/60.)*60. to be used if done in samples + endif else begin + fwhm=2*SQRT(2*ALOG(2))*A_fit[2]*60.0*cos(decs) ; in arcmin, if fit is done on coordinates in degrees ; *(speed/60.)*60. to be used if done in samples + endelse + x_peak=A_fit[1]+gpos ; coordinate corresponding to gaussian peak + PRINTF, Out3, " GAUSSIAN + LINEAR " + PRINTF, Out3, section+" Gaussian fit -- amplitude:", A_fit[0], " +- ", A_sig[0], " counts" + PRINTF, Out3, section+" Gaussian fit -- sigma: ", A_fit[2], " +- ", A_sig[2], " [s]" + PRINTF, Out3, " corresponding to FWHM: ", fwhm, " [arcmin]" + + calibresiduals, X,Y,datafit,Out3,A_fit[0],A_sig[0], x_peak, rasd, decsd, tipo, SNR, peak_cnt, err_cnt, Level, n_off, off, residual + + ; computing and accounting for the RMS of the residuals in the overall amplitude error estimate + nres=n_elements(residual) + resrange=ceil(nres/5.0) + rescut=[residual[0:resrange],residual[-1*resrange,-1]] ; avoiding the central part of the subscan, where artifacts can be present due to sidelobes + resstat=moment(rescut) + res_rms=sqrt(resstat[1]) + err_cnt=sqrt(err_cnt^2+res_rms^2+(0.03*peak_cnt)^2) ; updated error for the amplitude measurement, including a default 3% uncertainty on calibrator flux-amplitude + + getcnt2Jy, flux, peak_cnt, err_cnt, tau0, datael, cnt2Jy, err_cnt2Jy, Out3, Level + if doplot eq 'y' then dataplot, 1, X, Y, namefile, allpath, section, datafit, cnt2Jy, Level, SNR, residual, direflag, polyflag, plo, stacflag + + endif else begin + + ; estimates for fit with A0*exp(-(X-A1)^2/(2*A2))+A3+A4*X+A5*X^2+A6*X^3; ie cubic+gaussian + A = [Y[x_mid],0.0,sd_sub,(X[0]*Y[-1]-X[-1]*Y[0])/(X[0]-X[-1]),(Y[-1]-Y[0])/(X[-1]-X[0]),0.0,0.0] + A_GUESS = A + ; compute the parameters without weights (DUM_W not set) + datafit = CURVEFIT(X,Y,DUM_W, A, SIGMA, FITA=[1,1,1,1,1,1,1], FUNCTION_NAME='cal_gfunct', /DOUBLE, STATUS=suc_fit) +; print, ' A0=Ampl(cnt) A1=Peak_pos(deg) A2=sigma(deg) A3(cnt) A4(cnt/deg) A5(cnt/deg2) A6(cnt/deg3)' +; print, 'Function estimates: ', A_GUESS +; print, 'Function parameters: ', A +; print, 'Sigma parameters: ', SIGMA + if (tipo eq 0) then begin + fwhm=2*SQRT(2*ALOG(2))*A[2]*60.0 ; in arcmin, if fit is done on coordinates in degrees ; *(speed/60.)*60. to be used if done in samples + endif else begin + fwhm=2*SQRT(2*ALOG(2))*A[2]*60.0*cos(decs) ; in arcmin, if fit is done on coordinates in degrees ; *(speed/60.)*60. to be used if done in samples + endelse + x_peak=A[1]+gpos ; coordinate corresponding to gaussian peak + PRINTF, Out3, " GAUSSIAN + CUBIC " + PRINTF, Out3, section+" Gaussian fit -- amplitude:", A[0], " +- ", SIGMA[2], " counts" + PRINTF, Out3, section+" Gaussian fit -- sigma: ", A[0], " +- ", SIGMA[0], " [s]" + PRINTF, Out3, " corresponding to FWHM: ", fwhm, " [arcmin]" + + calibresiduals, X,Y,datafit,Out3,A[0],A[2], x_peak, rasd, decsd, tipo, SNR, peak_cnt, err_cnt, Level, n_off, off, residual + + ; computing and accounting for the RMS of the residuals in the overall amplitude error estimate + nres=n_elements(residual) + resrange=ceil(nres/5.0) + rescut_b=[residual[0:resrange],residual[-1*resrange,-1]] ; avoiding the central part of the subscan, where artifacts can be present due to sidelobes + resstat=moment(rescut_b) + res_rms=sqrt(resstat[1]) + err_cnt=sqrt(err_cnt^2+res_rms^2+(0.03*peak_cnt)^2) ; updated error for the amplitude measurement + + getcnt2Jy, flux, peak_cnt, err_cnt, tau0, datael, cnt2Jy, err_cnt2Jy, Out3, Level + if doplot eq 'y' then dataplot, 1, X, Y, namefile, allpath, section, datafit, cnt2Jy, Level, SNR, residual, direflag, polyflag, plo, stacflag + + endelse + + el_d=datael*180.0/!dpi + p=peak_cnt + e=err_cnt + c=cnt2Jy + d=err_cnt2Jy + + endelse + + return + +end + + diff --git a/dataplot.pro b/dataplot.pro new file mode 100644 index 0000000..f297c32 --- /dev/null +++ b/dataplot.pro @@ -0,0 +1,54 @@ +pro dataplot, scanflag, X, Y, allnamefile, mypath, section, datafit, cnt2Jy, Level, SNR, residual, direflag, tipofit, pdis, stac + + ; This procedure is in charge of the plotting operations + ; + ; Authors: Marcello Giroletti, Simona Righini + ; Last edited: Oct 5, 2017 + ; + + sep=path_sep() + if strmatch(allnamefile, '*'+sep+'*') then begin + splitted=strsplit(allnamefile,sep,/extract) + namefile=splitted[-1] + endif else begin + namefile=allnamefile + endelse + plot_title = namefile+" "+section+" "+direflag+' '+tipofit + + if (scanflag eq 0) then begin + if (section eq 'Ch_0') then begin + pdis=plot(X, Y, buffer=1, layout=[2,2,1], dimensions=[1024,768], title = plot_title, /NODATA) + pdis=plot(X, Y, buffer=1, layout=[2,2,2], dimensions=[1024,768], title = "RESIDUALS", /NODATA, /CURRENT) + endif else begin + pdis=plot(X, Y, buffer=1, layout=[2,2,3], dimensions=[1024,768], title = plot_title, /NODATA, /CURRENT) + pdis=plot(X, Y, buffer=1, layout=[2,2,4], dimensions=[1024,768], title = "RESIDUALS", /NODATA, /CURRENT) + endelse + endif else begin + if (section eq 'Ch_0') then begin + nw=1 + ; pdis=plot(X, Y, layout=[2,2,nw], dimensions=[1024,768],$ + pdis=plot(X, Y, buffer=1, layout=[2,2,nw], dimensions=[1024,768],$ + title = plot_title, xtitle = "x_J2000 [deg]", xrange = [X[0],X[-1]], ytitle = "signal [cnt]") + endif else begin + nw=3 + ; pdis=plot(X, Y, layout=[2,2,nw], /current, $ + pdis=plot(X, Y, buffer=1, layout=[2,2,nw], /current, $ + title = plot_title, xtitle = "x_J2000 [deg]", xrange = [X[0],X[-1]], ytitle = "signal [cnt]") + endelse + pdis.symbol=2 + pdis.sym_size=0.2 + + ; overplot of the gaussian fit + pdis=plot(X, datafit, /overplot) + ; plot residuals in side panel + ; pdis=plot(X, residual, layout=[2,2,nw+1], /current, title="RESIDUALS") + pdis=plot(X, residual,buffer=1, layout=[2,2,nw+1], /current, title="RESIDUALS") + pdis.xtitle = "x_J2000 [deg]" + pdis.ytitle = "signal [cnt]" + endelse + if stac eq 'single' then begin + if tipofit eq 'linear' then pdis.save, mypath+sep+namefile+'_lin_single.pdf', /LANDSCAPE, RESOLUTION=300 + if tipofit eq 'cubic' then pdis.save, mypath+sep+namefile+'_cub_single.pdf', /LANDSCAPE, RESOLUTION=300 + endif + return +end diff --git a/findh0key.pro b/findh0key.pro new file mode 100644 index 0000000..5eada70 --- /dev/null +++ b/findh0key.pro @@ -0,0 +1,33 @@ +pro findh0key, header, stringa, key, value, info, rflag + + ; Procedure to read keyword values from FITS headers + ; + ; Author: Simona Righini + ; Last edited: Oct 6, 2017 + + + rflag=1 + search=strmatch(header, stringa) + index=where(search eq 1) + sp=path_sep() + if index ne -1 and stringa eq '*Converted FITS*' then return + if index ne -1 and stringa eq '*Obtained from ESCS 0.1 TP FITS*' then return + if index ne -1 then begin + content=strsplit(header[index],'=',/extract) + key=content[0] + value_and_info=strsplit(string(content[1]),'/',/extract) + if stringa ne '*ScheduleName =*' then begin + value=value_and_info[0] + info=value_and_info[1] + endif else begin + steps=n_elements(value_and_info) + info=value_and_info[-1] + value=value_and_info[0] + for s=1, steps-2 do begin + value=value+sp+value_and_info[s] + endfor + endelse + endif else begin + rflag=0 + endelse +end \ No newline at end of file diff --git a/newflag.pro b/newflag.pro new file mode 100644 index 0000000..7d76763 --- /dev/null +++ b/newflag.pro @@ -0,0 +1,490 @@ +; Programma per l'ispezione visuale ed il flagging dei file FITS (ESCS 0.2-0.3) +; Creato il 6 febbraio 2013 da A.Mattana e S.Righini +; Debug: +; 7 feb - corretta mancanza variabile "app" in start, caso di creazione checkfile +; 12 feb - corretto problema con Logical Units, aggiungendo free_lun +; 21 feb - tentativo di introdurre tasti "BACK" e "SKIP", modificando il metodo +; con cui viene costruita la lista dei file ancora da flaggare +; 28 feb - corretto errore nel caso di "skip" su ultimo subscan +; Modifiche: +; 27/04/16 - la folder da selezionare all'inizio è ora quella "parent" che contiene +; tutte le subfolders coi dati da flaggare. Nell'interfaccia PLOT è ora +; possibile saltare alla folder precedente o successiva. Aggiunte varie +; labels e migliorata l'impaginazione. + +; Come usare il tool: +; 1) Lanciare con il nome flagging_widget_develop; +; 2) selezionare l'utente dal menu a tendina in alto; +; 3) cliccare sul pulsante BROWSE per selezionare la folder che contiene i FITS; +; 4) cliccare sul pulsante GO!; +; 5) se esiste già un checkfile per quella folder, comparirà una finestra +; nella quale scegliere "SI" se si vuole continuare il flagging (sui file non ancora +; flaggati) oppure "NO" per ripartire da capo, nel qual caso il precedente +; checkfile verrà perso! +; Se, comunque, si è risposto "SI" ma il checkfile pre-esistente era completo +; (cioè tutti i FITS erano stati flaggati) il programma chiede un'altra folder; +; 6) nella nuova finestra di PLOT compariranno uno alla volta i FITS, per ciascuno +; selezionare quali dati si vuole TENERE, usando i pulsanti; +; 7) col pulsante "Back" si torna indietro di un subscan, il cui precedente flagging +; viene cancellato. Il pulsante "Skip" salta il subscan in esame. Non è possibile +; invocare "Back" dopo uno "Skip" (il bottone, per sicurezza, viene disattivato); +; 8) per interrompere, cliccare sul pulsante "STOP and EXIT". +; + +pro WID_BASE_0, GROUP_LEADER=wGroup, _EXTRA=_VWBExtra_ + + common wid, WID_BASE_0, WID_BOX, wid_left_id, wid_right_id $ + ,WID_DRAW_LEFT, Button_BROWSE $ + ,TEXT_OBS_PATH, Button_START, WID_TAB $ + ,Button_LEFT,Button_RIGHT,Button_PERFECT,Button_ALL,Button_NONE,WID_BUTTON_STOP $ + ,WID_BUTTON_SKIP, WID_BUTTON_BACK, WID_BUTTON_PREVDIR, WID_BUTTON_NEXTDIR + + common cfg, folder, dirlist, dirindex, userlist, user, outname, sp, list, $ + flagged, tobeflagged, flagnext, flags, app, index + + sp=path_sep() + + Resolve_Routine, 'WID_BASE_0_event';,/COMPILE_FULL_FILE ; Load event callback routines + + WID_BASE_0 = Widget_Base( GROUP_LEADER=wGroup, UNAME='WID_BASE_0' $ + ,XOFFSET=5 ,YOFFSET=5 ,SCR_XSIZE=1200 ,SCR_YSIZE=750 $ + ,TITLE="FLAGGING: WHAT'S TO KEEP?" ,SPACE=3 ,XPAD=3 ,YPAD=3) + + WID_BOX = Widget_Base( GROUP_LEADER=boxes, UNAME='WID_BOX' $ + ,XOFFSET=5 ,YOFFSET=5 ,SCR_XSIZE=400 ,SCR_YSIZE=300 $ + ,TITLE="ENTER CHECKFILE CUSTOM NAME" ,SPACE=3 ,XPAD=3 ,YPAD=3) + + WID_TAB = Widget_Tab(WID_BASE_0, UNAME='WID_TAB' $ + ,XOFFSET=5 ,YOFFSET=5 ,SCR_XSIZE=1200 ,SCR_YSIZE=750 $ + ) + + WID_TAB_CFG = Widget_Base(WID_TAB, UNAME='WID_TAB_CFG' $ + ,XOFFSET=5 ,YOFFSET=5 $ + ,TITLE=" CONFIG ") + + WID_TAB_PLOT = Widget_Base(WID_TAB, UNAME='WID_TAB_PLOT' $ + ,XOFFSET=5 ,YOFFSET=5 $ + ,TITLE=" PLOT ") + + WID_DRAW_LEFT = Widget_Draw(WID_TAB_PLOT, UNAME='WID_DRAW_LEFT' $ + ,XOFFSET=83 ,YOFFSET=25 ,SCR_XSIZE=450 ,SCR_YSIZE=350) + + WID_DRAW_RIGHT = Widget_Draw(WID_TAB_PLOT, UNAME='WID_DRAW_RIGHT' $ + ,XOFFSET=566 ,YOFFSET=25 ,SCR_XSIZE=450 ,SCR_YSIZE=350) + + + ;Materiale nel TAB CFG + userlist=['NA','AT','EV','MG','SR','AZ','MM'] + flags = ['-1c0','-1c1','+0cx','+1cx','-1cx'] + user='NA' + + WID_USER_LABEL = Widget_Label(WID_TAB_CFG, UNAME='WID_USER_LABEL' $ + ,XOFFSET=300 ,YOFFSET=50 ,SCR_XSIZE=250 ,SCR_YSIZE=50 $ + ,/ALIGN_RIGHT ,VALUE='Select User ID (or leave default value)') + + DROPLIST_USERS = Widget_Droplist(WID_TAB_CFG, $ + UNAME='DROPLIST_USERS' ,XOFFSET=300 ,YOFFSET=100 $ + ,SCR_XSIZE=250 ,SCR_YSIZE=40, value=userlist) + + WID_FOLDER_LABEL = Widget_Label(WID_TAB_CFG, UNAME='WID_FOLDER_LABEL' $ + ,XOFFSET=300 ,YOFFSET=200 ,SCR_XSIZE=400 ,SCR_YSIZE=40 $ + ,/ALIGN_LEFT ,VALUE='Select main folder (containing the subfolders with FITS files)') + + Button_BROWSE = Widget_Button(WID_TAB_CFG, UNAME='Button_BROWSE' $ + ,XOFFSET=200 ,YOFFSET=240 ,SCR_XSIZE=100 ,SCR_YSIZE=50 $ + ,/ALIGN_CENTER ,VALUE='BROWSE') + + TEXT_OBS_PATH = Widget_Text(WID_TAB_CFG, UNAME='TEXT_OBS_PATH' $ + ,XOFFSET=300 ,YOFFSET=240 ,SCR_XSIZE=600 ,SCR_YSIZE=50 $ + ,XSIZE=40,YSIZE=1) + + Button_START = Widget_Button(WID_TAB_CFG, UNAME='Button_START' $ + ,XOFFSET=450 ,YOFFSET=400, SCR_XSIZE=150 ,SCR_YSIZE=150 $ + ,/ALIGN_CENTER ,VALUE='START!') + + WID_HEADER_LABEL = Widget_Label(WID_TAB_PLOT, UNAME='WID_HEADER_LABEL' $ + ,XOFFSET=375 ,YOFFSET=390 ,SCR_XSIZE=350 ,SCR_YSIZE=40 $ + ,/ALIGN_CENTER ,VALUE='Indicate the section(s) to be kept for analysis') + + Button_LEFT = Widget_Button(WID_TAB_PLOT, UNAME='Button_LEFT' $ + ,XOFFSET=200 ,YOFFSET=450 ,SCR_XSIZE=100 ,SCR_YSIZE=100 $ + ,/ALIGN_CENTER ,VALUE='LEFT') + + Button_RIGHT = Widget_Button(WID_TAB_PLOT, UNAME='Button_RIGHT' $ + ,XOFFSET=350 ,YOFFSET=450 ,SCR_XSIZE=100 ,SCR_YSIZE=100 $ + ,/ALIGN_CENTER ,VALUE='RIGHT') + + Button_PERFECT = Widget_Button(WID_TAB_PLOT, UNAME='Button_PERFECT' $ + ,XOFFSET=500 ,YOFFSET=450 ,SCR_XSIZE=100 ,SCR_YSIZE=100 $ + ,/ALIGN_CENTER ,VALUE='BOTH (V.GOOD)') + + Button_ALL = Widget_Button(WID_TAB_PLOT, UNAME='Button_ALL' $ + ,XOFFSET=650 ,YOFFSET=450 ,SCR_XSIZE=100 ,SCR_YSIZE=100 $ + ,/ALIGN_CENTER ,VALUE='BOTH') + + Button_NONE = Widget_Button(WID_TAB_PLOT, UNAME='Button_NONE' $ + ,XOFFSET=800 ,YOFFSET=450 ,SCR_XSIZE=100 ,SCR_YSIZE=100 $ + ,/ALIGN_CENTER ,VALUE='NONE') + + WID_BUTTON_PREVDIR = Widget_Button(WID_TAB_PLOT, $ + UNAME='WID_BUTTON_PREVDIR' ,XOFFSET=180 ,YOFFSET=600 $ + ,SCR_XSIZE=120 ,SCR_YSIZE=50 ,/ALIGN_CENTER ,VALUE='<< PREV DIR') + + WID_BUTTON_BACK = Widget_Button(WID_TAB_PLOT, $ + UNAME='WID_BUTTON_BACK' ,XOFFSET=350 ,YOFFSET=600 $ + ,SCR_XSIZE=100 ,SCR_YSIZE=50 ,/ALIGN_CENTER ,VALUE='< BACK') + + WID_BUTTON_STOP = Widget_Button(WID_TAB_PLOT, $ + UNAME='WID_BUTTON_STOP' ,XOFFSET=500 ,YOFFSET=600 $ + ,SCR_XSIZE=100 ,SCR_YSIZE=50 ,/ALIGN_CENTER ,VALUE='STOP and EXIT') + + WID_BUTTON_SKIP = Widget_Button(WID_TAB_PLOT, $ + UNAME='WID_BUTTON_SKIP' ,XOFFSET=650 ,YOFFSET=600 $ + ,SCR_XSIZE=100 ,SCR_YSIZE=50 ,/ALIGN_CENTER ,VALUE='SKIP >') + + WID_BUTTON_NEXTDIR = Widget_Button(WID_TAB_PLOT, $ + UNAME='WID_BUTTON_NEXTDIR' ,XOFFSET=800 ,YOFFSET=600 $ + ,SCR_XSIZE=120 ,SCR_YSIZE=50 ,/ALIGN_CENTER ,VALUE='NEXT DIR >>') + + + Widget_Control, /REALIZE, WID_BASE_0 + WIDGET_CONTROL, WID_DRAW_RIGHT, GET_VALUE=wid_right_id + WIDGET_CONTROL, WID_DRAW_LEFT, GET_VALUE=wid_left_id + widget_control, WID_BUTTON_PREVDIR, SENSITIVE=0 + widget_control, WID_BUTTON_BACK, SENSITIVE=0 + widget_control, WID_BUTTON_STOP, SENSITIVE=0 + widget_control, WID_BUTTON_SKIP, SENSITIVE=0 + widget_control, WID_BUTTON_NEXTDIR, SENSITIVE=0 + widget_control, Button_LEFT, SENSITIVE=0 + widget_control, Button_RIGHT, SENSITIVE=0 + widget_control, Button_PERFECT, SENSITIVE=0 + widget_control, Button_ALL, SENSITIVE=0 + widget_control, Button_NONE, SENSITIVE=0 + + + + XManager, 'WID_BASE_0', WID_BASE_0, /NO_BLOCK + +end + + +pro WID_BASE_0_event, Event + + common wid + + common cfg + + Device, decomposed=0 + + ;help,event,/str + wTarget = (widget_info(Event.id,/NAME) eq 'TREE' ? $ + widget_info(Event.id, /tree_root) : event.id) + + wWidget = Event.top + + case wTarget of + + ;Selezione utente + Widget_Info(wWidget, FIND_BY_UNAME='DROPLIST_USERS'): begin + if( Tag_Names(Event, /STRUCTURE_NAME) eq 'WIDGET_DROPLIST' )then begin $ + user=userlist(event.index) + endif + end + + + ;Selezione FITS folder + Widget_Info(wWidget, FIND_BY_UNAME='Button_BROWSE'): begin + if( Tag_Names(Event, /STRUCTURE_NAME) eq 'WIDGET_BUTTON' )then begin $ + folder=dialog_pickfile(/DIRECTORY) + Widget_Control, TEXT_OBS_PATH, set_value=folder + dirlist=file_search(folder,'2*',count=dirnum,/TEST_DIRECTORY,/FULLY_QUALIFY_PATH) + outname=strarr(dirnum) + for thisfol=0, dirnum-1 do begin + dirs=strsplit(dirlist[thisfol],sp,/extract) + outname[thisfol]='checkfile_'+strcompress(dirs[-1],/remove_all)+'.txt' + endfor + endif +end + +;Inizio del flagging +Widget_Info(wWidget, FIND_BY_UNAME='Button_START'): begin + if( Tag_Names(Event, /STRUCTURE_NAME) eq 'WIDGET_BUTTON' )then begin $ + if keyword_set(folder) then begin + list=file_search(dirlist[0]+sp+'2*.fits', count=numb) + flagged=intarr(numb) + flagged[*]=0 + if numb gt 0 then begin + dirindex=0 + start + endif else begin + choice=dialog_message("No FITS in first folder!",/INFO,/CENTER) + endelse + endif else begin + choice=dialog_message("No folder was selected!",/INFO,/CENTER) + endelse +endif +end + +;Accetta solo il canale left +Widget_Info(wWidget, FIND_BY_UNAME='Button_LEFT'): begin + if( Tag_Names(Event, /STRUCTURE_NAME) eq 'WIDGET_BUTTON' )then begin + widget_control, WID_BUTTON_BACK, SENSITIVE=1 + flagga, list[tobeflagged[index]], 0, app, tobeflagged[index] + endif +end + +;Accetta solo il canale right +Widget_Info(wWidget, FIND_BY_UNAME='Button_RIGHT'): begin + if( Tag_Names(Event, /STRUCTURE_NAME) eq 'WIDGET_BUTTON' )then begin + widget_control, WID_BUTTON_BACK, SENSITIVE=1 + flagga, list[tobeflagged[index]], 1, app, tobeflagged[index] + endif +end + +;Accetta entrambi e valutali "perfetti" +Widget_Info(wWidget, FIND_BY_UNAME='Button_PERFECT'): begin + if( Tag_Names(Event, /STRUCTURE_NAME) eq 'WIDGET_BUTTON' )then begin + widget_control, WID_BUTTON_BACK, SENSITIVE=1 + flagga, list[tobeflagged[index]], 2, app, tobeflagged[index] + endif +end + +;Accetta entrambi nonostante imperfezioni +Widget_Info(wWidget, FIND_BY_UNAME='Button_ALL'): begin + if( Tag_Names(Event, /STRUCTURE_NAME) eq 'WIDGET_BUTTON' )then begin + widget_control, WID_BUTTON_BACK, SENSITIVE=1 + flagga, list[tobeflagged[index]], 3, app, tobeflagged[index] + endif +end + +;Rifiuta entrambi i canali +Widget_Info(wWidget, FIND_BY_UNAME='Button_NONE'): begin + if( Tag_Names(Event, /STRUCTURE_NAME) eq 'WIDGET_BUTTON' )then begin + widget_control, WID_BUTTON_BACK, SENSITIVE=1 + flagga, list[tobeflagged[index]], 4, app, tobeflagged[index] + endif +end + +;Back to previous subscan +Widget_Info(wWidget, FIND_BY_UNAME='WID_BUTTON_BACK'): begin + if( Tag_Names(Event, /STRUCTURE_NAME) eq 'WIDGET_BUTTON' )then begin + if tobeflagged[index] eq 0 then begin + choice=dialog_message("Not applicable, this is the first subscan!",/INFORMATION,/CENTER) + return + endif + readcol, dirlist[dirindex]+sp+outname[dirindex], fname, flag, usid, format='(A,A,A)' + nitems=n_elements(fname) + openw, Unit, dirlist[dirindex]+sp+outname[dirindex], /get_lun + for i=0, nitems-2 do begin + printf, Unit, fname[i], flag[i], usid[i], format='(A,1X,A,1X,A)' + endfor + close, Unit + free_lun, Unit + index=index-1 + plot_ch + endif +end + +;Skip this subscan, go to the next one +Widget_Info(wWidget, FIND_BY_UNAME='WID_BUTTON_SKIP'): begin + if( Tag_Names(Event, /STRUCTURE_NAME) eq 'WIDGET_BUTTON' )then begin $ + if index eq n_elements(tobeflagged)-1 then begin + choice=dialog_message("Not applicable: this is the last subscan!",/INFORMATION,/CENTER) + return + endif + index=index+1 + widget_control, WID_BUTTON_BACK, SENSITIVE=0 + plot_ch +endif +end + +;Stop and exit +Widget_Info(wWidget, FIND_BY_UNAME='WID_BUTTON_STOP'): begin + if( Tag_Names(Event, /STRUCTURE_NAME) eq 'WIDGET_BUTTON' )then $ + Widget_Control, WID_BASE_0, /DESTROY +end + +;Back to previous folder +Widget_Info(wWidget, FIND_BY_UNAME='WID_BUTTON_PREVDIR'): begin + if( Tag_Names(Event, /STRUCTURE_NAME) eq 'WIDGET_BUTTON' )then begin + if dirindex eq 0 then begin + choice=dialog_message("Not applicable, this was the first folder!",/INFORMATION,/CENTER) + return + endif + dirindex=dirindex-1 + list=file_search(dirlist[dirindex]+sp+'2*.fits', count=numb) + flagged=intarr(numb) + flagged[*]=0 + if numb gt 0 then begin + start + endif else begin + choice=dialog_message("No FITS in first folder!",/INFO,/CENTER) + endelse + endif else begin + choice=dialog_message("No folder was selected!",/INFO,/CENTER) + endelse +end + + +;On to next folder +Widget_Info(wWidget, FIND_BY_UNAME='WID_BUTTON_NEXTDIR'): begin + if( Tag_Names(Event, /STRUCTURE_NAME) eq 'WIDGET_BUTTON' )then begin + if dirindex eq n_elements(dirlist)-1 then begin + choice=dialog_message("Not applicable, this was the last folder!",/INFORMATION,/CENTER) + return + endif + dirindex=dirindex+1 + list=file_search(dirlist[dirindex]+sp+'2*.fits', count=numb) + flagged=intarr(numb) + flagged[*]=0 + if numb gt 0 then begin + start + endif else begin + choice=dialog_message("No FITS in first folder!",/INFO,/CENTER) + endelse + endif else begin + choice=dialog_message("No folder was selected!",/INFO,/CENTER) + endelse +end + +else: +endcase + +end + + + + +; +; Empty stub procedure used for autoloading. +; +pro newflag, GROUP_LEADER=wGroup, _EXTRA=_VWBExtra_ + WID_BASE_0, GROUP_LEADER=wGroup, _EXTRA=_VWBExtra_ +end + +;############################################################## +pro Show_Message, comment + + common wid + + widget_control, WID_LABEL_COMMENT, Set_Value=comment + +end +;############################################################## + +pro plot_ch + common cfg + common wid + ;how many files to be flagged? + numb=n_elements(list) + ; load FITS file + data=MRDFITS(list[tobeflagged[index]],4,/silent) + ; extract the filename from full path + segments=strsplit(list[tobeflagged[index]],sp,/extract) + segcount=n_elements(segments) + filename=segments[segcount-1] + ; count the data samples and generate x-axis + nsample=n_elements(data.time) + xaxis=indgen(nsample) + ; ch0 plot + wset, wid_left_id + plot, xaxis[10:nsample-10], data[10:nsample-10].ch0, $ + ys=1, background=255, color=0, title=filename + ; ch1 plot + wset, wid_right_id + plot, xaxis[10:nsample-10], data[10:nsample-10].ch1, $ + ys=1, background=255, color=0, title='Dir'+strcompress(string(dirindex+1))+' out of'+strcompress(n_elements(dirlist))+' - FITS'+strcompress(string(tobeflagged[index]+1))+' out of'+strcompress(string(numb)) +end + + +pro start + common wid + common cfg + + widget_control, Button_LEFT, SENSITIVE=1 + widget_control, Button_RIGHT, SENSITIVE=1 + widget_control, Button_PERFECT, SENSITIVE=1 + widget_control, Button_ALL, SENSITIVE=1 + widget_control, Button_NONE, SENSITIVE=1 + widget_control, WID_BUTTON_PREVDIR, SENSITIVE=1 + widget_control, WID_BUTTON_BACK, SENSITIVE=1 + widget_control, WID_BUTTON_STOP, SENSITIVE=1 + widget_control, WID_BUTTON_SKIP, SENSITIVE=1 + widget_control, WID_BUTTON_NEXTDIR, SENSITIVE=1 + + if file_test(dirlist[dirindex]+sp+outname[dirindex]) then begin + choice=dialog_message("Checkfile already exists: resume and append? 'NO' means you'll start again and overwrite!!!",/QUESTION,/CENTER) + if choice eq 'No' then begin + app='No' + endif else begin + readcol, dirlist[dirindex]+sp+outname[dirindex], fname, flag, usid, format='(A,A,A)' + checklines=n_elements(fname) + match_index=0 + ; evidenzio i file già flaggati + for i=0, n_elements(list)-1 do begin + for j=0, checklines-1 do begin + if strmatch(list[i],'*'+fname[j]) then flagged[i]=1 + endfor + endfor + app='Yes' + if (total(flagged) ge n_elements(list)) then begin + choice=dialog_message("All FITS were already flagged. Change folder.",/INFORMATION,/CENTER) + return + endif + endelse + tobeflagged=where(flagged eq 0, howmany) + endif else begin + tobeflagged=where(flagged eq 0, howmany) + openw, Unit, dirlist[dirindex]+sp+outname[dirindex], /GET_LUN + close, Unit + free_lun, Unit + app='Yes' + endelse + Widget_Control, WID_TAB, SET_TAB_CURRENT=1 + index=0 + flagnext=tobeflagged[index] + plot_ch +end + +pro flagga, current, flag_id, appflag, ind + common cfg + common wid + slash=strpos(current,sp,/reverse_search) + current=strmid(current,slash+1,strlen(current)-slash-1) + if (appflag eq "No") and (ind eq 0) then begin + openw, Unit1, dirlist[dirindex]+sp+outname[dirindex], /get_lun + endif else begin + openw, Unit1, dirlist[dirindex]+sp+outname[dirindex], /get_lun, /append + endelse + printf, Unit1, current, flags[flag_id], user, format='(A,1X,A,1X,A)' + ; print, 'Writing in ', Unit1 + close, Unit1 + free_lun, Unit1 + if index lt n_elements(tobeflagged)-1 then begin + index=index+1 + flagnext=tobeflagged[index] + endif else begin + finale=dialog_message('No more files to flag!',/INFORMATION,/CENTER) + widget_control, Button_LEFT, SENSITIVE=0 + widget_control, Button_RIGHT, SENSITIVE=0 + widget_control, Button_PERFECT, SENSITIVE=0 + widget_control, Button_ALL, SENSITIVE=0 + widget_control, Button_NONE, SENSITIVE=0 + widget_control, WID_BUTTON_SKIP, SENSITIVE=0 + widget_control, WID_BUTTON_BACK, SENSITIVE=0 + endelse + if flagnext lt n_elements(list) then begin plot_ch +endif else begin + finale=dialog_message('No more files to flag!',/INFORMATION,/CENTER) + widget_control, Button_LEFT, SENSITIVE=0 + widget_control, Button_RIGHT, SENSITIVE=0 + widget_control, Button_PERFECT, SENSITIVE=0 + widget_control, Button_ALL, SENSITIVE=0 + widget_control, Button_NONE, SENSITIVE=0 + widget_control, WID_BUTTON_SKIP, SENSITIVE=0 + widget_control, WID_BUTTON_BACK, SENSITIVE=0 +endelse +end + + + diff --git a/runcalib.pro b/runcalib.pro new file mode 100644 index 0000000..22ce014 --- /dev/null +++ b/runcalib.pro @@ -0,0 +1,1247 @@ +pro runcalib, pickpath=pickpath, skypath=skypath, sub=sub, linear=linear, cubic=cubic, plot=plot, extlist=extlist, skipgc=skipgc + + ; This procedure is devoted to the reduction of cross scans acquired on flux calibrators and + ; achieve the conversion factors cnt->Jy. + ; It plots and analyses the data obtained by stacking all the subscans contained in any scan + ; (i.e. subfolder) stored in the working folder. + ; If users want to graphically choose the data folder, they must use: + ; + ; IDL> runcalib, /pickpath in order to select the path to the data parent folder (otherwise the chosen path is ./CALIBRATORS) + ; IDL> runcalib, /skypath in order to select the path to the folder containing skydips (otherwise the chosen path is ./SKYDIPS) + ; + ; Measurements for all the single subscans can be obtained by explicitly choosing this option: + ; + ; IDL> runcalib, /sub + ; + ; The gaussian fitting performed on the (sub)scans can be run using three different baseline-fitting options, + ; in particular: + ; + ; IDL> runcalib, /linear performs, and outputs, only the linear fitting of baselines + ; IDL> runcalib, /cubic performs, and outputs, only the cubic fitting of baselines + ; IDL> runcalib, /both performs, and outputs, both the cubic and linear fitting of baselines + ; + ; By default, the choice is 'both'. + ; + ; At present, graphic output is exceedingly slow, in particular for single subscans (if the /sub option is selected). + ; So it is disabled by default. + ; To enable it, explicitly set the option: + ; + ; IDL> runcalib, /plot + ; + ; If non-standard calibrators are to be considered, users must specify the /extlist option, + ; which allows them to pick a proper external list, where source names and nominal flux + ; densities (Jy) are provided in the simple format (A,D). + ; + ; The option + ; + ; IDL> runcalib, /skipgc + ; + ; overrides the use of the standard gain curves; it poses Gain=1 for all the elevations, + ; thus practically avoiding the elevation-dependent gain compensation. + ; + ; COMPATIBILITY WITH CONVERTED FILES + ; FITS files achieved with older versions of DISCOS (dating back to years 2008-2015) + ; need to be converted into the proper FITS format (ask for procedure named updatefits.pro). + ; The program automatically handles the TP-like files resulting from the + ; conversion of SARDARA acquisitions, which have raw counts levels in the + ; orders of magnitude of 10E+06..10E+07, producing properly-formatted output tables. + ; + ; Authors: Marcello Giroletti, Simona Righini + ; Last edited: Oct 18, 2017 + ; + + + common logistics, workpath, calpath, sep, skydippath, myextlist, fitchoice, dosingle, doplot + + close, /all ; to close possible pending text files + + sep=path_sep() + + if keyword_set (sub) then dosingle='y' else dosingle='n' + + if keyword_set(pickpath) then begin + workpath=dialog_pickfile(/DIRECTORY, TITLE='Please pick PARENT FOLDER containing CALIBRATION SCANS') + endif else begin + CD, CURRENT=curr + workpath=curr+sep+'CALIBRATORS'+sep + endelse + + if keyword_set(skypath) then begin + skydippath=dialog_pickfile(/DIRECTORY, TITLE='Please pick folder containing file Tau.txt') + endif else begin + CD, CURRENT=curr + workpath=curr+sep+'SKYDIPS'+sep + endelse + + if keyword_set(linear) then begin + fitchoice='linear' + endif else begin + if keyword_set(cubic) then fitchoice='cubic' else fitchoice='both' + endelse + + if keyword_set(plot) then doplot='y' else doplot='n' + + if keyword_set(extlist) then begin + myextlist=dialog_pickfile(TITLE='Please pick EXTERNAL LIST of calibrators') + endif else begin + myextlist=' ' + endelse + + cal_stack, freq=freq + return + +end + + +PRO cal_gfunct, X, A, F, pder + EX2 = A[0]*EXP(-(X-A[1])^2/(2*A[2]^2)) + F = EX2 + A[3] + A[4]*X +A[5]*X^2 +A[6]*X^3 + IF N_PARAMS() GE 4 THEN $ + pder = [[EX2/A[0]], [EX2*(X-A[1])/(A[2]^2)], [EX2*((X-A[1])^2)/(A[2]^3)], [replicate(1.0, N_ELEMENTS(X))], [X], [X^2], [X^3]] +END + + +function cal_wmean, val, dval, nan = nan, error = error + compile_opt idl2 + on_error, 2 + + ;- check inputs + if n_params() ne 2 then begin + return, !values.f_nan + endif + sz = n_elements(val) + if n_elements(dval) ne sz then $ + message, 'val and dval must contain the same number of elements' + + + ;- the calculation + ;- requires at least one finite value + if keyword_set(nan) then begin + hit = where(finite(val) and finite(dval), ct) + if ct eq 0 then begin + error = !values.f_nan + return, !values.f_nan + endif + endif + + dw = total(1D / dval^2, nan = keyword_set(nan)) + sum = total(1D * val / dval^2, nan = keyword_set(nan)) / dw + error = sqrt(1D / dw) + + return, sum +end + + +pro cal_stack, path=path, out=out, plot=plot, beam=beam, speed=speed, dt=dt, source=source, flux=flux, freq=freq + common logistics + + ; !EXCEPT=2 ; enabling the indications of location/line numbers in exception messages + !EXCEPT=0 ; silencing the printout of exception messages + + ; reads first file from first scan to set basic info in output files + list = FILE_SEARCH(workpath+'20*', COUNT=number, /TEST_DIRECTORY) + + timings=dblarr(number+1) + gap=dblarr(number+1) + timings[0]=systime(1) + gap[0]=0 + + sublist = FILE_SEARCH(list[0]+sep+'20*.fits', COUNT=subnumber) + ; XXX ToDO: ripetere la lettura per un file di ogni sub-folder e verificare che siano omogenei + data=MRDFITS(sublist[0],4,/SILENT) + RFinfo=MRDFITS(sublist[0],2,/SILENT) + Sectinfo=MRDFITS(sublist[0],1,/SILENT) + + ; read frequency and bandwidth from RF table + bw=RFinfo[0].bandWidth + + freq=RFinfo[0].frequency+bw/2.0 + dt=1.0/(Sectinfo[0].sampleRate*1e6) + + ; which telescope was used? Let's get to a 3-char code to be used afterwards + firstmainh=mrdfits(sublist[0],0,firsthead0,/silent) + findh0key, firsthead0, '*ANTENNA =*', keylab, sitename, infolab, siteflag + splitsitename=strsplit(sitename,"'",/extract) + cleansite=strupcase(strcompress(splitsitename[1],/remove_all)) + site=strmid(cleansite,0,3) ; this way, possible cases will be MED, SRT and NOT + + case site of + 'MED': begin + beam=38.7/(freq/1000.0) ; beam in arcmin + ; stop + ; frequency dependent normalized gain polynomials + ; http://www.med.ira.inaf.it/ManualeMedicina/English/5.%20Efficiency.htm + ; taken on 2015 September 17 + if (freq gt 4000 and freq lt 6000) then begin + ; this is C band + A_0=[7.9212981E-1,6.2403816E-3,-4.6834953E-5,0] + A_1=[7.9212981E-1,6.2403816E-3,-4.6834953E-5,0] + endif + if (freq gt 8000 and freq lt 10000) then begin + ; this is X band + A_0=[0.730487,0.00773332,-5.54743e-05,0] ; updated 2015, but it is not opacity-corrected + A_1=[0.730487,0.00773332,-5.54743e-05,0] ; updated 2015, but it is not opacity-corrected + ; A_0=[6.1059261E-01,1.0623634E-02,-7.2457279E-05,0] ; origin??? + ; A_1=[6.1059261E-01,1.0623634E-02,-7.2457279E-05,0] ; origin??? + endif + if (freq gt 18000 and freq le 18500) then begin + ; this is K-low band + A_0=[0.8373589,0.005140111,-4.058370E-5,0] ; as of Sept.30, 2015 + A_1=[0.8373589,0.005140111,-4.058370E-5,0] ; as of Sept.30, 2015 + endif + if (freq gt 18500 and freq le 26500) then begin + ; this is K-high band + A_0=[0.7929185,0.005900533,-4.203179E-5,0] ; as of Sept.30, 2015 + A_1=[0.7929185,0.005900533,-4.203179E-5,0] ; as of Sept.30, 2015 + endif + end + 'SRT': begin + ; frequency dependent normalized gain polynomials + ; as provided by Andrea Orlati (commissioning) + if (freq gt 4000 and freq lt 8000) then begin + ; this is C band + beam=18.937/(freq/1000.0) ; beam in arcmin + A_0=[0.939813,0.00139,-0.000009,0] ; XXX To be updated + A_1=[0.939813,0.00139,-0.000009,0] ; XXX To be updated + endif + if (freq ge 18000 and freq lt 27000) then begin + ; this is K band + beam=18.264/(freq/1000.0) ; beam in arcmin + if freq lt 21000 then begin + A_0=[0.508145,0.014839,-0.000111,0] ; from official SRT page + A_1=[0.508145,0.014839,-0.000111,0] ; from official SRT page + endif + if freq ge 21000 and freq le 23000 then begin + A_0=[0.606733,0.01132,-0.000081,0] ; from official SRT page + A_1=[0.606733,0.01132,-0.000081,0] ; from official SRT page + endif + if freq gt 23000 then begin + A_0=[0.533025,0.015561,-0.000129,0] ; from official SRT page + A_1=[0.533025,0.015561,-0.000129,0] ; from official SRT page + endif + endif + end + 'NOT': begin + ;beam=38.7/(freq/1000.0) ; beam in arcmin - CLONED FROM MEDICINA XXX + if (freq gt 4000 and freq lt 6000) then begin + ; this is C band + beam=7.5 ; fixed value! + ; A_0=[0.98463659,0.00017053038,1.9348601e-09,0] ; old + ; A_1=[0.98463659,0.00017053038,1.9348601e-09,0] ; old + A_0=[0.95227062926,0.00198726460815,-2.0685473288e-05,0] ; Cassaro, Oct 16th 2017 + A_1=[0.95227062926,0.00198726460815,-2.0685473288e-05,0] ; Cassaro, Oct 16th 2017 + endif + if (freq gt 8000 and freq lt 10000) then begin + ; this is X band + print, 'X Curve for Noto is not available, yet. Returning.' + return + endif + if (freq gt 18000 and freq le 26500) then begin + ; this is K band + print, 'K Curve for Noto is not available, yet. Returning.' + return + endif + end + else: begin + print, 'UNKNOWN SITE: cannot apply gain curves. Returning.' + return + end + endcase + beamd=beam/60. + sd=beamd/(2*SQRT(2*ALOG(2))) + + + ; open the main output files + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + openw, Unit0, workpath+'cal_stack_lin_final.txt', /GET_LUN + openw, Unit1, workpath+'cal_stack_lin_sep.txt', /GET_LUN + printf, Unit0, FORMAT = '(A7,1X,A3)', "Site =",site + endif + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + openw, Unit10, workpath+'cal_stack_cub_final.txt', /GET_LUN + openw, Unit11, workpath+'cal_stack_cub_sep.txt', /GET_LUN + endif + + if strmatch(sublist[0],'*TPlike*.*',/FOLD_CASE) then begin + TPlike=1 ; this is a SARDARA file converted/integrated in TotalPower-like format + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + printf, Unit0, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 c2J_0 c2J_1 e_c2J_0 e_c2J_1' + printf, Unit0, ' MHz deg count count count count Jy/cnt Jy/cnt Jy/cnt Jy/cnt' + printf, Unit1, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 c2J_0 c2J_1 e_c2J_0 e_c2J_1 SNR_0 SNR_1 Offset FW/HP Type' + printf, Unit1, ' MHz deg count count count count Jy/cnt Jy/cnt Jy/cnt Jy/cnt deg ' + endif + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + printf, Unit10, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 c2J_0 c2J_1 e_c2J_0 e_c2J_1' + printf, Unit10, ' MHz deg count count count count Jy/cnt Jy/cnt Jy/cnt Jy/cnt' + printf, Unit11, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 c2J_0 c2J_1 e_c2J_0 e_c2J_1 SNR_0 SNR_1 Offset FW/HP Type' + printf, Unit11, ' MHz deg count count count count Jy/cnt Jy/cnt Jy/cnt Jy/cnt deg ' + endif + endif else begin + TPlike=2 ; this is an original TotalPower FITS + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + printf, Unit0, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 c2J_0 c2J_1 e_c2J_0 e_c2J_1' + printf, Unit0, ' MHz deg count count count count Jy/cnt Jy/cnt Jy/cnt Jy/cnt' + printf, Unit1, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 c2J_0 c2J_1 e_c2J_0 e_c2J_1 SNR_0 SNR_1 Offset FW/HP Type' + printf, Unit1, ' MHz deg count count count count Jy/cnt Jy/cnt Jy/cnt Jy/cnt deg ' + endif + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + printf, Unit10, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 c2J_0 c2J_1 e_c2J_0 e_c2J_1' + printf, Unit10, ' MHz deg count count count count Jy/cnt Jy/cnt Jy/cnt Jy/cnt' + printf, Unit11, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 c2J_0 c2J_1 e_c2J_0 e_c2J_1 SNR_0 SNR_1 Offset FW/HP Type' + printf, Unit11, ' MHz deg count count count count Jy/cnt Jy/cnt Jy/cnt Jy/cnt deg ' + endif + endelse + + + + openw, Unit2, workpath+'cal_prints.txt', /GET_LUN + + + ; reads skydip information + ans=file_test(skydippath+'Tau.txt') + if ans eq 1 then begin + readcol, skydippath+'Tau.txt', tau_mjd, tau_freq, t30l, t90l, t30r, t90r, tau_empL, tau_empR, tau0L, tau0R, format='(d,d,d,d,d,d,d,d,d,d)', /SILENT + for t=0,n_elements(tau_mjd)-1 do begin + ; if for any reason tau_emp is considered more reliable than tau0, uncomment the following two lines + ; tau0L[t]=tau_empL[t] + ; tau0R[t]=tau_empR[t] + endfor + caltau = 1 + endif else begin + print, ' ' + print, '+++++++++++++++++++++++++++++++++++++++++++++' + print, 'Beware: no Tau.txt available, assuming tau0=0 + print, '+++++++++++++++++++++++++++++++++++++++++++++' + wait, 0.5 + tau0L=[0,0] + tau0R=[0,0] + caltau = 0 + endelse + + ; definition of arrays + p0=fltArr(2) + p1=fltArr(2) + e0=fltArr(2) + e1=fltArr(2) + c0=fltArr(2) ; count to jansky array + c1=fltArr(2) + d0=fltArr(2) ; error on count to jansky array + d1=fltArr(2) + SNR0=fltArr(2) + SNR1=fltArr(2) + offset=dblArr(2) + p0c=fltArr(2) + p1c=fltArr(2) + e0c=fltArr(2) + e1c=fltArr(2) + c0c=fltArr(2) ; count to jansky array + c1c=fltArr(2) + d0c=fltArr(2) ; error on count to jansky array + d1c=fltArr(2) + SNR0c=fltArr(2) + SNR1c=fltArr(2) + offsetc=dblArr(2) + fw_ratio0=fltarr(2) + fw_ratio1=fltarr(2) + fw_ratio0c=fltarr(2) + fw_ratio1c=fltarr(2) + + + ; scan-based analysis + for i =0,number-1 do begin + print, ' ' + print, '=========================================' + print, 'Now analysing scan ', list[i] + print, '=========================================' + sommaL=0 + sommaR=0 + sommaLA = 0 + sommaRA = 0 + numLA = 0 + numRA = 0 + iiA = 0 + ffA = 0 + sommaLB = 0 + sommaRB = 0 + numLB = 0 + numRB = 0 + iiB = 0 + ffB = 0 + + splitfoldname1=strsplit(list[i],sep,/extract) + foldname=splitfoldname1[-1] + splitfoldname2=strsplit(foldname,'-',/extract) + hhmmss = splitfoldname2[1] + sublist = FILE_SEARCH(list[i]+sep+'2*.fits', COUNT=subnumber) + subscan=file_basename(sublist) + + ; reads flagging info + err=0 + checkfile = list[i]+sep+'checkfile_'+foldname+'.txt' + openr, Unit9, checkfile, ERROR=err, /GET_LUN + flagcode=strarr(subnumber) + flagcode[*]='+1cx' + if (err lt 0) then begin + print, "No valid flagging information for this scan" + print, "Temporarily assuming all scans are good" + endif else begin + readcol, checkfile, F='a,a,a', subscan_t, flagcode_t, who, /SILENT + if (n_elements(subscan_t) lt subnumber) then begin + print, "Some subscans miss flagging information" + print, "Temporarily assuming those subscans are good" + endif + for n_flag = 0, n_elements(subscan_t)-1 do begin + for j_flag = 0, subnumber-1 do begin + if (subscan_t[n_flag] eq subscan[j_flag]) then begin + flagcode[j_flag]=flagcode_t[n_flag] + endif + endfor + endfor + close, Unit9 + free_lun, Unit9 + endelse + + + if dosingle eq 'y' then begin + ; creating output files for single-subscan analysis + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + openw, Unit3, list[i]+sep+'cal_single_linear.txt', /GET_LUN + printf, Unit3,' Name # MJD El(r) El(d) cnt_0 cnt_1 e_c_0 e_c_1 offL(d) offR(d) FW/HP scantype' + endif + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + openw, Unit4, list[i]+sep+'cal_single_cubic.txt', /GET_LUN + printf, Unit4,' Name # MJD El(r) El(d) cnt_0 cnt_1 e_c_0 e_c_1 offL(d) offR(d) FW/HP scantype' + endif + openw, Unit5, list[i]+sep+'cal_prints.txt', /GET_LUN + endif + + ; Dynamic arrays (for amplitudes and errors, for each section) + + Lampl=[ ] + Rampl=[ ] + Lerrs=[ ] + Rerrs=[ ] + offset_sub=[ ] + type=[ ] + mean_lat=[ ] + + ; subscan-based measurements + for j=0,subnumber-1 do begin + ; READ DATA TABLE (binary data table of the MegaTable FITS file) + data=MRDFITS(sublist[j],4,/SILENT) + RFinfo=MRDFITS(sublist[j],2,/SILENT) + Sectinfo=MRDFITS(sublist[j],1,/SILENT) + ; read frequency and bandwith from RF table + bw=RFinfo[0].bandWidth + freq=RFinfo[0].frequency+bw/2.0 + dt=1.0/(Sectinfo[0].sampleRate*1e6) + + ndat = (size(data.time)) + + ; takes source name and subscantype from header + mainh=mrdfits(sublist[j],0,head0,/silent) + + ; assessing whether this is an old converted FITS + findh0key, head0, '*Converted FITS*', keylab, converted, infolab, convertedflag + findh0key, head0, '*Obtained from ESCS 0.1 TP FITS*', keylab, version, infolab, oldestflag + + if convertedflag eq 0 then begin ; this IS NOT a converted file: ESCS 0.3 keyword names hold + findh0key, head0, '*SOURCE*', keylab, target, infolab, tarflag + findh0key, head0, '*SubScanType =*', keylab, orscantype, infolab, tarflag + findh0key, head0, '*Azimuth Offset =*', keylab, hAZoff, infolab, typeflag ; offset utente e di sistema + findh0key, head0, '*Elevation Offset =*', keylab, hELoff, infolab, typeflag ; offset utente e di sistema + findh0key, head0, '*RightAscension Offset =*', keylab, hRAoff, infolab, typeflag ; offset utente e di sistema + findh0key, head0, '*Declination Offset =*', keylab, hDECoff, infolab, typeflag ; offset utente e di sistema + findh0key, head0, '*RightAscension =*', keylab, ras, infolab, typeflag ; RA nominale del target, radianti + findh0key, head0, '*Declination =*', keylab, decs, infolab, typeflag ; DEC nominale del target, radianti + endif else begin ; this IS a converted file: custom 8-char keyword names hold according to original FITS date, and offsets are not present + findh0key, head0, '*SOURCE*', keylab, target, infolab, tarflag + findh0key, head0, '*SUBSCANT=*', keylab, orscantype, infolab, tarflag + hAZoff=0.0d + hELoff=0.0d + hRAoff=0.0d + hDECoff=0.0d + if oldestflag eq 0 then begin + findh0key, head0, '*RightAscension =*', keylab, ras, infolab, typeflag ; RA nominale del target, radianti + findh0key, head0, '*Declination =*', keylab, decs, infolab, typeflag ; DEC nominale del target, radianti + endif else begin ; these are the oldest FITS, where several more keywords have converted names + findh0key, head0, '*RIGHTASC=*', keylab, ras, infolab, typeflag ; RA nominale del target, radianti + findh0key, head0, '*DECLINAT=*', keylab, decs, infolab, typeflag ; DEC nominale del target, radianti + endelse + endelse + + + ; assigning labels and coefficients to calibrators + target=strcompress(target, /remove_all) + cleantarget=strsplit(target,"'",/extract) + target=cleantarget[0] + sourcename=target + method='none' + case strupcase(target) of + '3C48': begin + calflux, strupcase(target), freq/1000.0, intSnu, intSnu_err + method='P&B' + end + '3C123': begin + ; here follows a patch made necessary by a schedule mistake: 3C48 was observed instead of 3C123, using 3C123 as target name! + if abs(ras*180.0/!dpi - 24.4221) lt 0.1 and abs(decs*180.0/!dpi - 33.1598) lt 0.1 then begin + target='3C48' + sourcename=target + endif + calflux, strupcase(target), freq/1000.0, intSnu, intSnu_err + method='P&B' + end + '3C147': begin + calflux, strupcase(target), freq/1000.0, intSnu, intSnu_err + method='P&B' + end + '3C161': begin + a=1.250 + b=+0.726 + c=-0.2286 + low=1408 + high=10550 + flux=10^(a+b*alog10(freq)+c*alog10(freq)*alog10(freq)) + method='Ott' + end + '3C286': begin + calflux, strupcase(target), freq/1000.0, intSnu, intSnu_err + method='P&B' + end + '3C295': begin + calflux, strupcase(target), freq/1000.0, intSnu, intSnu_err + method='P&B' + end + 'NGC7027': begin + calflux, strupcase(target), freq/1000.0, intSnu, intSnu_err + method='P&B' + end + 'DR21': begin + a=1.810 + b=-0.122 + c=0.0 + low=7000 + high=31000 + flux=10^(a+b*alog10(freq)+c*alog10(freq)*alog10(freq)) + method='Ott' + end + 'MYSRC': begin + ; assuming flux density = 1 Jy + a=0.00 + b=0.00 + c=0.00 + flux=10^(a+b*alog10(freq)+c*alog10(freq)*alog10(freq)) + method='CustomSRC' + end + 'CAL01': begin + ; fake source for tests, having flux density=0.1 Jy + a=-1.00 + b=0.00 + c=0.00 + flux=10^(a+b*alog10(freq)+c*alog10(freq)*alog10(freq)) + method='CustomSRC' + end + 'CAL05': begin + ; fake source for tests, having flux density=0.5 Jy + a=-0.301030 + b=0.00 + c=0.00 + flux=10^(a+b*alog10(freq)+c*alog10(freq)*alog10(freq)) + method='CustomSRC' + end + 'CAL1': begin + ; fake source for tests, having flux density=1 Jy + a=0.00 + b=0.00 + c=0.00 + flux=10^(a+b*alog10(freq)+c*alog10(freq)*alog10(freq)) + method='CustomSRC' + end + 'CAL2': begin + ; fake source for tests, having flux density=2 Jy + a=0.30103 + b=0.00 + c=0.00 + flux=10^(a+b*alog10(freq)+c*alog10(freq)*alog10(freq)) + method='CustomSRC' + end + 'CAL3': begin + ; fake source for tests, having flux density=3 Jy + a=0.4771212 + b=0.00 + c=0.00 + flux=10^(a+b*alog10(freq)+c*alog10(freq)*alog10(freq)) + method='CustomSRC' + end + else: begin + ; reading external list + if myextlist ne ' ' then begin + readcol, myextlist, extcalname, extcalflux, format='(A,D)', /SILENT + foundit=where(strmatch(extcalname,target,/FOLD_CASE) eq 1) + if foundit eq -1 then begin + print, ' ' + print, target+' is an unknown calibrator, cannot retrieve a flux density to use:' + print, 'it will return dummy values in the measurements' + print, ' ' + a=4.00 + b=0.00 + c=0.00 + flux=10^(a+b*alog10(freq)+c*alog10(freq)*alog10(freq)) + method='CustomSRC' + endif else begin + a=alog10(extcalflux[foundit]) ; fixed value, as the external list does not contain polynomials + b=0.00 + c=0.00 + flux=10^(a+b*alog10(freq)+c*alog10(freq)*alog10(freq)) + method='CustomSRC' + endelse + endif else begin + print, ' ' + print, target+' is an unknown calibrator, cannot retrieve a flux density to use:' + print, 'it will return dummy values in the measurements' + print, ' ' + a=0.00 + b=0.00 + flux=10^(a+b*alog10(freq)+c*alog10(freq)*alog10(freq)) + method='CustomSRC' + endelse + endelse + endcase + + ; aligning the P&B estimates to the general nomenclature + if method eq 'P&B' then flux=intSnu + + if j eq 0 then begin + ; expliciting the flux density computation result, only for the first subscan + print, ' ' + print, '***********' + print, sourcename, 'Flux density', flux, ' Jy at', freq/1000.0, 'GHz', format='(A,1X,A,1X,D6.3,1X,A,1X,D5.2,1X,A3)' + print, 'Estimate by', method, format='(A,1X,A)' + print, '***********' + endif + + ; nominal target coordinates, converted in degrees + rasd=ras*180.0d/!dpi + decsd=decs*180.0d/!dpi + + elapsed=(double(data.time)-double(data[0].time))*24.0*3600.0 ; elapsed time from first sample + duration=dt*ndat[1]/60.0 ; duration of the whole subscan (approx) + ; RA and DEC span of the subscan + deltaDec=abs(data[0].decj2000-data[-1].decj2000)/!dpi*180.0 + deltaRA=abs(data[0].raj2000-data[-1].raj2000)*cos(mean(data.decj2000))/!dpi*180.0 + + ; determine speed (deg/min) along scan direction + ; if type is AZ or EL, treats is as DEC or RA depending on which quantity spans a larger range + splitscantype=strsplit(orscantype,"'",/extract) + cleanscantype=strcompress(splitscantype[1], /remove_all) ; cleaning the string so that it become more manageable + case cleanscantype of + 'DEC': begin + if (data[0].decj2000 lt data[-1].decj2000) then scantype = 1 else scantype = 2 + speed=deltaDec/duration + end + 'RA': begin + if (data[0].raj2000 lt data[-1].raj2000) then scantype = 3 else scantype = 4 + speed=deltaRA/duration + end + 'AZ': begin + if (deltaDec gt deltaRA) then begin + if (data[0].decj2000 lt data[-1].decj2000) then scantype = 1 else scantype = 2 + speed=deltaDec/duration + endif + if (deltaDec lt deltaRA) then begin + if (data[0].raj2000 lt data[-1].raj2000) then scantype = 3 else scantype = 4 + speed=deltaRA/duration + endif + end + 'EL': begin + if (deltaDec gt deltaRA) then begin + if (data[0].decj2000 lt data[-1].decj2000) then scantype = 1 else scantype = 2 + speed=deltaDec/duration + endif + if (deltaDec lt deltaRA) then begin + if (data[0].raj2000 lt data[-1].raj2000) then scantype = 3 else scantype = 4 + speed=deltaRA/duration + endif + end + 'UNKNOWN': begin + if (deltaDec gt deltaRA) then begin + if (data[0].decj2000 lt data[-1].decj2000) then scantype = 1 else scantype = 2 + speed=deltaDec/duration + endif + if (deltaDec lt deltaRA) then begin + if (data[0].raj2000 lt data[-1].raj2000) then scantype = 3 else scantype = 4 + speed=deltaRA/duration + endif + end + endcase + + ; check for spikes in X-band wide bandwidth data XXX Perché farlo solo per la banda X? + if (bw gt 500 and freq lt 10000) then begin + for k=3, ndat[1]-4 do begin + temp0=[data[k-3].ch0,data[k-2].ch0,data[k-1].ch0,data[k+1].ch0,data[k+2].ch0,data[k+3].ch0] + temp1=[data[k-3].ch1,data[k-2].ch1,data[k-1].ch1,data[k+1].ch1,data[k+2].ch1,data[k+3].ch1] + avech0=mean(temp0) + avech1=mean(temp1) + stdch0=stddev(temp0) + stdch1=stddev(temp1) + if ((data[k].ch0-avech0) gt 7*stdch0) then begin + data[k].ch0=avech0 + endif + if ((data[k].ch1-avech1) gt 7*stdch1) then begin + data[k].ch1=avech1 + endif + endfor + endif + + ; reverse arrays for scans in decreasing RA or Dec + if (scantype eq 2 or scantype eq 4) then begin + ; data=reverse(data) + data.raj2000=reverse(data.raj2000) + data.decj2000=reverse(data.decj2000) + data.ch0=reverse(data.ch0) ; original line + ; data.ch0=reverse(data.ch2) ; ABC + data.ch1=reverse(data.ch1) + endif + + ; computing the sample offset with respect to the source position + xdec=(data.decj2000-decs)/!dpi*180.0 + xra=(data.raj2000-ras)*cos(mean(data.decj2000))/!dpi*180.0 + + scanname = sublist[j] + case flagcode[j] of + '-1c0': begin + k0=1 + k1=0 + end + '-1c1': begin + k0=0 + k1=1 + end + '+0cx': begin + k0=1 + k1=1 + end + '+1cx': begin + k0=1 + k1=1 + end + '-1cx': begin + k0=0 + k1=0 + end + endcase + + ; shift scans of type 'Dec' so they align with first one; sets start and end value of array with data from every scan + if (scantype eq 1 or scantype eq 2) then begin + if ((numLA+numRA) eq 0) then begin + dataA=data + ndatA=ndat + xdecA=xdec + iiA=0 + if (k0+k1 gt 0) then ffA=ndatA[1]-1 + decoffA = min(xdecA,xdec0A,/ABSOLUTE) ; xdec0A is the number of sample at the source declination for scan A + endif else begin + decoff = min(xdec,xdec0,/ABSOLUTE) + ; shift data by (xdec0-xdec0B) samples + ndelta=xdec0-xdec0A + if (ndelta gt 0) then begin + iiA=iiA + if (k0+k1 gt 0) then ffA=min([ffA,ndat[1]-1-ndelta]) + for k=iiA,ffA do begin + if (k0+k1 gt 0) then data[k]=data[k+ndelta] + endfor + endif + if (ndelta lt 0) then begin + if (k0+k1 gt 0) then iiA=max([iiA,0-ndelta]) + ; here we have a problem: we would like to start from ffA and go down to iiA, since ffA can be filled by shifting "data" right. but data only has dimension ndat, not ndat-ndelta (ndelta is <0) + ; ffA=min([ffA,ndat[1]-1-ndelta]) + if (k0+k1 gt 0) then ffA=min([ffA,ndat[1]-1]) + for k=ffA,iiA,-1 do begin + if (k0+k1 gt 0) then data[k]=data[k+ndelta] + endfor + endif + if (ndelta eq 0) then begin + if (k0+k1 gt 0) then ffA=min([ffA,ndat[1]-1]) + endif + endelse + endif + + ; shifts scans of type 'RA' so they align with first one; sets start and end value of array with data from every scan + if (scantype eq 3 or scantype eq 4) then begin + if ((numLB+numRB) eq 0) then begin + dataB=data + ndatB=ndat + xraB=xra + iiB=0 + if (k0+k1 gt 0) then ffB=ndatB[1]-1 + raoffB = min(xraB,xra0B,/ABSOLUTE) ; xra0B is the number of sample at the source right ascension for scan B + endif else begin + raoff = min(xra,xra0,/ABSOLUTE) + ; shift data by (xra0-xra0B) samples + ndelta=xra0-xra0B + if (ndelta gt 0) then begin + iiB=iiB + if (k0+k1 gt 0) then ffB=min([ffB,ndat[1]-1-ndelta]) + for k=iiB,ffB do begin + if (k0+k1 gt 0) then data[k]=data[k+ndelta] + endfor + endif + if (ndelta lt 0) then begin + if (k0+k1 gt 0) then iiB=max([iiB,0-ndelta]) + ; here we have a problem: we would like to start from ffB and go down to iiB, since ffB can be filled by shifting "data" right. but data only has dimension ndat, not ndat-ndelta (ndelta is <0) + ; ffB=min([ffB,ndat[1]-1-ndelta]) + if (k0+k1 gt 0) then ffB=min([ffB,ndat[1]-1]) + for k=ffB,iiB,-1 do begin + if (k0+k1 gt 0) then data[k]=data[k+ndelta] + endfor + endif + if (ndelta eq 0) then begin + if (k0+k1 gt 0) then ffB=min([ffB,ndat[1]-1]) + endif + endelse + endif + + time_s = (data.time-data[round(ndat[1]/2)].time)*24.d0*3600.d0 + + midscan=FIX(ndat[1]/2) + beamspan= beam/speed ; duration of acquisition for one beamwidth + Nsamples=ROUND(beamspan/dt) ; number of samples in one beamwidth + + el_d=data[midscan].el*180.0/!dpi + + g_0=A_0[0]+A_0[1]*el_d+A_0[2]*el_d^2+A_0[3]*el_d^3 + g_1=A_1[0]+A_1[1]*el_d+A_1[2]*el_d^2+A_1[3]*el_d^3 + + if keyword_set(skipgc) then begin + ; Overriding the gain computation and applying a flat gain curve (constant value = 1.0) + g_0=1.0 + g_1=1.0 + endif + + ; actually sums counts only for scans flagged as good and sets inputs for single scan fit + if ( scantype eq 1 or scantype eq 2) then begin + ; if (k0 ne 0) then sommaLA = sommaLA + k0*data.ch2/g_0 ; ABC + if (k0 ne 0) then sommaLA = sommaLA + k0*data.ch0/g_0 ; original line + if (k1 ne 0) then sommaRA = sommaRA + k1*data.ch1/g_1 + numLA = numLA + k0*1 + numRA = numRA + k1*1 + ascissa=(data.decj2000-hDECoff)*180d/!dpi + gpos=decsd + endif else begin + ; if (k0 ne 0) then sommaLB = sommaLB + k0*data.ch2/g_0 ; ABC + if (k0 ne 0) then sommaLB = sommaLB + k0*data.ch0/g_0 ; original line + if (k1 ne 0) then sommaRB = sommaRB + k1*data.ch1/g_1 + numLB = numLB + k0*1 + numRB = numRB + k1*1 + ascissa=(data.raj2000-hRAoff/cos(mean(data.decj2000)))*180d/!dpi + gpos=rasd + endelse + + + if dosingle eq 'y' then begin + + if (fitchoice eq 'linear') or (fitchoice eq 'both') then begin + ; fit single subscans with linear + gaussian + Lfwhm=0 + Rfwhm=0 + ; yy0=data.ch2/g_0 ; ABC + yy0=data.ch0/g_0 ; original line + yy1=data.ch1/g_1 + if (scantype eq 1 or scantype eq 2) then tipo=0 else tipo=1 + calibfit, k0,'single','linear','Ch_0',tipo,list[i],subscan[j],Unit5,flux,data[0].el,tau0L,ascissa,yy0,0,ndat[1]-1,midscan,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Lfwhm, nL, offL, peak_cnt_0, err_cnt_0, dum3, dum4, dum5, psingle, doplot + calibfit, k1,'single','linear','Ch_1',tipo,list[i],subscan[j],Unit5,flux,data[0].el,tau0R,ascissa,yy1,0,ndat[1]-1,midscan,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Rfwhm, nR, offR, peak_cnt_1, err_cnt_1, dum3, dum4, dum5, psingle, doplot + if doplot eq 'y' then begin + splitscan=strsplit(subscan[j],'.',/extract) + scanroot=splitscan[0] + if (j eq subnumber-1) then begin + ; psingle.save, list[i]+'_lin_single.pdf', /LANDSCAPE, RESOLUTION=300, /APPEND, /CLOSE + endif else begin + ; psingle.save, list[i]+sep+scanroot+'_single.pdf', /LANDSCAPE, RESOLUTION=300, /APPEND + endelse + endif + if TPlike eq 2 then begin + printf, Unit3, FORMAT = '(a15,D12.5,9(" ",D08.4),i)', sourcename, data[midscan].time,$ + data[midscan].el, el_d, peak_cnt_0, peak_cnt_1, err_cnt_0, err_cnt_1, offL, offR, Lfwhm/beam, scantype + endif else begin + printf, Unit3, FORMAT = '(a15,1X,D12.5,1X,D08.4,1X,D08.4,1X,E12.6,1X,E12.6,1X,E12.6,1X,E12.6,1X,D08.4,1X,D08.4,1X,D08.4,1X,i)', sourcename, data[midscan].time,$ + data[midscan].el, el_d, peak_cnt_0, peak_cnt_1, err_cnt_0, err_cnt_1, offL, offR, Lfwhm/beam, scantype + endelse + endif + + if (fitchoice eq 'cubic') or (fitchoice eq 'both') then begin + ; fit single subscans with cubic + gaussian + Lfwhm=0 + Rfwhm=0 + ; yy0=data.ch2/g_0 ;ABC + yy0=data.ch0/g_0 ; original line + yy1=data.ch1/g_1 + if (scantype eq 1 or scantype eq 2) then tipo=0 else tipo=1 + calibfit, k0,'single','cubic','Ch_0',tipo,list[i],subscan[j],Unit5,flux,data[0].el,tau0L,ascissa,yy0,0,ndat[1]-1,midscan,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Lfwhm, nL, offL, peak_cnt_0, err_cnt_0, dum3, dum4, dum5, psingle, doplot + calibfit, k1,'single','cubic','Ch_1',tipo,list[i],subscan[j],Unit5,flux,data[0].el,tau0R,ascissa,yy1,0,ndat[1]-1,midscan,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Rfwhm, nR, offR, peak_cnt_1, err_cnt_1, dum3, dum4, dum5, psingle, doplot + ; NOTA: modificata la gestione del caso /sub; ora il salvataggio del PDF avviene direttamente dentro alla procedura dataplot, + ; che salva un file per ogni chiamata. + + if TPlike eq 2 then begin + printf, Unit4, FORMAT = '(a15,D12.5,9(" ",D08.4),i)', sourcename, data[midscan].time,$ + data[midscan].el, el_d, peak_cnt_0, peak_cnt_1, err_cnt_0, err_cnt_1, offL, offR, Lfwhm/beam, scantype + endif else begin + printf, Unit4, FORMAT = '(a15,1X,D12.5,1X,D08.4,1X,D08.4,1X,E12.6,1X,E12.6,1X,E12.6,1X,E12.6,1X,D08.4,1X,D08.4,1X,D08.4,1X,i)', sourcename, data[midscan].time,$ + data[midscan].el, el_d, peak_cnt_0, peak_cnt_1, err_cnt_0, err_cnt_1, offL, offR, Lfwhm/beam, scantype + endelse + endif + + if (nL+nR gt 0) then offset_sub = [offset_sub,(offL+offR)/(nL+nR)] + + endif + + endfor + + if dosingle eq 'y' then begin + ; closing the output files and freeing memory + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + close, Unit3 + FREE_LUN, Unit3 + endif + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + close, Unit4 + FREE_LUN, Unit4 + endif + close, Unit5 + free_lun, Unit5, /force + endif + + ; stacked analysis begins + + for tipo = 0,1 do begin + ascissa=time_s ; default if all subscans are void + + ; set some inputs + case tipo of + 0: begin + if (numLA+numRA gt 0) then begin + if (numLA gt 0) then sommaL = sommaLA/numLA else sommaL = sommaLA + if (numRA gt 0) then sommaR = sommaRA/numRA else sommaR = sommaRA + numL = numLA + numR = numRA + ii=iiA + ff=ffA + ; tipo=0 means scan declination, ie scantype Dec, ie constant RA + ascissa=(dataA.decj2000-hDECoff)*180d/!dpi + xdiffA = min(ascissa-decsd,x_mid,/ABSOLUTE) ; x_mid contains the sample number nearest to the source nominal position + gpos=decsd + endif else begin + numL=0 + numR=0 + endelse + end + 1: begin + if (numLB+numRB gt 0) then begin + if (numLB gt 0) then sommaL = sommaLB/numLB else sommaL = sommaLB + if (numRB gt 0) then sommaR = sommaRB/numRB else sommaR = sommaRB + numL = numLB + numR = numRB + ii=iiA + ff=ffB + ; tipo=1 means scan right asc., ie scantype RA, ie constant Dec + ascissa=(dataB.raj2000-hRAoff/cos(mean(dataB.decj2000)))*180d/!dpi + xdiffB = min(ascissa-rasd,x_mid,/ABSOLUTE) ; x_mid contains the sample number nearest to the source nominal position + gpos=rasd + endif else begin + numL=0 + numR=0 + endelse + end + endcase + + PRINTF, Unit2, " " + PRINTF, Unit2, "****** "+scanname+" ******" + PRINTF, Unit2, FORMAT = '("* MJD = ", D014.8)', data[midscan].time + PRINTF, Unit2, "* nsamples = ", ndat[1] + PRINTF, Unit2, "* Beamspan = ", beamspan, " [s]" + PRINTF, Unit2, "* Beam samples = ", Nsamples + + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + + ; LINEAR FIT OF STACKED SCAN + Lfwhm=0 + Rfwhm=0 + + calibfit, numL,'stacked','linear','Ch_0',tipo,list[i],list[i],Unit2,flux,data[0].el,tau0L,ascissa,sommaL,ii,ff,x_mid,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Lfwhm, n_offL, offL, dum1, dum2, dum3, dum4, dum5, pstack, doplot + p0[tipo]=dum1 + e0[tipo]=dum2 + c0[tipo]=dum3 + d0[tipo]=dum4 + SNR0[tipo]=dum5 + fw_ratio0[tipo]=Lfwhm/beam + + calibfit, numR,'stacked','linear','Ch_1',tipo,list[i],list[i],Unit2,flux,data[0].el,tau0R,ascissa,sommaR,ii,ff,x_mid,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Rfwhm, n_offR, offR, dum1, dum2, dum3, dum4, dum5, pstack, doplot + p1[tipo]=dum1 + e1[tipo]=dum2 + c1[tipo]=dum3 + d1[tipo]=dum4 + SNR1[tipo]=dum5 + fw_ratio1[tipo]=Rfwhm/beam + + if doplot eq 'y' then begin + if (tipo eq 1 and fitchoice eq 'linear') then begin + pstack.save, list[i]+'_stacked.pdf', /LANDSCAPE, RESOLUTION=300, /APPEND, /CLOSE + endif else begin + pstack.save, list[i]+'_stacked.pdf', /LANDSCAPE, RESOLUTION=300, /APPEND + endelse + endif + + if (n_offL+n_offR gt 0) then offset[tipo]=(offL+offR)/(n_offL+n_offR) else offset[tipo] = -99. + + el_d=data[0].el*180.0/!dpi + + if tipo eq 0 then sclab='DEC' + if tipo eq 1 then sclab='RA' + + if TPlike eq 2 then begin + printf, Unit1, FORMAT = '(a15,a7,D12.5,d8.1,13(" ",D8.4),A7)', sourcename, hhmmss, data[0].time, freq, $ + el_d, p0[tipo],p1[tipo],e0[tipo],e1[tipo],c0[tipo],c1[tipo],d0[tipo],d1[tipo],SNR0[tipo], SNR1[tipo], offset[tipo], fw_ratio0[tipo], sclab + endif else begin + printf, Unit1, FORMAT = '(a15,a7,D12.5,d8.1,d8.4,8(" ",E12.6),4(" ",D8.4),A7)', sourcename, hhmmss, data[0].time, freq, $ + el_d, p0[tipo],p1[tipo],e0[tipo],e1[tipo],c0[tipo],c1[tipo],d0[tipo],d1[tipo],SNR0[tipo], SNR1[tipo], offset[tipo], fw_ratio0[tipo], sclab + endelse + + endif + + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + + ; CUBIC FIT OF STACKED SCAN + Lfwhm=0 + Rfwhm=0 + + calibfit, numL,'stacked','cubic','Ch_0',tipo,list[i],list[i],Unit2,flux,data[0].el,tau0L,ascissa,sommaL,ii,ff,x_mid,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Lfwhm, n_offL, offL, dum1, dum2, dum3, dum4, dum5, pstack, doplot + p0c[tipo]=dum1 + e0c[tipo]=dum2 + c0c[tipo]=dum3 + d0c[tipo]=dum4 + SNR0c[tipo]=dum5 + fw_ratio0c[tipo]=Lfwhm/beam + + calibfit, numR,'stacked','cubic','Ch_1',tipo,list[i],list[i],Unit2,flux,data[0].el,tau0R,ascissa,sommaR,ii,ff,x_mid,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Rfwhm, n_offR, offR, dum1, dum2, dum3, dum4, dum5, pstack, doplot + p1c[tipo]=dum1 + e1c[tipo]=dum2 + c1c[tipo]=dum3 + d1c[tipo]=dum4 + SNR1c[tipo]=dum5 + fw_ratio1c[tipo]=Rfwhm/beam + + if doplot eq 'y' then begin + if (tipo eq 1) then begin + pstack.save, list[i]+'_stacked.pdf', /LANDSCAPE, RESOLUTION=300, /APPEND, /CLOSE + endif else begin + pstack.save, list[i]+'_stacked.pdf', /LANDSCAPE, RESOLUTION=300, /APPEND + endelse + endif + + if (n_offL+n_offR gt 0) then offsetc[tipo]=(offL+offR)/(n_offL+n_offR) else offsetc[tipo] = -99. + + el_d=data[0].el*180.0/!dpi + + if tipo eq 0 then sclab='DEC' + if tipo eq 1 then sclab='RA' + + if TPlike eq 2 then begin + printf, Unit11, FORMAT = '(a15,a7,D12.5,d8.1,13(" ",D8.4),A7)', sourcename, hhmmss, data[0].time, freq, $ + el_d, p0c[tipo],p1c[tipo],e0c[tipo],e1c[tipo],c0c[tipo],c1c[tipo],d0c[tipo],d1c[tipo],SNR0c[tipo], SNR1c[tipo], offsetc[tipo], Lfwhm/beam, sclab + endif else begin + printf, Unit11, FORMAT = '(a15,a7,D12.5,d8.1,d8.4,8(" ",E12.6),4(" ",D8.4),A7)', sourcename, hhmmss, data[0].time, freq, $ + el_d, p0c[tipo],p1c[tipo],e0c[tipo],e1c[tipo],c0c[tipo],c1c[tipo],d0c[tipo],d1c[tipo],SNR0c[tipo], SNR1c[tipo], offsetc[tipo], Lfwhm/beam, sclab + endelse + endif + endfor + + ; If the pointing offsets are available in both directions, the amplitude measurement is compensated for them. + ; Otherwise, the measurement is rejected. + ; In this position also the FWHM check against HPBW - i.e. beam - is performed. A tolerance of FWHM/beam is here defined, the + ; measurements outside the consequent boundaries are rejected as well. + + fw_tol=0.2 ; we accept measurements where FWHM is within 20% of the nominal HPBW + + if (fitchoice eq 'linear') or (fitchoice eq 'both') then begin + + if (offset[0] eq -99. or offset[1] eq -99. or (abs(fw_ratio0[0]-1.0) ge fw_tol and abs(fw_ratio0[1]-1.0) ge fw_tol) or (abs(fw_ratio1[0]-1.0) ge fw_tol and abs(fw_ratio1[1]-1.0) ge fw_tol)) then begin + peak_cnt_0=-99.0 + peak_cnt_1=-99.0 + err_cnt_0=-99.0 + err_cnt_1=-99.0 + cnt2Jy_0=-99.0 + cnt2Jy_1=-99.0 + err_cnt2Jy_0=-99.0 + err_cnt2Jy_1=-99.0 + endif else begin + ; riscalo le ampiezze e gli errori della polarizzazione 0 per l'offset incrociato + p0[0]=p0[0]/exp(-(offset[1]*1.66/beamd)^2.) + p0[1]=p0[1]/exp(-(offset[0]*1.66/beamd)^2.) + e0[0]=e0[0]/exp(-(offset[1]*1.66/beamd)^2.) + e0[1]=e0[1]/exp(-(offset[0]*1.66/beamd)^2.) + ; riscalo le ampiezze e gli errori della polarizzazione 1 per l'offset incrociato + p1[0]=p1[0]/exp(-(offset[1]*1.66/beamd)^2.) + p1[1]=p1[1]/exp(-(offset[0]*1.66/beamd)^2.) + e1[0]=e1[0]/exp(-(offset[1]*1.66/beamd)^2.) + e1[1]=e1[1]/exp(-(offset[0]*1.66/beamd)^2.) + + peak_cnt_0 = cal_wmean(p0, e0, /nan) + peak_cnt_1 = cal_wmean(p1, e1, /nan) + w0=1/(e0/p0) + w1=1/(e1/p1) + err_cnt_0 = sqrt((w0[0]*e0[0])^2+(w0[1]*e0[1])^2)/(w0[0]+w0[1]) + err_cnt_1 = sqrt((w1[0]*e1[0])^2+(w1[1]*e1[1])^2)/(w1[0]+w1[1]) + ; err_cnt_0 = stddev(p0, /nan) + ; err_cnt_1 = stddev(p1, /nan) + + ; riscalo i ct2Jy e gli errori della polarizzazione 0 per l'offset incrociato + c0[0]=c0[0]*exp(-(offset[1]*1.66/beamd)^2.) + c0[1]=c0[1]*exp(-(offset[0]*1.66/beamd)^2.) + d0[0]=d0[0]*exp(-(offset[1]*1.66/beamd)^2.) + d0[1]=d0[1]*exp(-(offset[0]*1.66/beamd)^2.) + ; riscalo i ct2Jy e gli errori della polarizzazione 1 per l'offset incrociato + c1[0]=c1[0]*exp(-(offset[1]*1.66/beamd)^2.) + c1[1]=c1[1]*exp(-(offset[0]*1.66/beamd)^2.) + d1[0]=d1[0]*exp(-(offset[1]*1.66/beamd)^2.) + d1[1]=d1[1]*exp(-(offset[0]*1.66/beamd)^2.) + + cnt2Jy_0 = cal_wmean(c0, d0, /nan) + cnt2Jy_1 = cal_wmean(c1, d1, /nan) + w0=1/(d0/c0) + w1=1/(d1/c1) + err_cnt2Jy_0 = sqrt((w0[0]*d0[0])^2+(w0[1]*d0[1])^2)/(w0[0]+w0[1]) + err_cnt2Jy_1 = sqrt((w1[0]*d1[0])^2+(w1[1]*d1[1])^2)/(w1[0]+w1[1]) + + ; last check: when calibrators are not known, the resulting dummy cnt2Jy values are to be reset to -99.00 + ; as the above offset compensation affects even them (and steers them a little bit from the dummy value) + if (c0[0] lt 0) or (c0[1] lt 0) then begin + cnt2Jy_0 = -99.0 + err_cnt2Jy_0 = -99.0 + endif + if (c1[0] lt 0) or (c1[1] lt 0) then begin + cnt2Jy_1 = -99.0 + err_cnt2Jy_1 = -99.0 + endif + endelse + + + if TPlike eq 2 then begin + printf, Unit0, FORMAT = '(a15,a7,D12.5,d8.1,9(" ",D8.4))', sourcename, hhmmss, data[0].time, freq, $ + el_d, peak_cnt_0, peak_cnt_1,$ + err_cnt_0, err_cnt_1, cnt2Jy_0, cnt2Jy_1, $ + err_cnt2Jy_0, err_cnt2Jy_1;, SNR0, SNR1 + endif else begin + printf, Unit0, FORMAT = '(a15,a7,D12.5,d8.1,D8.4,8(" ",E12.6))', sourcename, hhmmss, data[0].time, freq, $ + el_d, peak_cnt_0, peak_cnt_1,$ + err_cnt_0, err_cnt_1, cnt2Jy_0, cnt2Jy_1, $ + err_cnt2Jy_0, err_cnt2Jy_1;, SNR0, SNR1 + endelse + + ; note that SNR is not reported as it was not recalculated for the combined types + endif + + + if (fitchoice eq 'cubic') or (fitchoice eq 'both') then begin + + if (offsetc[0] eq -99. or offsetc[1] eq -99. or (abs(fw_ratio0c[0]-1.0) ge fw_tol and abs(fw_ratio0c[1]-1.0) ge fw_tol) or (abs(fw_ratio1c[0]-1.0) ge fw_tol and abs(fw_ratio1c[1]-1.0) ge fw_tol)) then begin + peak_cnt_0=-99.0 + peak_cnt_1=-99.0 + err_cnt_0=-99.0 + err_cnt_1=-99.0 + cnt2Jy_0=-99.0 + cnt2Jy_1=-99.0 + err_cnt2Jy_0=-99.0 + err_cnt2Jy_1=-99.0 + endif else begin + ; riscalo le ampiezze e gli errori della polarizzazione 0 per l'offset incrociato + p0c[0]=p0c[0]/exp(-(offsetc[1]*1.66/beamd)^2.) + p0c[1]=p0c[1]/exp(-(offsetc[0]*1.66/beamd)^2.) + e0c[0]=e0c[0]/exp(-(offsetc[1]*1.66/beamd)^2.) + e0c[1]=e0c[1]/exp(-(offsetc[0]*1.66/beamd)^2.) + ; riscalo le ampiezze e gli errori della polarizzazione 1 per l'offset incrociato + p1c[0]=p1c[0]/exp(-(offsetc[1]*1.66/beamd)^2.) + p1c[1]=p1c[1]/exp(-(offsetc[0]*1.66/beamd)^2.) + e1c[0]=e1c[0]/exp(-(offsetc[1]*1.66/beamd)^2.) + e1c[1]=e1c[1]/exp(-(offsetc[0]*1.66/beamd)^2.) + + peak_cnt_0 = cal_wmean(p0c, e0c, /nan) + peak_cnt_1 = cal_wmean(p1c, e1c, /nan) + w0=1/(e0c/p0c) + w1=1/(e1c/p1c) + err_cnt_0 = sqrt((w0[0]*e0[0])^2+(w0[1]*e0[1])^2)/(w0[0]+w0[1]) + err_cnt_1 = sqrt((w1[0]*e1[0])^2+(w1[1]*e1[1])^2)/(w1[0]+w1[1]) + ; err_cnt_0 = stddev(p0, /nan) + ; err_cnt_1 = stddev(p1, /nan) + + ; riscalo i ct2Jy e gli errori della polarizzazione 0 per l'offset incrociato + c0c[0]=c0c[0]*exp(-(offsetc[1]*1.66/beamd)^2.) + c0c[1]=c0c[1]*exp(-(offsetc[0]*1.66/beamd)^2.) + d0c[0]=d0c[0]*exp(-(offsetc[1]*1.66/beamd)^2.) + d0c[1]=d0c[1]*exp(-(offsetc[0]*1.66/beamd)^2.) + ; riscalo i ct2Jy e gli errori della polarizzazione 1 per l'offset incrociato + c1c[0]=c1c[0]*exp(-(offsetc[1]*1.66/beamd)^2.) + c1c[1]=c1c[1]*exp(-(offsetc[0]*1.66/beamd)^2.) + d1c[0]=d1c[0]*exp(-(offsetc[1]*1.66/beamd)^2.) + d1c[1]=d1c[1]*exp(-(offsetc[0]*1.66/beamd)^2.) + + cnt2Jy_0 = cal_wmean(c0c, d0c, /nan) + cnt2Jy_1 = cal_wmean(c1c, d1c, /nan) + w0=1/(d0/c0) + w1=1/(d1/c1) + err_cnt2Jy_0 = sqrt((w0[0]*d0[0])^2+(w0[1]*d0[1])^2)/(w0[0]+w0[1]) + err_cnt2Jy_1 = sqrt((w1[0]*d1[0])^2+(w1[1]*d1[1])^2)/(w1[0]+w1[1]) + + ; last check: when calibrators are not known, the resulting dummy cnt2Jy values are to be reset to -99.00 + ; as the above offset compensation affects even them (and steers them a little bit from the dummy value) + if (c0[0] lt 0) or (c0[1] lt 0) then begin + cnt2Jy_0 = -99.0 + err_cnt2Jy_0 = -99.0 + endif + if (c1[0] lt 0) or (c1[1] lt 0) then begin + cnt2Jy_1 = -99.0 + err_cnt2Jy_1 = -99.0 + endif + + endelse + + if TPlike eq 2 then begin + printf, Unit10, FORMAT = '(a15,a7,D12.5,d8.1,9(" ",D8.4))', sourcename, hhmmss, data[0].time, freq, $ + el_d, peak_cnt_0, peak_cnt_1,$ + err_cnt_0, err_cnt_1, cnt2Jy_0, cnt2Jy_1, $ + err_cnt2Jy_0, err_cnt2Jy_1;, SNR0, SNR1 + endif else begin + printf, Unit10, FORMAT = '(a15,a7,D12.5,d8.1,D8.4,8(" ",E12.6))', sourcename, hhmmss, data[0].time, freq, $ + el_d, peak_cnt_0, peak_cnt_1,$ + err_cnt_0, err_cnt_1, cnt2Jy_0, cnt2Jy_1, $ + err_cnt2Jy_0, err_cnt2Jy_1;, SNR0, SNR1 + endelse + ; note that SNR is not reported as it was not recalculated for the combined types + endif + + endfor + + ; closing all txt output files + close, /ALL + + ; freeing logical units + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + free_lun, Unit0, /force + free_lun, Unit1, /force + endif + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + free_lun, Unit10, /force + free_lun, Unit11, /force + endif + + free_lun, Unit2, /force + + + print, ' ' + print, ' ********************************' + print, ' ******* RUNCALIB is DONE *******' + print, ' *** Next step is CJ2TIMELINE ***' + print, ' ********************************' + print, ' ' + + return +end + diff --git a/runtarget.pro b/runtarget.pro new file mode 100644 index 0000000..d9b7ec4 --- /dev/null +++ b/runtarget.pro @@ -0,0 +1,1242 @@ +pro runtarget, pickpath=pickpath, skypath=skypath, sub=sub, linear=linear, cubic=cubic, plot=plot, skipgc=skipgc + + ; This procedure is devoted to the reduction of cross scans acquired on the observed target sources and + ; apply the previously-measured conversion factors (cnt->Jy). + ; It plots and analyses the data obtained by stacking all the subscans contained in any scan + ; (i.e. subfolder) present inside the working folder. + ; If users want to graphically choose the data folder, they must use: + ; + ; IDL> runtarget, /pickpath in order to select the path to the data parent folders (otherwise the chosen paths are ./CALIBRATORS and ./TARGETS) + ; IDL> runtarget, /skypath in order to select the path to the folder containing skydips (otherwise the chosen path is ./SKYDIPS) + ; + ; Measurements for all the single subscans can be obtained by explicitly choosing this option: + ; + ; IDL> runtarget, /sub + ; + ; The gaussian fitting performed on the (sub)scans can be run using three different baseline-fitting options, + ; in particular: + ; + ; IDL> runtarget, /linear performs, and outputs, only the linear fitting of baselines + ; IDL> runtarget, /cubic performs, and outputs, only the cubic fitting of baselines + ; IDL> runtarget, /both performs, and outputs, both the cubic and linear fitting of baselines + ; + ; By default, the choice is 'both'. + ; + ; At present, graphic output is exceedingly slow, in particular for single subscans (if the /sub option is selected). + ; So it is disabled by default. + ; To enable it, explicitly set the option: + ; + ; IDL> runtarget, /plot + ; + ; + ; The option + ; + ; IDL> runcalib, /skipgc + ; + ; overrides the use of the standard gain curves; it poses Gain=1 for all the elevations, + ; thus practically avoiding the elevation-dependent gain compensation. + ; + ; COMPATIBILITY WITH CONVERTED FILES + ; FITS files achieved with older versions of DISCOS (dating back to years 2008-2015) + ; need to be converted into the proper FITS format (ask for procedure named updatefits.pro). + ; The program automatically handles the TP-like files resulting from the + ; conversion of SARDARA acquisitions, which have raw counts levels in the + ; orders of magnitude of 10E+06..10E+07, producing properly-formatted output tables. + ; + ; Authors: Marcello Giroletti, Simona Righini + ; Last edited: Oct 18, 2017 + ; + + + common logistics, workpath, calpath, sep, skydippath, myextlist, fitchoice, dosingle, doplot + + close, /all ; to close possible pending text files + + sep=path_sep() + + if keyword_set (sub) then dosingle='y' else dosingle='n' + + if keyword_set(pickpath) then begin + workpath=dialog_pickfile(/DIRECTORY, TITLE='Pick the DATA parent folder') + endif else begin + CD, CURRENT=curr + workpath=curr+sep+'TARGETS'+sep + endelse + + if keyword_set(pickpath) then begin + calpath=dialog_pickfile(/DIRECTORY, TITLE='Pick the CALIBRATION FILES parent folder') + endif else begin + calpath=curr+sep+'CALIBRATORS'+sep ; assuming a rigidly-defined folder hierarchy and naming! + endelse + + if keyword_set(skypath) then begin + skydippath=dialog_pickfile(/DIRECTORY, TITLE='Please pick folder containing the Tau.txt file') + endif else begin + CD, CURRENT=curr + skydippath=curr+sep+'SKYDIP'+sep ; assuming a rigidly-defined folder hierarchy and naming! + endelse + + if keyword_set(linear) then begin + fitchoice='linear' + endif else begin + if keyword_set(cubic) then fitchoice='cubic' else fitchoice='both' + endelse + + if keyword_set (plot) then doplot='y' else doplot='n' + + tar_stack, freq=freq + return + +end + + +PRO tar_gfunct, X, A, F, pder + + EX2 = A[0]*EXP(-(X-A[1])^2/(2*A[2]^2)) + F = EX2 + A[3] + A[4]*X +A[5]*X^2 +A[6]*X^3 + IF N_PARAMS() GE 4 THEN $ + pder = [[EX2/A[0]], [EX2*(X-A[1])/(A[2]^2)], [EX2*((X-A[1])^2)/(A[2]^3)], [replicate(1.0, N_ELEMENTS(X))], [X], [X^2], [X^3]] +END + + +function tar_wmean, val, dval, nan = nan, error = error + COMPILE_OPT idl2 + on_error, 2 + + ;- check inputs + if n_params() ne 2 then begin + return, !values.f_nan + endif + sz = n_elements(val) + if n_elements(dval) ne sz then $ + message, 'val and dval must contain the same number of elements' + + + ;- the calculation + ;- requires at least one finite value + if keyword_set(nan) then begin + hit = where(finite(val) and finite(dval), ct) + if ct eq 0 then begin + error = !values.f_nan + return, !values.f_nan + endif + endif + + dw = total(1D / dval^2, nan = keyword_set(nan)) + sum = total(1D * val / dval^2, nan = keyword_set(nan)) / dw + error = sqrt(1D / dw) + return, sum +end + + +pro tar_stack, path=path, out=out, plot=plot, beam=beam, speed=speed, dt=dt, source=source, flux=flux, freq=freq + common logistics + + + ; !EXCEPT=2 ; enabling the indications of location/line numbers in exception messages + !EXCEPT=0 ; silencing the printout of exception messages + + fw_tol=0.3 ; we accept measurements where FWHM is within 30% of the nominal HPBW + + ; reads first file from first scan to set basic info in output files + list = FILE_SEARCH(workpath+'20*', COUNT=number, /TEST_DIRECTORY) + + timings=dblarr(number+1) + gap=dblarr(number+1) + timings[0]=systime(1) + gap[0]=0 + + sublist = FILE_SEARCH(list[0]+sep+'20*.fits', COUNT=subnumber) + ; XXX ripetere la lettura per un file di ogni sub-folder e verificare che siano omogenei + data=MRDFITS(sublist[0],4,/SILENT) + RFinfo=MRDFITS(sublist[0],2,/SILENT) + Sectinfo=MRDFITS(sublist[0],1,/SILENT) + + + ; read frequency and bandwidth from RF table + bw=RFinfo[0].bandWidth + + freq=RFinfo[0].frequency+bw/2.0 + dt=1.0/(Sectinfo[0].sampleRate*1e6) + + ; which telescope was used? Let's get to a 3-char code to be used afterwards + firstmainh=mrdfits(sublist[0],0,firsthead0,/silent) + findh0key, firsthead0, '*ANTENNA =*', keylab, sitename, infolab, siteflag + splitsitename=strsplit(sitename,"'",/extract) + cleansite=strupcase(strcompress(splitsitename[1],/remove_all)) + site=strmid(cleansite,0,3) ; this way, possible cases will be MED, SRT and NOT + + case site of + 'MED': begin + beam=38.7/(freq/1000.0) ; beam in arcmin + ; frequency-dependent normalized gain polynomials + ; http://www.med.ira.inaf.it/ManualeMedicina/English/5.%20Efficiency.htm + ; taken on 2015 September 17 + if (freq gt 4000 and freq lt 6000) then begin + ; this is C band + A_0=[7.9212981E-1,6.2403816E-3,-4.6834953E-5,0] + A_1=[7.9212981E-1,6.2403816E-3,-4.6834953E-5,0] + endif + if (freq gt 8000 and freq lt 10000) then begin + ; this is X band + A_0=[0.730487,0.00773332,-5.54743e-05,0] ; updated 2015, but it is not opacity-corrected + A_1=[0.730487,0.00773332,-5.54743e-05,0] ; updated 2015, but it is not opacity-corrected + ; A_0=[6.1059261E-01,1.0623634E-02,-7.2457279E-05,0] ; origin??? + ; A_1=[6.1059261E-01,1.0623634E-02,-7.2457279E-05,0] ; origin??? + endif + if (freq gt 18000 and freq le 18500) then begin + ; this is K-low band + A_0=[0.8373589,0.005140111,-4.058370E-5,0] ; as of Sept.30, 2015 + A_1=[0.8373589,0.005140111,-4.058370E-5,0] ; as of Sept.30, 2015 + endif + if (freq gt 18500 and freq le 26500) then begin + ; this is K-high band + A_0=[0.7929185,0.005900533,-4.203179E-5,0] ; as of Sept.30, 2015 + A_1=[0.7929185,0.005900533,-4.203179E-5,0] ; as of Sept.30, 2015 + endif + end + 'SRT': begin + ; frequency dependent normalized gain polynomials + ; as provided by Andrea Orlati (commissioning) + if (freq gt 4000 and freq lt 8000) then begin + ; this is C band + beam=18.937/(freq/1000.0) ; beam in arcmin + A_0=[0.939813,0.00139,-0.000009,0] ; XXX To be updated + A_1=[0.939813,0.00139,-0.000009,0] ; XXX To be updated + endif + if (freq ge 18000 and freq lt 27000) then begin + ; this is K band + beam=18.264/(freq/1000.0) ; beam in arcmin + if freq lt 21000 then begin + A_0=[0.508145,0.014839,-0.000111,0] ; from official SRT page + A_1=[0.508145,0.014839,-0.000111,0] ; from official SRT page + endif + if freq ge 21000 and freq le 23000 then begin + A_0=[0.606733,0.01132,-0.000081,0] ; from official SRT page + A_1=[0.606733,0.01132,-0.000081,0] ; from official SRT page + endif + if freq gt 23000 then begin + A_0=[0.533025,0.015561,-0.000129,0] ; from official SRT page + A_1=[0.533025,0.015561,-0.000129,0] ; from official SRT page + endif + endif + end + 'NOT': begin + ;beam=38.7/(freq/1000.0) ; beam in arcmin - CLONED FROM MEDICINA XXX + if (freq gt 4000 and freq lt 6000) then begin + ; this is C band + beam=7.5 ; fixed value + ; A_0=[0.98463659,0.00017053038,1.9348601e-09,0] ; old + ; A_1=[0.98463659,0.00017053038,1.9348601e-09,0] ; old + A_0=[0.95227062926,0.00198726460815,-2.0685473288e-05,0] ; Cassaro, Oct 16th 2017 + A_1=[0.95227062926,0.00198726460815,-2.0685473288e-05,0] ; Cassaro, Oct 16th 2017 + endif + if (freq gt 8000 and freq lt 10000) then begin + ; this is X band + print, 'X Curve for Noto is not available, yet. Returning.' + return + endif + if (freq gt 18000 and freq le 26500) then begin + ; this is K band + print, 'K Curve for Noto is not available, yet. Returning.' + return + endif + end + else: begin + print, 'UNKNOWN SITE: cannot apply gain curves. Returning.' + return + end + endcase + beamd=beam/60. + sd=beamd/(2*SQRT(2*ALOG(2))) + + + + ; open the main output files + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + openw, Unit0, workpath+'tar_stack_lin_final.txt', /GET_LUN + openw, Unit1, workpath+'tar_stack_lin_sep.txt', /GET_LUN + printf, Unit0, FORMAT = '(A7,1X,A3)', "Site =",site + endif + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + openw, Unit10, workpath+'tar_stack_cub_final.txt', /GET_LUN + openw, Unit11, workpath+'tar_stack_cub_sep.txt', /GET_LUN + endif + + if strmatch(sublist[0],'*TPlike*.*',/FOLD_CASE) then begin + TPlike=1 ; this is a SARDARA file converted/integrated in TotalPower-like format + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + printf, Unit0, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 flux_0 flux_1 e_flux_0 e_flux_1' + printf, Unit0, ' MHz deg count count count count Jy Jy Jy Jy' + printf, Unit1, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 flux_0 flux_1 e_flux_0 e_flux_1 SNR_0 SNR_1 Offset FW/HP Type' + printf, Unit1, ' MHz deg count count count count Jy Jy Jy Jy deg ' + endif + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + printf, Unit10, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 flux_0 flux_1 e_flux_0 e_flux_1' + printf, Unit10, ' MHz deg count count count count Jy Jy Jy Jy' + printf, Unit11, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 flux_0 flux_1 e_flux_0 e_flux_1 SNR_0 SNR_1 Offset FW/HP Type' + printf, Unit11, ' MHz deg count count count count Jy Jy Jy Jy deg ' + endif + endif else begin + TPlike=2 ; this is an original TotalPower FITS + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + printf, Unit0, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 flux_0 flux_1 e_flux_0 e_flux_1' + printf, Unit0, ' MHz deg count count count count Jy Jy Jy Jy' + printf, Unit1, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 flux_0 flux_1 e_flux_0 e_flux_1 SNR_0 SNR_1 Offset FW/HP Type' + printf, Unit1, ' MHz deg count count count count Jy Jy Jy Jy deg ' + endif + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + printf, Unit10, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 flux_0 flux_1 e_flux_0 e_flux_1' + printf, Unit10, ' MHz deg count count count count Jy Jy Jy Jy' + printf, Unit11, ' Name HHMMSS # MJD Freq Elev Amp_0 Amp_1 e_Amp_0 e_Amp_1 flux_0 flux_1 e_flux_0 e_flux_1 SNR_0 SNR_1 Offset FW/HP Type' + printf, Unit11, ' MHz deg count count count count Jy Jy Jy Jy deg ' + endif + endelse + + openw, Unit2, workpath+'tar_prints.txt', /GET_LUN + + + ; reads skydip information + ans=file_test(skydippath+'Tau.txt') + if ans eq 1 then begin + readcol, skydippath+'Tau.txt', tau_mjd, tau_freq, t30l, t90l, t30r, t90r, tau_empL, tau_empR, tau0L, tau0R, format='(d,d,d,d,d,d,d,d,d,d)', /SILENT + for t=0,n_elements(tau_mjd)-1 do begin + ; if for any reason tau_emp is considered more reliable than tau0, uncomment the following two lines + ; tau0L[t]=tau_empL[t] + ; tau0R[t]=tau_empR[t] + endfor + caltau = 1 + print, '++++++++++++++++++++++++++++++++' + print, 'Reading tau0 values from Tau.txt + print, '++++++++++++++++++++++++++++++++' + + endif else begin + print, ' ' + print, '+++++++++++++++++++++++++++++++++++++++++++++' + print, 'Beware: no Tau.txt available, assuming tau0=0 + print, '+++++++++++++++++++++++++++++++++++++++++++++' + wait, 0.5 + tau0L=[0,0] + tau0R=[0,0] + caltau = 0 + endelse + + ; definition of arrays + p0=fltArr(2) + p1=fltArr(2) + e0=fltArr(2) + e1=fltArr(2) + c0=fltArr(2) ; count to jansky array + c1=fltArr(2) + d0=fltArr(2) ; error on count to jansky array + d1=fltArr(2) + SNR0=fltArr(2) + SNR1=fltArr(2) + offset=dblArr(2) + p0c=fltArr(2) + p1c=fltArr(2) + e0c=fltArr(2) + e1c=fltArr(2) + c0c=fltArr(2) ; count to jansky array + c1c=fltArr(2) + d0c=fltArr(2) ; error on count to jansky array + d1c=fltArr(2) + SNR0c=fltArr(2) + SNR1c=fltArr(2) + offsetc=dblArr(2) + fw_ratio0=fltarr(2) + fw_ratio1=fltarr(2) + fw_ratio0c=fltarr(2) + fw_ratio1c=fltarr(2) + + + ; scan-based analysis + for i =0,number-1 do begin + print, ' ' + print, '=========================================' + print, 'Now analysing scan ', list[i] + print, '=========================================' + sommaL=0 + sommaR=0 + sommaLA = 0 + sommaRA = 0 + numLA = 0 + numRA = 0 + iiA = 0 + ffA = 0 + sommaLB = 0 + sommaRB = 0 + numLB = 0 + numRB = 0 + iiB = 0 + ffB = 0 + + + splitfoldname1=strsplit(list[i],sep,/extract) + foldname=splitfoldname1[-1] + splitfoldname2=strsplit(foldname,'-',/extract) + hhmmss = splitfoldname2[1] + ; splitfoldname=strsplit(list[i],'-',/extract) + ; hhmmss = splitfoldname[1] + sublist = FILE_SEARCH(list[i]+sep+'2*.fits', COUNT=subnumber) + scanname=file_basename(list[i]) + + ; reads flagging info + err=0 + checkfile = list[i]+sep+'checkfile_'+scanname+'.txt' + openr, Unit9, checkfile, ERROR=err, /GET_LUN + subscan=file_basename(sublist) + flagcode=strarr(subnumber) + flagcode[*]='+1cx' + if (err lt 0) then begin + print, "No valid flagging information for this scan" + print, "Temporarily assuming all scans are good" + endif else begin + readcol, checkfile, F='a,a,a', subscan_t, flagcode_t, who, /SILENT + if (n_elements(subscan_t) lt subnumber) then begin + print, "Some subscans miss flagging information" + print, "Temporarily assuming those subscans are good" + endif + for n_flag = 0, n_elements(subscan_t)-1 do begin + for j_flag = 0, subnumber-1 do begin + if (subscan_t[n_flag] eq subscan[j_flag]) then begin + flagcode[j_flag]=flagcode_t[n_flag] + endif + endfor + endfor + close, Unit9 + free_lun, Unit9 + endelse + + + if dosingle eq 'y' then begin + ; creating output files for single-subscan analysis + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + openw, Unit3, list[i]+sep+scanname+'_single_linear.txt', /GET_LUN + printf, Unit3,' Name # MJD El(r) El(d) cnt_0 cnt_1 e_c_0 e_c_1 offL(d) offR(d) FW/HP scantype' + endif + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + openw, Unit4, list[i]+sep+scanname+'_single_cubic.txt', /GET_LUN + printf, Unit4,' Name # MJD El(r) El(d) cnt_0 cnt_1 e_c_0 e_c_1 offL(d) offR(d) FW/HP scantype' + endif + openw, Unit5, list[i]+sep+'tar_prints.txt', /GET_LUN + endif + + ; Dynamic arrays (for amplitudes and errors, for each section) + + Lampl=[ ] + Rampl=[ ] + Lerrs=[ ] + Rerrs=[ ] + offset_sub=[ ] + type=[ ] + mean_lat=[ ] + + ; readout of the previously produced files, containing the list of counts-to-Jy factors achieved on calibrators + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + lincal=file_search(calpath+sep+'cnt2jy_exlin.txt', count=lincalnum) + if lincalnum eq 0 then begin + print, '++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++' + print, 'cnt2Jy_lin factors not found. Assuming cnt2Jy=1 and returning all amplitudes in raw counts + print, '++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++' + endif else begin + readcol, calpath+sep+'cnt2jy_exlin.txt', lincalname, linfre, lintime_i, lintime_f, linm_0, linc2j_0, err_linm_0, err_linc2j_0, linm_1, linc2j_1, err_linm_1, err_linc2j_1, format='(A,D,D,D,D,D,D,D,D,D,D,D)', /SILENT + endelse + endif + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + cubcal=file_search(calpath+sep+'cnt2jy_excub.txt', count=cubcalnum) + if cubcalnum eq 0 then begin + print, '++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++' + print, 'cnt2Jy_cub factors not found. Assuming cnt2Jy=1 and returning all amplitudes in raw counts + print, '++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++' + endif else begin + readcol, calpath+sep+'cnt2jy_excub.txt', cubcalname, cubfre, cubtime_i, cubtime_f, cubm_0, cubc2j_0, err_cubm_0, err_cubc2j_0, cubm_1, cubc2j_1, err_cubm_1, err_cubc2j_1, format='(A,D,D,D,D,D,D,D,D,D,D,D)', /SILENT + endelse + endif + + + ; subscan-based measurements + for j=0,subnumber-1 do begin + ; READ DATA TABLE (binary data tables) + data=MRDFITS(sublist[j],4,/SILENT) + RFinfo=MRDFITS(sublist[j],2,/SILENT) + Sectinfo=MRDFITS(sublist[j],1,/SILENT) + ; read frequency and bandwith from RF table + bw=RFinfo[0].bandWidth + freq=RFinfo[0].frequency+bw/2.0 + dt=1.0/(Sectinfo[0].sampleRate*1e6) + ndat = (size(data.time)) + + ;assigning an MJD value to this subscan (first sample) + mymjd=data[0].time + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + if lincalnum eq 0 then begin + ; forcing cnt2Jy values to 1.0 --> final "fluxes" will be expressed in raw counts + mylinc2j_0=1.0 + mylinc2j_1=1.0 + err_mylinc2j_0=1.0 + err_mylinc2j_1=1.0 + endif else begin + ; find which cnt2Jy interval contains such MJD - LINEAR-FIT-BASED data + my_lin_interval=where(((mymjd-lintime_i) ge 0) and ((mymjd-lintime_f) le 0)) + if my_lin_interval eq -1 then begin + ; the target observation time does not belong to any of the available intervals + ; thus some degree of approximation is needed + if mymjd le lintime_i[0] then begin + ; my target was observed before the first available cnt2Jy interval: I'll use the first cnt2Jy value + my_lin_interval=0 + endif else begin + if mymjd ge lintime_f[-1] then begin + ; my target was observed after the last available cnt2Jy interval: I'll use the last cnt2Jy value + my_lin_interval=-1 + endif else begin + ; my target was observed in between calibrators, but outside any of the intervals available: + ; going to the nearest value/interval as a compromise! + min_dist_to_init=min(abs(mymjd-lintime_i)) + min_dist_to_final=min(abs(mymjd-lintime_f)) + if min_dist_to_final lt min_dist_to_init then begin + my_lin_interval=where(abs(mymjd-lintime_f) eq min(abs(mymjd-lintime_f))) + endif else begin + my_lin_interval=where(abs(mymjd-lintime_i) eq min(abs(mymjd-lintime_i))) + endelse + endelse + endelse + endif + ; assigning the identified best-matching cnt2Jy values to the final attributes, for the subsequent analysis + mylinc2j_0=linc2j_0[my_lin_interval]+mymjd*linm_0[my_lin_interval] + mylinc2j_1=linc2j_1[my_lin_interval]+mymjd*linm_1[my_lin_interval] + err_mylinc2j_0=mylinc2j_0*sqrt(err_linc2j_0[my_lin_interval]^2+mymjd*err_linm_0[my_lin_interval]^2) + err_mylinc2j_1=mylinc2j_1*sqrt(err_linc2j_1[my_lin_interval]^2+mymjd*err_linm_1[my_lin_interval]^2) + endelse + endif + + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + if cubcalnum eq 0 then begin + ; forcing cnt2Jy values to 1.0 --> final "fluxes" will be expressed in raw counts + mycubc2j_0=1.0 + mycubc2j_1=1.0 + err_mycubc2j_0=1.0 + err_mycubc2j_1=1.0 + endif else begin + ; find which cnt2Jy interval contains such MJD - CUBIC-FIT-BASED data + my_cub_interval=where(((mymjd-cubtime_i) ge 0) and ((mymjd-cubtime_f) le 0)) + if my_cub_interval eq -1 then begin + ; the target observation time does not belong to any of the available intervals + ; thus some degree of approximation is needed + if mymjd le cubtime_i[0] then begin + ; my target was observed before the first available cnt2Jy interval: I'll use the first cnt2Jy value + my_cub_interval=0 + endif else begin + if mymjd ge cubtime_f[-1] then begin + ; my target was observed after the last available cnt2Jy interval: I'll use the last cnt2Jy value + my_cub_interval=-1 + endif else begin + ; my target was observed in between calibrators, but outside any of the intervals available: + ; going to the nearest value/interval as a compromise! + min_dist_to_init=min(abs(mymjd-cubtime_i)) + min_dist_to_final=min(abs(mymjd-cubtime_f)) + if min_dist_to_final lt min_dist_to_init then begin + my_cub_interval=where(abs(mymjd-cubtime_f) eq min(abs(mymjd-cubtime_f))) + endif else begin + my_cub_interval=where(abs(mymjd-cubtime_i) eq min(abs(mymjd-cubtime_i))) + endelse + endelse + endelse + endif + ; assigning the identified best-matching cnt2Jy values to the final attributes, for the subsequent analysis + mycubc2j_0=cubc2j_0[my_cub_interval]+mymjd*cubm_0[my_cub_interval] + mycubc2j_1=cubc2j_1[my_cub_interval]+mymjd*cubm_1[my_cub_interval] + err_mycubc2j_0=mycubc2j_0*sqrt(err_cubc2j_0[my_cub_interval]^2+mymjd*err_cubm_0[my_cub_interval]^2) + err_mycubc2j_1=mycubc2j_1*sqrt(err_cubc2j_1[my_cub_interval]^2+mymjd*err_cubm_1[my_cub_interval]^2) + endelse + endif + + ; UNCOMMENT IF WANTING TO PRINT THE cnt2Jy choice at runtime + ; print, '***********************************************' + ; if fitchoice eq 'linear' or fitchoice eq 'both' then begin + ; print, '** Chosen cnt2Jy values for LIN (pos, v0, v1): ', my_lin_interval+1, mylinc2j_0, mylinc2j_1 + ; endif + ; if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + ; print, '** Chosen cnt2Jy values for CUB (pos, v0, v1): ', my_cub_interval+1, mycubc2j_0, mycubc2j_1 + ; endif + ; print, '***********************************************' + + ; takes source name and subscantype from header + mainh=mrdfits(sublist[j],0,head0,/silent) + + ; assessing whether this is an old converted FITS + findh0key, head0, '*Converted FITS*', keylab, converted, infolab, convertedflag + findh0key, head0, '*Obtained from ESCS 0.1 TP FITS*', keylab, version, infolab, oldestflag + + if convertedflag eq 0 then begin ; this IS NOT a converted file: ESCS 0.3 keyword names hold + findh0key, head0, '*SOURCE*', keylab, target, infolab, tarflag + findh0key, head0, '*SubScanType =*', keylab, orscantype, infolab, tarflag + findh0key, head0, '*Azimuth Offset =*', keylab, hAZoff, infolab, typeflag ; offset utente e di sistema + findh0key, head0, '*Elevation Offset =*', keylab, hELoff, infolab, typeflag ; offset utente e di sistema + findh0key, head0, '*RightAscension Offset =*', keylab, hRAoff, infolab, typeflag ; offset utente e di sistema + findh0key, head0, '*Declination Offset =*', keylab, hDECoff, infolab, typeflag ; offset utente e di sistema + findh0key, head0, '*RightAscension =*', keylab, ras, infolab, typeflag ; RA nominale del target, radianti + findh0key, head0, '*Declination =*', keylab, decs, infolab, typeflag ; DEC nominale del target, radianti + endif else begin ; this IS a converted file: custom 8-char keyword names hold according to original FITS date, and offsets are not present + findh0key, head0, '*SOURCE*', keylab, target, infolab, tarflag + findh0key, head0, '*SUBSCANT=*', keylab, orscantype, infolab, tarflag + hAZoff=0.0d + hELoff=0.0d + hRAoff=0.0d + hDECoff=0.0d + if oldestflag eq 0 then begin + findh0key, head0, '*RightAscension =*', keylab, ras, infolab, typeflag ; RA nominale del target, radianti + findh0key, head0, '*Declination =*', keylab, decs, infolab, typeflag ; DEC nominale del target, radianti + endif else begin ; these are the oldest FITS, where several more keywords have converted names + findh0key, head0, '*RIGHTASC=*', keylab, ras, infolab, typeflag ; RA nominale del target, radianti + findh0key, head0, '*DECLINAT=*', keylab, decs, infolab, typeflag ; DEC nominale del target, radianti + endelse + endelse + + + target=strcompress(target, /remove_all) + cleantarget=strsplit(target,"'",/extract) + target=cleantarget[0] + sourcename=target + + ; nominal target coordinates, converted in degrees + rasd=ras*180.0d/!dpi + decsd=decs*180.0d/!dpi + + elapsed=(double(data.time)-double(data[0].time))*24.0*3600.0 ; elapsed time from first sample + duration=dt*ndat[1]/60.0 ; duration of the whole subscan (approx) + ; RA and DEC span of the subscan + deltaDec=abs(data[0].decj2000-data[-1].decj2000)/!dpi*180.0 + deltaRA=abs(data[0].raj2000-data[-1].raj2000)*cos(mean(data.decj2000))/!dpi*180.0 + + ; determine speed (deg/min) along scan direction + ; if type is AZ or EL, treats is as DEC or RA depending on which quantity spans a larger range + splitscantype=strsplit(orscantype,"'",/extract) + cleanscantype=strcompress(splitscantype[1], /remove_all) ; cleaning the string so that it become more manageable + case cleanscantype of + 'DEC': begin + if (data[0].decj2000 lt data[-1].decj2000) then scantype = 1 else scantype = 2 + speed=deltaDec/duration + end + 'RA': begin + if (data[0].raj2000 lt data[-1].raj2000) then scantype = 3 else scantype = 4 + speed=deltaRA/duration + end + 'AZ': begin + if (deltaDec gt deltaRA) then begin + if (data[0].decj2000 lt data[-1].decj2000) then scantype = 1 else scantype = 2 + speed=deltaDec/duration + endif + if (deltaDec lt deltaRA) then begin + if (data[0].raj2000 lt data[-1].raj2000) then scantype = 3 else scantype = 4 + speed=deltaRA/duration + endif + end + 'EL': begin + if (deltaDec gt deltaRA) then begin + if (data[0].decj2000 lt data[-1].decj2000) then scantype = 1 else scantype = 2 + speed=deltaDec/duration + endif + if (deltaDec lt deltaRA) then begin + if (data[0].raj2000 lt data[-1].raj2000) then scantype = 3 else scantype = 4 + speed=deltaRA/duration + endif + end + 'UNKNOWN': begin + if (deltaDec gt deltaRA) then begin + if (data[0].decj2000 lt data[-1].decj2000) then scantype = 1 else scantype = 2 + speed=deltaDec/duration + endif + if (deltaDec lt deltaRA) then begin + if (data[0].raj2000 lt data[-1].raj2000) then scantype = 3 else scantype = 4 + speed=deltaRA/duration + endif + end + endcase + + ; check for spikes in X-band wide bandwidth data XXX Perché farlo solo per la banda X? + if (bw gt 500 and freq lt 10000) then begin + for k=3, ndat[1]-4 do begin + temp0=[data[k-3].ch0,data[k-2].ch0,data[k-1].ch0,data[k+1].ch0,data[k+2].ch0,data[k+3].ch0] + temp1=[data[k-3].ch1,data[k-2].ch1,data[k-1].ch1,data[k+1].ch1,data[k+2].ch1,data[k+3].ch1] + avech0=mean(temp0) + avech1=mean(temp1) + stdch0=stddev(temp0) + stdch1=stddev(temp1) + if ((data[k].ch0-avech0) gt 7*stdch0) then begin + data[k].ch0=avech0 + endif + if ((data[k].ch1-avech1) gt 7*stdch1) then begin + data[k].ch1=avech1 + endif + endfor + endif + + ; reverse arrays for scans in decreasing RA or Dec + if (scantype eq 2 or scantype eq 4) then begin + ; data=reverse(data) + data.raj2000=reverse(data.raj2000) + data.decj2000=reverse(data.decj2000) + data.ch0=reverse(data.ch0) + data.ch1=reverse(data.ch1) + endif + + ; computing the sample offset with respect to the source position + xdec=(data.decj2000-decs)/!dpi*180.0 + xra=(data.raj2000-ras)*cos(mean(data.decj2000))/!dpi*180.0 + + scanname = sublist[j] + case flagcode[j] of + '-1c0': begin + k0=1 + k1=0 + end + '-1c1': begin + k0=0 + k1=1 + end + '+0cx': begin + k0=1 + k1=1 + end + '+1cx': begin + k0=1 + k1=1 + end + '-1cx': begin + k0=0 + k1=0 + end + endcase + + ; shift scans of type 'Dec' so they align with first one; sets start and end value of array with data from every scan + if (scantype eq 1 or scantype eq 2) then begin + if ((numLA+numRA) eq 0) then begin + dataA=data + ndatA=ndat + xdecA=xdec + iiA=0 + if (k0+k1 gt 0) then ffA=ndatA[1]-1 + decoffA = min(xdecA,xdec0A,/ABSOLUTE) ; xdec0A is the number of sample at the source declination for scan A + endif else begin + decoff = min(xdec,xdec0,/ABSOLUTE) + ; shift data by (xdec0-xdec0B) samples + ndelta=xdec0-xdec0A + if (ndelta gt 0) then begin + iiA=iiA + if (k0+k1 gt 0) then ffA=min([ffA,ndat[1]-1-ndelta]) + for k=iiA,ffA do begin + if (k0+k1 gt 0) then data[k]=data[k+ndelta] + endfor + endif + if (ndelta lt 0) then begin + if (k0+k1 gt 0) then iiA=max([iiA,0-ndelta]) + ; here we have a problem: we would like to start from ffA and go down to iiA, since ffA can be filled by shifting "data" right. but data only has dimension ndat, not ndat-ndelta (ndelta is <0) + ; ffA=min([ffA,ndat[1]-1-ndelta]) + if (k0+k1 gt 0) then ffA=min([ffA,ndat[1]-1]) + for k=ffA,iiA,-1 do begin + if (k0+k1 gt 0) then data[k]=data[k+ndelta] + endfor + endif + if (ndelta eq 0) then begin + if (k0+k1 gt 0) then ffA=min([ffA,ndat[1]-1]) + endif + endelse + endif + + ; shifts scans of type 'RA' so they align with first one; sets start and end value of array with data from every scan + if (scantype eq 3 or scantype eq 4) then begin + if ((numLB+numRB) eq 0) then begin + dataB=data + ndatB=ndat + xraB=xra + iiB=0 + if (k0+k1 gt 0) then ffB=ndatB[1]-1 + raoffB = min(xraB,xra0B,/ABSOLUTE) ; xra0B is the number of sample at the source right ascension for scan B + endif else begin + raoff = min(xra,xra0,/ABSOLUTE) + ; shift data by (xra0-xra0B) samples + ndelta=xra0-xra0B + if (ndelta gt 0) then begin + iiB=iiB + if (k0+k1 gt 0) then ffB=min([ffB,ndat[1]-1-ndelta]) + for k=iiB,ffB do begin + if (k0+k1 gt 0) then data[k]=data[k+ndelta] + endfor + endif + if (ndelta lt 0) then begin + if (k0+k1 gt 0) then iiB=max([iiB,0-ndelta]) + ; here we have a problem: we would like to start from ffB and go down to iiB, since ffB can be filled by shifting "data" right. but data only has dimension ndat, not ndat-ndelta (ndelta is <0) + ; ffB=min([ffB,ndat[1]-1-ndelta]) + if (k0+k1 gt 0) then ffB=min([ffB,ndat[1]-1]) + for k=ffB,iiB,-1 do begin + if (k0+k1 gt 0) then data[k]=data[k+ndelta] + endfor + endif + if (ndelta eq 0) then begin + if (k0+k1 gt 0) then ffB=min([ffB,ndat[1]-1]) + endif + endelse + endif + + time_s = (data.time-data[round(ndat[1]/2)].time)*24.d0*3600.d0 + + midscan=FIX(ndat[1]/2) + beamspan= beam/speed ; duration of acquisition for one beamwidth + Nsamples=ROUND(beamspan/dt) ; number of samples in one beamwidth + + el_d=data[midscan].el*180.0/!dpi + + g_0=A_0[0]+A_0[1]*el_d+A_0[2]*el_d^2+A_0[3]*el_d^3 + g_1=A_1[0]+A_1[1]*el_d+A_1[2]*el_d^2+A_1[3]*el_d^3 + + if keyword_set(skipgc) then begin + ; Overriding the gain computation and applying a flat gain curve (constant value = 1.0) + g_0=1.0 + g_1=1.0 + endif + + ; actually sums counts only for scans flagged as good and sets inputs for single scan fit + if ( scantype eq 1 or scantype eq 2) then begin + if (k0 ne 0) then sommaLA = sommaLA + k0*data.ch0/g_0 + if (k1 ne 0) then sommaRA = sommaRA + k1*data.ch1/g_1 + numLA = numLA + k0*1 + numRA = numRA + k1*1 + ascissa=(data.decj2000-hDECoff)*180d/!dpi + gpos=decsd + endif else begin + if (k0 ne 0) then sommaLB = sommaLB + k0*data.ch0/g_0 + if (k1 ne 0) then sommaRB = sommaRB + k1*data.ch1/g_1 + numLB = numLB + k0*1 + numRB = numRB + k1*1 + ascissa=(data.raj2000-hRAoff/cos(mean(data.decj2000)))*180d/!dpi + gpos=rasd + endelse + + + + if dosingle eq 'y' then begin + + if (fitchoice eq 'linear') or (fitchoice eq 'both') then begin + ; fit single subscans with linear + gaussian + Lfwhm=0 + Rfwhm=0 + yy0=data.ch0/g_0 + yy1=data.ch1/g_1 + if (scantype eq 1 or scantype eq 2) then tipo=0 else tipo=1 + + ; pro targetfit, scanflag,stacflag,polyflag,section,tipo,namefile,Out3,cnt2Jy,err_cnt2Jy,datael,tau0, xx,yy,ii,ff,x_mid,Nsamples,sd_sub,gpos,decsd,rasd, fwhm, n_off, off, p, e, c, d, SNR, plo, doplot + targetfit, k0,'single','linear','Ch_0',tipo,subscan[j],Unit5,mylinc2j_0,err_mylinc2j_0,data[0].el,tau0L,ascissa,yy0,0,ndat[1]-1,midscan,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Lfwhm, nL, offL, peak_cnt_0, err_cnt_0, dum3, dum4, dum5, psingle, doplot + targetfit, k1,'single','linear','Ch_1',tipo,subscan[j],Unit5,mylinc2j_1,err_mylinc2j_1,data[0].el,tau0R,ascissa,yy1,0,ndat[1]-1,midscan,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Rfwhm, nR, offR, peak_cnt_1, err_cnt_1, dum3, dum4, dum5, psingle, doplot + ; NOTA: modificata la gestione del caso /sub; ora il salvataggio del PDF avviene direttamente dentro alla procedura dataplot, + ; che salva un file per ogni chiamata. + + + if TPlike eq 2 then begin + printf, Unit3, FORMAT = '(a15,D12.5,9(" ",D08.4),i)', sourcename, data[midscan].time,$ + data[midscan].el, el_d, peak_cnt_0, peak_cnt_1, err_cnt_0, err_cnt_1, offL, offR, Lfwhm/beam, scantype + endif else begin + printf, Unit3, FORMAT = '(a15,1X,D12.5,1X,D08.4,1X,D08.4,1X,E12.6,1X,E12.6,1X,E12.6,1X,E12.6,1X,D08.4,1X,D08.4,1X,D08.4,1X,i)', sourcename, data[midscan].time,$ + data[midscan].el, el_d, peak_cnt_0, peak_cnt_1, err_cnt_0, err_cnt_1, offL, offR, Lfwhm/beam, scantype + endelse + endif + + if (fitchoice eq 'cubic') or (fitchoice eq 'both') then begin + ; fit single subscans with cubic + gaussian + Lfwhm=0 + Rfwhm=0 + yy0=data.ch0/g_0 + yy1=data.ch1/g_1 + if (scantype eq 1 or scantype eq 2) then tipo=0 else tipo=1 + targetfit, k0,'single','cubic','Ch_0',tipo,list[i],subscan[j],Unit5,mycubc2j_0,err_mycubc2j_0,data[0].el,tau0L,ascissa,yy0,0,ndat[1]-1,midscan,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Lfwhm, nL, offL, peak_cnt_0, err_cnt_0, dum3, dum4, dum5, psingle, doplot + targetfit, k1,'single','cubic','Ch_1',tipo,list[i],subscan[j],Unit5,mycubc2j_1,err_mycubc2j_1,data[0].el,tau0R,ascissa,yy1,0,ndat[1]-1,midscan,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Rfwhm, nR, offR, peak_cnt_1, err_cnt_1, dum3, dum4, dum5, psingle, doplot + ; NOTA: modificata la gestione del caso /sub; ora il salvataggio del PDF avviene direttamente dentro alla procedura dataplot, + ; che salva un file per ogni chiamata. + + + if TPlike eq 2 then begin + printf, Unit4, FORMAT = '(a15,D12.5,9(" ",D08.4),i)', sourcename, data[midscan].time,$ + data[midscan].el, el_d, peak_cnt_0, peak_cnt_1, err_cnt_0, err_cnt_1, offL, offR, Lfwhm/beam, scantype + endif else begin + printf, Unit4, FORMAT = '(a15,1X,D12.5,1X,D08.4,1X,D08.4,1X,E12.6,1X,E12.6,1X,E12.6,1X,E12.6,1X,D08.4,1X,D08.4,1X,D08.4,1X,i)', sourcename, data[midscan].time,$ + data[midscan].el, el_d, peak_cnt_0, peak_cnt_1, err_cnt_0, err_cnt_1, offL, offR, Lfwhm/beam, scantype + endelse + endif + + if (nL+nR gt 0) then offset_sub = [offset_sub,(offL+offR)/(nL+nR)] + + endif + + endfor + + if dosingle eq 'y' then begin + ; closing the output files and freeing memory + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + close, Unit3 + FREE_LUN, Unit3 + endif + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + close, Unit4 + FREE_LUN, Unit4 + endif + close, Unit5 + free_lun, Unit5, /force + endif + + ; stacked analysis begins + + for tipo = 0,1 do begin + ascissa=time_s ; default if all subscans are void + + ; set some inputs + case tipo of + 0: begin + if (numLA+numRA gt 0) then begin + if (numLA gt 0) then sommaL = sommaLA/numLA else sommaL = sommaLA + if (numRA gt 0) then sommaR = sommaRA/numRA else sommaR = sommaRA + numL = numLA + numR = numRA + ii=iiA + ff=ffA + ; tipo=0 means scan declination, ie scantype Dec, ie constant RA + ascissa=(dataA.decj2000-hDECoff)*180d/!dpi + xdiffA = min(ascissa-decsd,x_mid,/ABSOLUTE) ; x_mid contains the sample number nearest to the source nominal position + gpos=decsd + endif else begin + numL=0 + numR=0 + endelse + end + 1: begin + if (numLB+numRB gt 0) then begin + if (numLB gt 0) then sommaL = sommaLB/numLB else sommaL = sommaLB + if (numRB gt 0) then sommaR = sommaRB/numRB else sommaR = sommaRB + numL = numLB + numR = numRB + ii=iiA + ff=ffB + ; tipo=1 means scan right asc., ie scantype RA, ie constant Dec + ascissa=(dataB.raj2000-hRAoff/cos(mean(dataB.decj2000)))*180d/!dpi + xdiffB = min(ascissa-rasd,x_mid,/ABSOLUTE) ; x_mid contains the sample number nearest to the source nominal position + gpos=rasd + endif else begin + numL=0 + numR=0 + endelse + end + endcase + + PRINTF, Unit2, " " + PRINTF, Unit2, "****** "+scanname+" ******" + PRINTF, Unit2, FORMAT = '("* MJD = ", D014.8)', data[midscan].time + PRINTF, Unit2, "* nsamples = ", ndat[1] + PRINTF, Unit2, "* Beamspan = ", beamspan, " [s]" + PRINTF, Unit2, "* Beam samples = ", Nsamples + + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + + ; LINEAR FIT OF STACKED SCAN + Lfwhm=0 + Rfwhm=0 + + targetfit, numL,'stacked','linear','Ch_0',tipo,list[i],list[i],Unit2,mylinc2j_0,err_mylinc2j_0,data[0].el,tau0L,ascissa,sommaL,ii,ff,x_mid,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Lfwhm, n_offL, offL, dum1, dum2, dum3, dum4, dum5, pstack, doplot + p0[tipo]=dum1 + e0[tipo]=dum2 + c0[tipo]=dum3 + d0[tipo]=dum4 + SNR0[tipo]=dum5 + fw_ratio0[tipo]=abs(Lfwhm/beam) + + targetfit, numR,'stacked','linear','Ch_1',tipo,list[i],list[i],Unit2,mylinc2j_1,err_mylinc2j_1,data[0].el,tau0R,ascissa,sommaR,ii,ff,x_mid,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Rfwhm, n_offR, offR, dum1, dum2, dum3, dum4, dum5, pstack, doplot + p1[tipo]=dum1 + e1[tipo]=dum2 + c1[tipo]=dum3 + d1[tipo]=dum4 + SNR1[tipo]=dum5 + fw_ratio1[tipo]=abs(Rfwhm/beam) + + if doplot eq 'y' then begin + if (tipo eq 1 and fitchoice eq 'linear') then begin + pstack.save, list[i]+'_stacked.pdf', /LANDSCAPE, RESOLUTION=300, /APPEND, /CLOSE + endif else begin + pstack.save, list[i]+'_stacked.pdf', /LANDSCAPE, RESOLUTION=300, /APPEND + endelse + endif + + if (n_offL+n_offR gt 0) then offset[tipo]=(offL+offR)/(n_offL+n_offR) else offset[tipo] = -99. + + el_d=data[0].el*180.0/!dpi + + if tipo eq 0 then sclab='DEC' + if tipo eq 1 then sclab='RA' + + if TPlike eq 2 then begin + printf, Unit1, FORMAT = '(a15,a7,D12.5,d8.1,13(" ",D8.4),A7)', target, hhmmss, data[0].time, freq, $ + el_d, p0[tipo],p1[tipo],e0[tipo],e1[tipo],c0[tipo],c1[tipo],d0[tipo],d1[tipo],SNR0[tipo], SNR1[tipo], offset[tipo], fw_ratio0[tipo], sclab + endif else begin + if c0[tipo] lt 0 or c1[tipo] lt 0 then begin ; here there are -99.000 values + printf, Unit1, FORMAT = '(a15,a7,D12.5,d8.1,d8.4,4(" ",D12.6),8(" ",D8.4),A7)', target, hhmmss, data[0].time, freq, $ + el_d, p0[tipo],p1[tipo],e0[tipo],e1[tipo],c0[tipo],c1[tipo],d0[tipo],d1[tipo],SNR0[tipo], SNR1[tipo], offset[tipo], fw_ratio0[tipo], sclab + endif else begin + printf, Unit1, FORMAT = '(a15,a7,D12.5,d8.1,d8.4,4(" ",E12.6),8(" ",D8.4),A7)', target, hhmmss, data[0].time, freq, $ + el_d, p0[tipo],p1[tipo],e0[tipo],e1[tipo],c0[tipo],c1[tipo],d0[tipo],d1[tipo],SNR0[tipo], SNR1[tipo], offset[tipo], fw_ratio0[tipo], sclab + endelse + endelse + endif + + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + + ; CUBIC FIT OF STACKED SCAN + Lfwhm=0 + Rfwhm=0 + + targetfit, numL,'stacked','cubic','Ch_0',tipo,list[i],list[i],Unit2,mycubc2j_0,err_mycubc2j_0,data[0].el,tau0L,ascissa,sommaL,ii,ff,x_mid,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Lfwhm, n_offL, offL, dum1, dum2, dum3, dum4, dum5, pstack, doplot + p0c[tipo]=dum1 + e0c[tipo]=dum2 + c0c[tipo]=dum3 + d0c[tipo]=dum4 + SNR0c[tipo]=dum5 + fw_ratio0c[tipo]=abs(Lfwhm/beam) + + targetfit, numR,'stacked','cubic','Ch_1',tipo,list[i],list[i],Unit2,mycubc2j_1,err_mycubc2j_1,data[0].el,tau0R,ascissa,sommaR,ii,ff,x_mid,Nsamples,sd*(1.+tipo*(1./cos(decs)-1.)),gpos,decsd,rasd, Rfwhm, n_offR, offR, dum1, dum2, dum3, dum4, dum5, pstack, doplot + p1c[tipo]=dum1 + e1c[tipo]=dum2 + c1c[tipo]=dum3 + d1c[tipo]=dum4 + SNR1c[tipo]=dum5 + fw_ratio1c[tipo]=abs(Rfwhm/beam) + + if doplot eq 'y' then begin + if (tipo eq 1) then begin + pstack.save, list[i]+'_stacked.pdf', /LANDSCAPE, RESOLUTION=300, /APPEND, /CLOSE + endif else begin + pstack.save, list[i]+'_stacked.pdf', /LANDSCAPE, RESOLUTION=300, /APPEND + endelse + endif + + if (n_offL+n_offR gt 0) then offsetc[tipo]=(offL+offR)/(n_offL+n_offR) else offsetc[tipo] = -99. + + el_d=data[0].el*180.0/!dpi + + if tipo eq 0 then sclab='DEC' + if tipo eq 1 then sclab='RA' + + if TPlike eq 2 then begin + printf, Unit11, FORMAT = '(a15,a7,D12.5,d8.1,13(" ",D8.4),A7)', target, hhmmss, data[0].time, freq, $ + el_d, p0c[tipo],p1c[tipo],e0c[tipo],e1c[tipo],c0c[tipo],c1c[tipo],d0c[tipo],d1c[tipo],SNR0c[tipo], SNR1c[tipo], offsetc[tipo], Lfwhm/beam, sclab + endif else begin + if c0[tipo] lt 0 or c1[tipo] lt 0 then begin ; here there are -99.000 values + printf, Unit11, FORMAT = '(a15,a7,D12.5,d8.1,d8.4,4(" ",D12.6),8(" ",D8.4),A7)', target, hhmmss, data[0].time, freq, $ + el_d, p0[tipo],p1[tipo],e0[tipo],e1[tipo],c0[tipo],c1[tipo],d0[tipo],d1[tipo],SNR0[tipo], SNR1[tipo], offset[tipo], fw_ratio0[tipo], sclab + endif else begin + printf, Unit11, FORMAT = '(a15,a7,D12.5,d8.1,d8.4,4(" ",E12.6),8(" ",D8.4),A7)', target, hhmmss, data[0].time, freq, $ + el_d, p0[tipo],p1[tipo],e0[tipo],e1[tipo],c0[tipo],c1[tipo],d0[tipo],d1[tipo],SNR0[tipo], SNR1[tipo], offset[tipo], fw_ratio0[tipo], sclab + endelse + endelse + endif + endfor + + ; If the pointing offsets are available in both directions, the amplitude measurement is compensated for them. + ; Otherwise, the measurement is rejected. + ; In this position also the FWHM check against HPBW - i.e. beam - is performed. A tolerance of FWHM/beam is here defined, the + ; measurements outside the consequent boundaries are rejected as well. + + fw_tol=0.3 ; we accept measurements where FWHM is within 30% of the nominal HPBW + + if (fitchoice eq 'linear') or (fitchoice eq 'both') then begin + + if (offset[0] eq -99. or offset[1] eq -99. or (abs(fw_ratio0[0]-1.0) ge fw_tol and abs(fw_ratio0[1]-1.0) ge fw_tol) or (abs(fw_ratio1[0]-1.0) ge fw_tol and abs(fw_ratio1[1]-1.0) ge fw_tol)) then begin + peak_cnt_0=-99.0 + peak_cnt_1=-99.0 + err_cnt_0=-99.0 + err_cnt_1=-99.0 + fluxdensity_0=-99.0 + fluxdensity_1=-99.0 + err_fluxdensity_0=-99.0 + err_fluxdensity_1=-99.0 + endif else begin + ; riscalo le ampiezze e gli errori della polarizzazione 0 per l'offset incrociato + p0[0]=p0[0]/exp(-(offset[1]*1.66/beamd)^2.) + p0[1]=p0[1]/exp(-(offset[0]*1.66/beamd)^2.) + e0[0]=e0[0]/exp(-(offset[1]*1.66/beamd)^2.) + e0[1]=e0[1]/exp(-(offset[0]*1.66/beamd)^2.) + ; riscalo le ampiezze e gli errori della polarizzazione 1 per l'offset incrociato + p1[0]=p1[0]/exp(-(offset[1]*1.66/beamd)^2.) + p1[1]=p1[1]/exp(-(offset[0]*1.66/beamd)^2.) + e1[0]=e1[0]/exp(-(offset[1]*1.66/beamd)^2.) + e1[1]=e1[1]/exp(-(offset[0]*1.66/beamd)^2.) + + peak_cnt_0 = tar_wmean(p0, e0, /nan) + peak_cnt_1 = tar_wmean(p1, e1, /nan) + w0=1/(e0/p0) + w1=1/(e1/p1) + err_cnt_0 = sqrt((w0[0]*e0[0])^2+(w0[1]*e0[1])^2)/(w0[0]+w0[1]) + err_cnt_1 = sqrt((w1[0]*e1[0])^2+(w1[1]*e1[1])^2)/(w1[0]+w1[1]) + ; err_cnt_0 = stddev(p0, /nan) + ; err_cnt_1 = stddev(p1, /nan) + + ; riscalo i flussi e gli errori della polarizzazione 0 per l'offset incrociato + c0[0]=c0[0]/exp(-(offset[1]*1.66/beamd)^2.) + c0[1]=c0[1]/exp(-(offset[0]*1.66/beamd)^2.) + d0[0]=d0[0]/exp(-(offset[1]*1.66/beamd)^2.) + d0[1]=d0[1]/exp(-(offset[0]*1.66/beamd)^2.) + ; riscalo i flussi e gli errori della polarizzazione 1 per l'offset incrociato + c1[0]=c1[0]/exp(-(offset[1]*1.66/beamd)^2.) + c1[1]=c1[1]/exp(-(offset[0]*1.66/beamd)^2.) + d1[0]=d1[0]/exp(-(offset[1]*1.66/beamd)^2.) + d1[1]=d1[1]/exp(-(offset[0]*1.66/beamd)^2.) + + ; XXX + fluxdensity_0 = tar_wmean(c0, d0, /nan) + fluxdensity_1 = tar_wmean(c1, d1, /nan) + w0=1/(d0/c0) + w1=1/(d1/c1) + err_fluxdensity_0 = sqrt((w0[0]*d0[0])^2+(w0[1]*d0[1])^2)/(w0[0]+w0[1]) + err_fluxdensity_1 = sqrt((w1[0]*d1[0])^2+(w1[1]*d1[1])^2)/(w1[0]+w1[1]) + + ; last check: when calibrators are not known, the resulting dummy flux values are to be reset to -99.00 + ; as the above offset compensation affects even them (and steers them a little bit from the dummy value) + if (c0[0] lt 0) or (c0[1] lt 0) then begin + fluxdensity_0 = -99.0 + err_fluxdensity_0 = -99.0 + endif + if (c1[0] lt 0) or (c1[1] lt 0) then begin + fluxdensity_1 = -99.0 + err_fluxdensity_1 = -99.0 + endif + endelse + + if TPlike eq 2 then begin + printf, Unit0, FORMAT = '(a15,a7,D12.5,d8.1,9(" ",D8.4))', sourcename, hhmmss, data[0].time, freq, $ + el_d, peak_cnt_0, peak_cnt_1,$ + err_cnt_0, err_cnt_1, fluxdensity_0, fluxdensity_1, $ + err_fluxdensity_0, err_fluxdensity_1 + ; note that SNR is not reported as it was not recalculated for the combined types + endif else begin + if fluxdensity_0 lt 0 or fluxdensity_1 lt 0 then begin ; here there are -99.000 values + printf, Unit0, FORMAT = '(a15,a7,D12.5,d8.1,D8.4,4(" ",D12.6),4(" ",D8.4))', sourcename, hhmmss, data[0].time, freq, $ + el_d, peak_cnt_0, peak_cnt_1,$ + err_cnt_0, err_cnt_1, fluxdensity_0, fluxdensity_1, $ + err_fluxdensity_0, err_fluxdensity_1 + endif else begin + printf, Unit0, FORMAT = '(a15,a7,D12.5,d8.1,D8.4,4(" ",E12.6),4(" ",D8.4))', sourcename, hhmmss, data[0].time, freq, $ + el_d, peak_cnt_0, peak_cnt_1,$ + err_cnt_0, err_cnt_1, fluxdensity_0, fluxdensity_1, $ + err_fluxdensity_0, err_fluxdensity_1 + endelse + endelse + endif + + + if (fitchoice eq 'cubic') or (fitchoice eq 'both') then begin + + if (offsetc[0] eq -99. or offsetc[1] eq -99. or (abs(fw_ratio0c[0]-1.0) ge fw_tol and abs(fw_ratio0c[1]-1.0) ge fw_tol) or (abs(fw_ratio1c[0]-1.0) ge fw_tol and abs(fw_ratio1c[1]-1.0) ge fw_tol)) then begin + peak_cnt_0=-99.0 + peak_cnt_1=-99.0 + err_cnt_0=-99.0 + err_cnt_1=-99.0 + fluxdensity_0=-99.0 + fluxdensity_1=-99.0 + err_fluxdensity_0=-99.0 + err_fluxdensity_1=-99.0 + endif else begin + ; riscalo le ampiezze e gli errori della polarizzazione 0 per l'offset incrociato + p0c[0]=p0c[0]/exp(-(offsetc[1]*1.66/beamd)^2.) + p0c[1]=p0c[1]/exp(-(offsetc[0]*1.66/beamd)^2.) + e0c[0]=e0c[0]/exp(-(offsetc[1]*1.66/beamd)^2.) + e0c[1]=e0c[1]/exp(-(offsetc[0]*1.66/beamd)^2.) + ; riscalo le ampiezze e gli errori della polarizzazione 1 per l'offset incrociato + p1c[0]=p1c[0]/exp(-(offsetc[1]*1.66/beamd)^2.) + p1c[1]=p1c[1]/exp(-(offsetc[0]*1.66/beamd)^2.) + e1c[0]=e1c[0]/exp(-(offsetc[1]*1.66/beamd)^2.) + e1c[1]=e1c[1]/exp(-(offsetc[0]*1.66/beamd)^2.) + + peak_cnt_0 = tar_wmean(p0c, e0c, /nan) + peak_cnt_1 = tar_wmean(p1c, e1c, /nan) + w0=1/(e0c/p0c) + w1=1/(e1c/p1c) + err_cnt_0 = sqrt((w0[0]*e0[0])^2+(w0[1]*e0[1])^2)/(w0[0]+w0[1]) + err_cnt_1 = sqrt((w1[0]*e1[0])^2+(w1[1]*e1[1])^2)/(w1[0]+w1[1]) + ; err_cnt_0 = stddev(p0, /nan) + ; err_cnt_1 = stddev(p1, /nan) + + ; riscalo i flussi e gli errori della polarizzazione 0 per l'offset incrociato + c0c[0]=c0c[0]/exp(-(offsetc[1]*1.66/beamd)^2.) + c0c[1]=c0c[1]/exp(-(offsetc[0]*1.66/beamd)^2.) + d0c[0]=d0c[0]/exp(-(offsetc[1]*1.66/beamd)^2.) + d0c[1]=d0c[1]/exp(-(offsetc[0]*1.66/beamd)^2.) + ; riscalo i flussi e gli errori della polarizzazione 1 per l'offset incrociato + c1c[0]=c1c[0]/exp(-(offsetc[1]*1.66/beamd)^2.) + c1c[1]=c1c[1]/exp(-(offsetc[0]*1.66/beamd)^2.) + d1c[0]=d1c[0]/exp(-(offsetc[1]*1.66/beamd)^2.) + d1c[1]=d1c[1]/exp(-(offsetc[0]*1.66/beamd)^2.) + + fluxdensity_0 = tar_wmean(c0c, d0c, /nan) + fluxdensity_1 = tar_wmean(c1c, d1c, /nan) + w0=1/(d0/c0) + w1=1/(d1/c1) + err_fluxdensity_0 = sqrt((w0[0]*d0[0])^2+(w0[1]*d0[1])^2)/(w0[0]+w0[1]) + err_fluxdensity_1 = sqrt((w1[0]*d1[0])^2+(w1[1]*d1[1])^2)/(w1[0]+w1[1]) + + ; last check: when calibrators are not known, the resulting dummy flux values are to be reset to -99.00 + ; as the above offset compensation affects even them (and steers them a little bit from the dummy value) + if (c0[0] lt 0) or (c0[1] lt 0) then begin + fluxdensity_0 = -99.0 + err_fluxdensity_0 = -99.0 + endif + if (c1[0] lt 0) or (c1[1] lt 0) then begin + fluxdensity_1 = -99.0 + err_fluxdensity_1 = -99.0 + endif + + endelse + + if TPlike eq 2 then begin + printf, Unit10, FORMAT = '(a15,a7,D12.5,d8.1,9(" ",D8.4))', sourcename, hhmmss, data[0].time, freq, $ + el_d, peak_cnt_0, peak_cnt_1,$ + err_cnt_0, err_cnt_1, fluxdensity_0, fluxdensity_1, $ + err_fluxdensity_0, err_fluxdensity_1 + ; note that SNR is not reported as it was not recalculated for the combined types + endif else begin + if fluxdensity_0 lt 0 or fluxdensity_1 lt 0 then begin ; here there are -99.000 values + printf, Unit10, FORMAT = '(a15,a7,D12.5,d8.1,D8.4,4(" ",D12.6),4(" ",D8.4))', sourcename, hhmmss, data[0].time, freq, $ + el_d, peak_cnt_0, peak_cnt_1,$ + err_cnt_0, err_cnt_1, fluxdensity_0, fluxdensity_1, $ + err_fluxdensity_0, err_fluxdensity_1 + endif else begin + printf, Unit10, FORMAT = '(a15,a7,D12.5,d8.1,D8.4,4(" ",E12.6),4(" ",D8.4))', sourcename, hhmmss, data[0].time, freq, $ + el_d, peak_cnt_0, peak_cnt_1,$ + err_cnt_0, err_cnt_1, fluxdensity_0, fluxdensity_1, $ + err_fluxdensity_0, err_fluxdensity_1 + endelse + endelse + endif + + endfor + + ; closing all txt output files + close, /ALL + + ; freeing logical units + if fitchoice eq 'linear' or fitchoice eq 'both' then begin + free_lun, Unit0, /force + free_lun, Unit1, /force + endif + if fitchoice eq 'cubic' or fitchoice eq 'both' then begin + free_lun, Unit10, /force + free_lun, Unit11, /force + endif + + free_lun, Unit2, /force + + + print, ' ' + print, '**************************' + print, '*** RUNTARGET is DONE ****' + print, '*** GO LOOK at RESULTS ***' + print, '**************************' + print, ' ' + + return +end + diff --git a/separate.pro b/separate.pro new file mode 100644 index 0000000..05f2d56 --- /dev/null +++ b/separate.pro @@ -0,0 +1,197 @@ + +; This procedure separates the acquisitions into sub-folders, depending on the frequency band and the +; calibrator/skydip/target nature of the scan. +; +; The usable RF bands are hard-coded. +; The original version of this procedure considers the following bands and frequency boundaries (MHz): +; +; pbounds=[300,1000] +; lbounds=[1000,2000] +; sbounds=[2000,4000] +; cbounds=[4000,8000] +; xbounds=[8000,10000] +; kbounds=[18000,26500] +; +; If observing with more than one sub-band in the K-band RF range, specify /multiK as in +; +; IDL> separa, /multiK (in this case, kbounds=[18000,20000,23000,26500]) +; +; Default flux density calibrators are considered to be 3C48, 3C123, 3C286, 3C295, NGC7027. +; When wanting to specify a custom list of calibrators, provide an external text file +; containing as columns the source names and flux densities (the same to be used when running "runcalib"), +; and then use +; +; IDL> separa, /extlist +; +; Latest update: Oct 18th, 2017 + + +pro separate, multiK=multiK, extlist=extlist + + sp=path_sep() + pin=dialog_pickfile(/DIRECTORY, title='Select main folder, containing the YYYMMDD data folders') + + if not keyword_set(multiK) then Kband='single' else Kband='multi' + + homecals=['3C48','3C123','3C286','3C295','NGC7027'] ; default flux density calibrators + + if keyword_set(extlist) then begin + myextlist=dialog_pickfile(TITLE='Please pick EXTERNAL LIST of calibrators') + readcol, myextlist, extcalname, extcalflux, format='(A,D)', /SILENT + calname=strupcase(extcalname) + endif else begin + myextlist=' ' + calname=strupcase(homecals) + endelse + + ; In this module the frequency bands are defined + ; Users might want to add/edit some + ; *********************************************** + + pbounds=[300,1000] + lbounds=[1000,2000] + sbounds=[2000,4000] + cbounds=[4000,8000] + xbounds=[8000,10000] + + ; NOTICE: in case of editing, edit the boundaries implementation as well! Search for string "edit here" in this code. + + if Kband eq 'single' then begin + bands=['P','L','S','C','X','K'] + kbounds=[18000,26500] + endif else begin + bands=['P','L','S','C','X','K-low','K-mid','K-hi'] + kbounds=[18000,20000,23000,26500] ; definition of sub-bands boundaries for K-band + endelse + + ; *********************************************** + + bandnum=n_elements(bands) + mybounds=make_array(1,bandnum,/integer) + + + ; Checking whether the folder for DAT files (produced by the CalibrationTool) exists. + ; If not, it is created. + ans0=file_test(pin+'DATS',/DIRECTORY) + if ans0 eq 0 then begin + print, 'Creating folder for dat files' + file_mkdir, pin+sp+'DATS' + endif + + for b=0, bandnum-1 do begin + + ; Checking whether the band-dependant and type-dependant folders alreasy exist. + ; If not, they are created. + + ans1=file_test(pin+bands[b],/DIRECTORY) + + if ans1 eq 0 then begin + print, 'Creating folder for '+bands[b]+'-band files' + file_mkdir, pin+bands[b] + endif + + ans2=file_test(pin+bands[b]+sp+'CALIBRATORS',/DIRECTORY) + ans3=file_test(pin+bands[b]+sp+'TARGETS',/DIRECTORY) + ans4=file_test(pin+bands[b]+sp+'SKYDIPS',/DIRECTORY) + + if ans2 eq 0 then begin + ; print, 'Creating folder for '+bands[b]+'-band CALIBRATORS files' + file_mkdir, pin+bands[b]+sp+'CALIBRATORS' + endif + if ans3 eq 0 then begin + ; print, 'Creating folder for '+bands[b]+'-band TARGETS files' + file_mkdir, pin+bands[b]+sp+'TARGETS' + endif + if ans4 eq 0 then begin + ; print, 'Creating folder for '+bands[b]+'-band SKYDIP files' + file_mkdir, pin+bands[b]+sp+'SKYDIPS' + endif + + endfor + + + list = FILE_SEARCH(pin+'20*', COUNT=number, /TEST_DIRECTORY) + + for i =0,number-1 do begin ; loop over 'number' days + print, list[i] + datlist = file_search(list[i]+sp+'*.dat', COUNT=datnum) + if datnum ne 0 then begin + for j=0,datnum-1 do begin + file_move, datlist[j], pin+sp+'DATS'+sp, /OVERWRITE + endfor + endif + + sublist = FILE_SEARCH(list[i]+sp+'20*', /TEST_DIRECTORY, COUNT=subnum) + for j=0,subnum-1 do begin ; loop over 'subnum' scans + ; print, sublist[j] + subsublist = FILE_SEARCH(sublist[j]+sp+'*.fits', COUNT=subsubnum) + ; obtaining info from the first subscan (i.e. the first FITS in the folder) + RFinfo=mrdfits(subsublist[0],2,/silent) + sky_bandwidth=RFinfo[0].bandWidth + sky_freq=RFinfo[0].frequency+sky_bandwidth/2.0 + mainh=mrdfits(subsublist[0],0,head0,/silent) + findh0key, head0, '*SOURCE*', keylab, target, infolab, tarflag + + ; elaborating the target name string + target=strcompress(target, /remove_all) + cleantarget=strsplit(target,"'",/extract) + target=cleantarget[0] + + ; directing the file to the proper location + type='TARGETS'+sp ; default choice + + checkcal=where(target eq calname, isitcal) + if isitcal ne 0 then begin + type = 'CALIBRATORS'+sp ; overriding the default choice: the source is a calibrator + if subsubnum eq 1 then begin ; to handle cases when SKYDIPS have been mis-named + data=mrdfits(subsublist[0],4,/SILENT) + span=abs(data[0].el-data[-1].el) + if span gt 45/180*!dpi then type = 'SKYDIPS'+sp + endif + endif + + if (strmatch(sublist[j],'*sky*',/FOLD_CASE)) then type = 'SKYDIPS'+sp ; overriding the default choice: the file is a SKYDIP + + ; actually moving the folders + for b=0, bandnum-1 do begin + case bands[b] of + ; edit here in case you have added more bands + 'P': mybounds=pbounds + 'L': mybounds=lbounds + 'S': mybounds=sbounds + 'C': mybounds=cbounds + 'X': mybounds=xbounds + 'K': mybounds=kbounds + endcase + intervals=n_elements(mybounds)-1 + for int=0, intervals-1 do begin + if (sky_freq ge mybounds[int] and sky_freq lt mybounds[int+1]) then begin + outband = pin+bands[b]+sp + file_move, sublist[j], outband+type, /OVERWRITE + break + endif + endfor + endfor + + endfor ; ending cycle on scans + + endfor; ending cycle on days + + ; removing empty/unused folders (the non-empty ones will be spared!) + for b=0, bandnum-1 do begin + file_delete, pin+bands[b]+sp+'CALIBRATORS',/quiet + file_delete, pin+bands[b]+sp+'TARGETS',/quiet + file_delete, pin+bands[b]+sp+'SKYDIPS',/quiet + file_delete, pin+bands[b],/quiet + endfor + file_delete, pin+sp+'DATS',/quiet + file_delete, pin+sp+'20*',/quiet + + + print, '***********' + print, 'PROGRAM END' + print, '***********' + return +end + + diff --git a/srcaver.pro b/srcaver.pro new file mode 100644 index 0000000..2ad3581 --- /dev/null +++ b/srcaver.pro @@ -0,0 +1,68 @@ +pro srcaver + +; Program to average the LCP and RCP measurements produced by CAP, performing a source-based weighted mean. +; Just select runtarget's output file from the graphic interface appearing after launching: +; +; IDL> srcaver +; +; Notice: results are written in a file named "Combined_measurements.txt" +; In case of multiple runs, the generated files are automatically numbered with a sequential suffix. +; The produced file includes the full path to the original file given in input. +; +; Last edited: Oct 16, 2017 + +fin=dialog_pickfile() + +close, /all + +readcol, fin, src, hhmmss, mjd, freq, elev, A0, A1, eA0, eA1, f0, f1, ef0, ef1, format='(A,D,D,D,D,D,D,D,D,D,D,D,D)' + +src_list=src[UNIQ(src, sort(src))] + +num=n_elements(src_list) + +aver_f0=dblarr(num) +aver_f1=dblarr(num) +aver_ef0=dblarr(num) +aver_ef1=dblarr(num) +combined_f=dblarr(num) +combined_ef=dblarr(num) + +outf=file_search('Combined_measurements*.txt', count=numf) +if numf eq 0 then begin + outfile='Combined_measurements.txt' +endif else begin + outfile='Combined_measurements_'+strcompress(string(numf+1),/remove_all)+'.txt' +endelse + +openw, 1, outfile +printf, 1, 'Elaboration of CAP results contained in: ' +printf, 1, fin +printf, 1, ' ' +printf, 1, ' Name Freq F0 eF0 F1 eF1 F eF' +printf, 1, ' MHz Jy Jy Jy Jy Jy Jy' + + +for i=0, num-1 do begin + + thissrc=where(src eq src_list[i]) + w0=1./ef0[thissrc]^2 + w1=1./ef1[thissrc]^2 + aver_f0[i]=total(f0[thissrc]/ef0[thissrc]^2)/total(w0) + aver_f1[i]=total(f1[thissrc]/ef1[thissrc]^2)/total(w1) + aver_ef0[i]=sqrt(1/total(w0)) + aver_ef1[i]=sqrt(1/total(w1)) + + combined_f[i]=(aver_f0[i]/aver_ef0[i]^2+aver_f1[i]/aver_ef1[i]^2)/(1./aver_ef0[i]^2+1./aver_ef1[i]^2) + combined_ef[i]=sqrt(1/(1./aver_ef0[i]^2+1./aver_ef1[i]^2)) + + printf, 1, src_list[i], freq[0], aver_f0[i], aver_ef0[i], aver_f1[i], aver_ef1[i], combined_f[i], combined_ef[i], format='(A13,1X,I5,1X,D6.3,1X,D5.3,1X,D6.3,1X,D5.3,1X,D6.3,1X,D5.3)' + +endfor + +close, 1 + +print, '****** DONE ******' + +return +end \ No newline at end of file diff --git a/targetfit.pro b/targetfit.pro new file mode 100644 index 0000000..03d9c3d --- /dev/null +++ b/targetfit.pro @@ -0,0 +1,172 @@ +pro targetresiduals, X,Y,datafit,Unit2,peak,err, x_peak, rasd, decsd, tipo, SNR, peak_cnt, err_cnt, Level, n_off, off, residual + + ; minor section to compute the rms-noise and the signal-to-noise ratio + + fit = LINFIT(X[0:20], Y[0:20], YFIT=datalinfit) + subtracted = Y[0:20]-datalinfit + mom = moment(subtracted) + noise = sqrt(mom[1]) + residual = Y-datafit + mom = moment(residual) + resnoise = sqrt(mom[1]) + PRINTF, Unit2, "noise: ", noise, " counts; noise residuals: ", resnoise, " counts" + ; ;;print, "noise: ", noise, " counts; noise residuals: ", resnoise, " counts" + SNR=peak/noise + PRINTF, Unit2, "SNR: ", SNR + peak_cnt = peak + err_cnt = err + Level = (Y[-1] + Y[0])/2. + + case tipo of + 0: begin + meanlat=x_peak/180d*!dpi ; radians + off=(x_peak-decsd) ; degrees + end + 1: begin + meanlat=decsd ; radians + off=(x_peak-rasd)*cos(meanlat) ; degrees + end + endcase + n_off=1 + + return +end + + + +pro getflux, peak_cnt, err_cnt, tau0, datael, cnt2Jy, err_cnt2Jy, outflux, err_outflux, Unit2, Level + + ; Minor section to computate flux density (application of cnt2jy factor) + + peak_cnt=peak_cnt*exp(mean(tau0)/sin(datael)) + err_cnt=err_cnt*exp(mean(tau0)/sin(datael)) + outflux = peak_cnt*cnt2Jy + err_outflux = outflux*sqrt((err_cnt/peak_cnt)^2+(err_cnt2Jy/cnt2Jy)^2) + Level=Level*cnt2Jy + PRINTF, Unit2, " " + PRINTF, Unit2, "outflux = (", outflux, "+-", err_outflux, ") [Jy]" + PRINTF, Unit2, " " + PRINTF, Unit2, "Level = ", Level, " [Jy]" + PRINTF, Unit2, " " + return +end + + + +pro targetfit, scanflag,stacflag,polyflag,section,tipo,allpath,namefile,Out3,cnt2Jy,err_cnt2Jy,datael,tau0, xx,yy,ii,ff,x_mid,Nsamples,sd_sub,gpos,decsd,rasd, fwhm, n_off, off, p, e, c, d, SNR, plo, doplot + +; Main procedure, devoted to the fitting operations, both in the linear and cubic cases, for calibrators. +; +; Authors: Marcello Giroletti, Simona Righini +; Last edited: Oct 6, 2017 +; + + ; scanflag = scan buono (=1) o cattivo (=0) + ; stacflag = single o stacked scan + ; polyflag = gaussiana più grado del polinomio (primo='linear', terzo='cubic') + ; section = Ch_0 o Ch_1 + ; tipo = flag di direzione - convertito qua sotto in direflag + if (tipo eq 1) then direflag = 'RA' else direflag = 'DEC' + + ; common logistics + sep=path_sep() + if scanflag eq 0 then begin + + n_off = 0 + off = 0. + p = -99. + e = -99. + c = -99. + d = -99. + SNR = -99. + if doplot eq 'y' then dataplot, 0, [0,1], [0,1], namefile, allpath, section, 0, 0, 0, 0, 0, direflag, polyflag, plo, stacflag + + endif else begin + + decs=decsd*!dpi/180.0 + media = (yy[ii]+yy[ff])/2. + ; define new coordinates in which the source is centred at x=0 and the baseline is about 0 at the source position + X=xx[ii:ff]-gpos + Y=yy[ii:ff]-media + + if (polyflag eq 'linear') then begin + + datafit = GAUSSFIT(X, Y, A_fit, ESTIMATES=[(Y[x_mid]),0.0,sd_sub,(X[0]*Y[-1]-X[-1]*Y[0])/(X[0]-X[-1]),(Y[-1]-Y[0])/(X[-1]-X[0])], NTERMS=5, SIGMA=A_sig) +; print, ' Ampl(cnt) Peak_pos(deg) sigma(deg) q(cnt) m(cnt/deg)' +; print, 'Estimates =',(Y[x_mid]),0.0,sd_sub,(X[0]*Y[-1]-X[-1]*Y[0])/(X[0]-X[-1]),(Y[-1]-Y[0])/(X[-1]-X[0]) +; print, 'Fit =',A_fit +; print, 'Sigma =',A_sig + if (tipo eq 0) then begin + fwhm=2*SQRT(2*ALOG(2))*A_fit[2]*60.0 ; in arcmin, if fit is done on coordinates in degrees ; *(speed/60.)*60. to be used if done in samples + endif else begin + fwhm=2*SQRT(2*ALOG(2))*A_fit[2]*60.0*cos(decs) ; in arcmin, if fit is done on coordinates in degrees ; *(speed/60.)*60. to be used if done in samples + endelse + x_peak=A_fit[1]+gpos ; coordinate corresponding to gaussian peak + PRINTF, Out3, " GAUSSIAN + LINEAR " + PRINTF, Out3, section+" Gaussian fit -- amplitude:", A_fit[0], " +- ", A_sig[0], " counts" + PRINTF, Out3, section+" Gaussian fit -- sigma: ", A_fit[2], " +- ", A_sig[2], " [s]" + PRINTF, Out3, " corresponding to FWHM: ", fwhm, " [arcmin]" + + targetresiduals, X,Y,datafit,Out3,A_fit[0],A_sig[0], x_peak, rasd, decsd, tipo, SNR, peak_cnt, err_cnt, Level, n_off, off, residual + + ; computing and accounting for the RMS of the residuals in the overall amplitude error estimate + nres=n_elements(residual) + resrange=ceil(nres/5.0) + rescut=[residual[0:resrange],residual[-1*resrange,-1]] ; avoiding the central part of the subscan, where artifacts can be present due to sidelobes + resstat=moment(rescut) + res_rms=sqrt(resstat[1]) + err_cnt=sqrt(err_cnt^2+res_rms^2) ; updated error for the amplitude measurement + + getflux, peak_cnt, err_cnt, tau0, datael, cnt2Jy, err_cnt2Jy, outflux, err_outflux, Out3, Level + if doplot eq 'y' then dataplot, 1, X, Y, namefile, allpath, section, datafit, cnt2Jy, Level, SNR, residual, direflag, polyflag, plo, stacflag + + endif else begin + + ; estimates for fit with A0*exp(-(X-A1)^2/(2*A2))+A3+A4*X+A5*X^2+A6*X^3; ie cubic+gaussian + A = [Y[x_mid],0.0,sd_sub,(X[0]*Y[-1]-X[-1]*Y[0])/(X[0]-X[-1]),(Y[-1]-Y[0])/(X[-1]-X[0]),0.0,0.0] + A_GUESS = A + ; compute the parameters without weights (DUM_W not set) + datafit = CURVEFIT(X,Y,DUM_W, A, SIGMA, FITA=[1,1,1,1,1,1,1], FUNCTION_NAME='tar_gfunct', /DOUBLE, STATUS=suc_fit) +; print, ' A0=Ampl(cnt) A1=Peak_pos(deg) A2=sigma(deg) A3(cnt) A4(cnt/deg) A5(cnt/deg2) A6(cnt/deg3)' +; print, 'Function estimates: ', A_GUESS +; print, 'Function parameters: ', A +; print, 'Sigma parameters: ', SIGMA + if (tipo eq 0) then begin + fwhm=2*SQRT(2*ALOG(2))*A[2]*60.0 ; in arcmin, if fit is done on coordinates in degrees ; *(speed/60.)*60. to be used if done in samples + endif else begin + fwhm=2*SQRT(2*ALOG(2))*A[2]*60.0*cos(decs) ; in arcmin, if fit is done on coordinates in degrees ; *(speed/60.)*60. to be used if done in samples + endelse + x_peak=A[1]+gpos ; coordinate corresponding to gaussian peak + PRINTF, Out3, " GAUSSIAN + CUBIC " + PRINTF, Out3, section+" Gaussian fit -- amplitude:", A[0], " +- ", SIGMA[2], " counts" + PRINTF, Out3, section+" Gaussian fit -- sigma: ", A[0], " +- ", SIGMA[0], " [s]" + PRINTF, Out3, " corresponding to FWHM: ", fwhm, " [arcmin]" + + targetresiduals, X,Y,datafit,Out3,A[0],A[2], x_peak, rasd, decsd, tipo, SNR, peak_cnt, err_cnt, Level, n_off, off, residual + + ; computing and accounting for the RMS of the residuals in the overall amplitude error estimate + nres=n_elements(residual) + resrange=ceil(nres/5.0) + rescut_b=[residual[0:resrange],residual[-1*resrange,-1]] ; avoiding the central part of the subscan, where artifacts can be present due to sidelobes + resstat=moment(rescut_b) + res_rms=sqrt(resstat[1]) + err_cnt=sqrt(err_cnt^2+res_rms^2) ; updated error for the amplitude measurement + + getflux, peak_cnt, err_cnt, tau0, datael, cnt2Jy, err_cnt2Jy, outflux, err_outflux, Out3, Level + if doplot eq 'y' then dataplot, 1, X, Y, namefile, allpath, section, datafit, cnt2Jy, Level, SNR, residual, direflag, polyflag, plo, stacflag + + endelse + + el_d=datael*180.0/!dpi + p=peak_cnt + e=err_cnt + c=outflux + d=err_outflux + + endelse + + return + +end + +