From 776c4383509e724b08cf1f8431d93a43574c2e36 Mon Sep 17 00:00:00 2001 From: SimonaRighini Date: Tue, 24 Oct 2017 14:35:33 +0200 Subject: [PATCH] First release, shared risk --- c2jtimeline.pro | 419 ++++++++++++++++ calflux.pro | 337 +++++++++++++ calibfit.pro | 179 +++++++ dataplot.pro | 54 ++ findh0key.pro | 33 ++ newflag.pro | 490 +++++++++++++++++++ runcalib.pro | 1247 +++++++++++++++++++++++++++++++++++++++++++++++ runtarget.pro | 1242 ++++++++++++++++++++++++++++++++++++++++++++++ separate.pro | 197 ++++++++ srcaver.pro | 68 +++ targetfit.pro | 172 +++++++ 11 files changed, 4438 insertions(+) create mode 100644 c2jtimeline.pro create mode 100644 calflux.pro create mode 100644 calibfit.pro create mode 100644 dataplot.pro create mode 100644 findh0key.pro create mode 100644 newflag.pro create mode 100644 runcalib.pro create mode 100644 runtarget.pro create mode 100644 separate.pro create mode 100644 srcaver.pro create mode 100644 targetfit.pro 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 + +