; The only required parameter is the scan number. This can be either scan ; in the sequence of two nod scans. Arguments to identify the IF number, ; polarization number, and feed number are optional. The default feed ; number (0) is the lowest numbered FEED found in the data. This is ; assumed to have been the tracking (source) feed for the first scan. ; This can be overridden by specifying fdnum=1 which ; will use the second FEED in the data as the tracking feed. ; Data with more than two feeds will be rejected by the procedure. ;
; The procedure calculates Tsys based on the Tcal values, and cal_on/cal_off phases. ; The user can override this calculation by entering a zenith system temperature. ; The procedure will then correct the Tsys for the observed elevation. ; If the data are calibrated to Ta* or Jy, then additional parameters are used. ; A zenith tau value may be specified, and an aperture efficiency may be specified. ; The user is *strongly* encouraged to enter values for these calibration parameters, ; but they will be estimated if none is provided. ;
; A parameter called smthoff can be used to smooth the reference spectrum prior ; to calibration. In certain cases this can improve the signal to noise ratio, ; but it may degrade baseline shapes and artificially emphasize spectrometer ; glitches. Use with care. A value of smthoff=16 is often a good choice. ; ;
; All averaging is done using dcaccum ; and accumave. This means that integrations ; are weighted using the default weight given in dcaccum. To turn off this weighting ; so that all integrations have equal weight in the average, use the /noweight flag. ; ; @param scan {in}{required}{type=integer} M&C scan number ; @keyword ifnum {in}{optional}{type=integer} IF number (starting with 0) ; @keyword intnum {in}{optional}{type=integer} Integration number (default=all) ; @keyword plnum {in}{optional}{type=integer} polarization number (default 0) ; @keyword fdnum {in}{optional}{type=integer} feed number (default 0) ; @keyword tau {in}{optional}{type=float} tau at zenith ; @keyword tsys {in}{optional}{type=float} tsys at zenith ; @keyword ap_eff {in}{optional}{type=float} aperture efficiency ; @keyword smthoff {in}{optional}{type=integer} smooth factor for reference spectrum ; @keyword units {in}{optional}{type=string} takes the value 'Jy', 'Ta', or 'Ta*' (default Ta) ; @keyword noweight {in}{optional}{type=flag} if set, the default ; weighting is turned off and all integrations are averaged with the ; same weight (1.0). ; ; @version $Id: getnod.pro,v 1.16 2005/05/30 04:03:43 bgarwood Exp $ ;- pro getnod,scan,ifnum=ifnum,intnum=intnum,plnum=plnum,fdnum=fdnum,tau=tau,$ tsys=tsys,ap_eff=ap_eff,smthoff=smthoff,units=units,noweight=noweight compile_opt idl2 if (n_elements(scan) eq 0) then begin message, 'The scan number is required', /info return endif if not !g.lineio->is_data_loaded() then begin message,'No line data is attached yet, use filein or dirin.',/info return endif ; set defaults if n_elements(ifnum) eq 0 then ifnum = 0 if n_elements(plnum) eq 0 then plnum = 0 if n_elements(fdnum) eq 0 then fdnum = 0 if n_elements(tsys) eq 0 then tsys = 1 if n_elements(smthoff) eq 0 then smthoff = 1 if n_elements(units) eq 0 then units = 'Ta' if keyword_set(noweight) then weight = 1.0 ; else undefined and use default weight ; Check if scan number is valid validscans = get_scan_numbers() if total(validscans eq scan, /integer) gt 1 then $ message,"Warning: More than one scan with the scan number is in the data file.",/info if total(validscans eq scan, /integer) eq 0 then begin message,"That scan is not available.",/info return end info = scan_info(scan) if info.procedure ne 'Nod' then begin message,"Cannot handle this scan: Procedure = " + strcompress(info.procedure,/remove_all),/info return end if info.subscan ne 0 then begin if info.subscan ne 1 then begin sSub = strcompress(string(info.subscan),/remove_all) message,"More than two subscans in this procedure, at least : " + sSub, /info return endif scan = scan - 1 if total(validscans eq scan,/integer) gt 1 then $ message,"Warning: First scan in procedure appears more than once in the data file.",/info if total(validscans eq scan,/integer) eq 0 then begin message,"First scan in this procedure is missing.",/info return endif info = scan_info(scan) if info.subscan ne 0 then begin message,"The first scan in this procedure can not be found.",/info return endif endif else begin ; check that the other scan exists if total(validscans eq (scan+1),/integer) gt 1 then $ message,"Warning: second scan number in procedure appears more than once in the data file.",/info if total(validscans eq (scan+1),/integer) eq 0 then begin message,"The second scan in the procedure is not available.",/info return endif info2 = scan_info(scan+1) if info2.subscan ne 1 then begin message,"The second scan in the procedure can not be found.",/info return endif endelse if info.n_switching_states ne 2 then begin message,string("This does not appeat to be total power data."),/info return endif if ifnum lt 0 or ifnum gt (info.n_ifs-1) then begin sifnum = strcompress(string(ifnum),/remove_all) snif = strcompress(string(info.n_ifs),/remove_all) message,"Illegal IF identifier: " + sifnum + ". This scan has " + snif + " IFs, zero-indexed.", /info return endif if info.n_feeds ne 2 then begin message,"This scan does not have two feeds. n_feeds = "+strcompress(string(info.n_feeds),/remove_all), /info return endif if fdnum lt 0 or fdnum gt 1 then begin message,"Invalid feed: " + strcompress(string(fdnum),/remove_all) + ". Feed must be 0 or 1",/info return endif fdnum2 = (fdnum eq 0) ? 1 : 0 if plnum lt 0 or plnum gt 1 then begin spol = strcompress(string(plnum),/remove_all) snpol = strcompress(string(info.n_polarizations),/remove_all) message, "Invalid polarization identifier: " + spol + ". This scan has " + snpol + " polarizations, zero-indexed.", /info message, string("Invalid polarization identifier: ",plnum), /info return endif if tsys lt 0.0 then begin message, 'Invalid tsys value', /info return endif if smthoff lt 1 or smthoff gt info.n_channels/4 then begin message, 'Invalid smthoff value', /info return endif thispol = info.polarizations[plnum] thisfeed = info.feeds[fdnum] thisfeed2 = info.feeds[fdnum2] ; Retrieve all the data necessary to satisfy this request singleInt = n_elements(intnum) eq 1 expectedCount = singleInt ? 1 : info.n_integrations if singleInt then begin if intnum ge 0 and intnum le (info.n_integrations-1) then begin data = !g.lineio->get_spectra(count,scan=[scan,scan+1],feed=[thisfeed,thisfeed2],ifnum=ifnum,pol=thispol,int=intnum) endif else begin message,"Integration number out of range", /info return endelse endif else begin data = !g.lineio->get_spectra(count,scan=[scan,scan+1],feed=[thisfeed,thisfeed2],ifnum=ifnum,pol=thispol) endelse if count le 0 then begin message,"No data found, this should never happen, can not continue.",/info return endif s1_b1_off = where(data.scan_number eq scan and data.feed eq thisfeed and data.cal_state eq 0,s1_b1_off_count) s1_b1_on = where(data.scan_number eq scan and data.feed eq thisfeed and data.cal_state eq 1,s1_b1_on_count) s2_b1_off = where(data.scan_number eq (scan+1) and data.feed eq thisfeed and data.cal_state eq 0,s2_b1_off_count) s2_b1_on = where(data.scan_number eq (scan+1) and data.feed eq thisfeed and data.cal_state eq 1,s2_b1_on_count) s1_b2_off = where(data.scan_number eq scan and data.feed eq thisfeed2 and data.cal_state eq 0,s1_b2_off_count) s1_b2_on = where(data.scan_number eq scan and data.feed eq thisfeed2 and data.cal_state eq 1,s1_b2_on_count) s2_b2_off = where(data.scan_number eq (scan+1) and data.feed eq thisfeed2 and data.cal_state eq 0,s2_b2_off_count) s2_b2_on = where(data.scan_number eq (scan+1) and data.feed eq thisfeed2 and data.cal_state eq 1,s2_b2_on_count) if (s1_b1_off_count ne expectedCount or $ s1_b1_off_count ne s1_b1_on_count or s1_b1_off_count ne s2_b1_off_count or s1_b1_off_count ne s2_b1_on_count or $ s1_b1_off_count ne s1_b2_on_count or s1_b1_off_count ne s2_b2_off_count or s1_b1_off_count ne s2_b2_on_count or $ s1_b1_off_count ne s1_b2_off_count) then begin message,"Unexpected number of spectra retrieved for some or all of the switching phases, can not continue.",/info data_free, data return endif ; copy first element into !g.s[0] to hold result old_frozen = !g.frozen freeze set_data_container,data[0] if n_elements(ap_eff) eq 0 then ap_eff = get_ap_eff(data[0].observed_frequency/1.0e9,data[0].elevation) if ap_eff lt 0.0 or ap_eff gt 1.0 then begin message, 'Invalid ap_eff value - it should be between 0 and 1', /info if old_frozen eq 1 then freeze else unfreeze data_free, data return endif if n_elements(tau) eq 0 then tau = get_tau(data[0].observed_frequency/1.0e9) if tau lt 0.0 or tau gt 1.0 then begin message, 'Invalid tau value - it should be between 0 and 1', /info if old_frozen eq 1 then freeze else unfreeze data_free, data return endif if singleInt then begin donodint,data[s1_b1_off],data[s1_b1_on], data[s2_b1_off],data[s2_b1_on],$ data[s1_b2_off],data[s1_b2_on], data[s2_b2_off],data[s2_b2_on],$ tau,tsys,ap_eff,smthoff,units,ret_tsys,ret_tau endif else begin thisaccum = {accum_struct} for i = 0,(info.n_integrations-1) do begin donodint,data[s1_b1_off[i]],data[s1_b1_on[i]],data[s2_b1_off[i]],data[s2_b1_on[i]],$ data[s1_b2_off[i]],data[s1_b2_on[i]],data[s2_b2_off[i]],data[s2_b2_on[i]],$ tau,tsys,ap_eff,smthoff,units,ret_tsys,ret_tau dcaccum,thisaccum,!g.s[0],weight=weight end naccum1 = thisaccum.n accumave,thisaccum, thisave, /quiet if naccum1 ne info.n_integrations then begin message,'unexpected problems in obtaining the average over integrations, can not continue',/info ; clean up if old_frozen eq 1 then freeze else unfreeze data_free, thisave accumclear, thisaccum data_free, data return endif set_data_container,thisave accumclear, thisaccum data_free, thisave endelse !g.s[0].exposure *= 2.0 ; required for nod procedure !g.s[0].duration *= 2.0 ; required for nod procedure if units eq 'Jy' then $ print,scan,ret_tau,ap_eff,ret_tsys[0],ret_tsys[1],ret_tsys[2],ret_tsys[3], $ format='(i5," units: Jy tau: ",f5.3," ap_eff: ",f5.3, " Tsys:",f7.2, f7.2, f7.2, f7.2)' if units eq 'Ta*' then $ print,scan,ret_tau,ret_tsys[0],ret_tsys[1],ret_tsys[2],ret_tsys[3], $ format='("Scan: ",i5," units: Ta* (K) tau: ",f5.3, " Tsys: ",f7.2, f7.2, f7.2, f7.2)' if units eq 'Ta' then $ print,scan,ret_tsys[0],ret_tsys[1],ret_tsys[2],ret_tsys[3], $ format='("Scan: ",i5," units: Ta (K) Tsys: ",f7.2, f7.2, f7.2, f7.2)' data_free, data if old_frozen eq 1 then freeze else unfreeze if !g.frozen eq 0 then show end ;+ ; This procedure calibrates a single integration from a Nod scan pair. ; This is currently an internal function used only by getnod ; ; @param s1_b1_off {in}{required}{type=spectrum} An uncalibrated ; spectrum from feed 0 (source/tracking feed) in scan 1 with the cal off. ; @param s1_b1_on {in}{required}{type=spectrum} An uncalibrated ; spectrum from feed 0 (source/tracking feed) in scan 1 with the cal on. ; @param s2_b1_off {in}{required}{type=spectrum} An uncalibrated ; spectrum from feed 0 (reference/non-tracking feed) in scan 2 with the cal off. ; @param s2_b1_on {in}{required}{type=spectrum} An uncalibrated ; spectrum from feed 0 (reference/non-tracking feed) in scan 2 with the cal on. ; @param s1_b2_off {in}{required}{type=spectrum} An uncalibrated ; spectrum from feed 1 (reference/non-tracking feed) in scan 1 with the cal off. ; @param s1_b2_on {in}{required}{type=spectrum} An uncalibrated ; spectrum from feed 1 (reference/non-tracking feed) in scan 1 with the cal on. ; @param s2_b2_off {in}{required}{type=spectrum} An uncalibrated ; spectrum from feed 1 (source/tracking feed) in scan 2 with the cal off. ; @param s2_b2_on {in}{required}{type=spectrum} An uncalibrated ; spectrum from feed 1 (source/tracking feed) in scan 2 with the cal on. ; @param tau {in}{required}{type=float} tau at zenith ; @param tsys {in}{required}{type=float} tsys at zenith ; @param ap_eff {in}{required}{type=float} aperture efficiency ; @param smthoff {in}{required}{type=integer} smooth factor for reference spectrum ; @param units {in}{required}{type=string} takes the value 'Jy', 'Ta', or 'Ta*' ; @param ret_tsys {out}{optional}{type=float} internal use ; @param ret_tau {out}{optional}{type=float} internal use ; ; @private ; ; @version $Id: getnod.pro,v 1.16 2005/05/30 04:03:43 bgarwood Exp $ ;- pro donodint,s1_b1_off,s1_b1_on,s2_b1_off,s2_b1_on,s1_b2_off,s1_b2_on,s2_b2_off,s2_b2_on,tau,tsys,ap_eff,smthoff,units,ret_tsys,ret_tau compile_opt idl2 elevation = s1_b1_off.elevation * !pi/180.0 ; Use the inner 80% of data to calculate mean Tsys pct10 = n_elements(*s1_b1_off.data_ptr)/10 pct90 = n_elements(*s1_b1_off.data_ptr)-pct10 mean_Tsys_b1_s1 = mean((*s1_b1_off.data_ptr)[pct10:pct90]) / $ mean(((*s1_b1_on.data_ptr) - (*s1_b1_off.data_ptr))[pct10:pct90]) * s1_b1_off.mean_tcal mean_Tsys_b1_s2 = mean((*s2_b1_off.data_ptr)[pct10:pct90]) / $ mean(((*s2_b1_on.data_ptr) - (*s2_b1_off.data_ptr))[pct10:pct90]) * s2_b1_off.mean_tcal mean_Tsys_b2_s1 = mean((*s1_b2_off.data_ptr)[pct10:pct90]) / $ mean(((*s1_b2_on.data_ptr) - (*s1_b2_off.data_ptr))[pct10:pct90]) * s1_b2_off.mean_tcal mean_Tsys_b2_s2 = mean((*s2_b2_off.data_ptr)[pct10:pct90]) / $ mean(((*s2_b2_on.data_ptr) - (*s2_b2_off.data_ptr))[pct10:pct90]) * s2_b2_off.mean_tcal sigs1 = (*s1_b1_off.data_ptr + *s1_b1_on.data_ptr)/2.0 refs1 = (*s1_b2_off.data_ptr + *s1_b2_on.data_ptr)/2.0 refs2 = (*s2_b1_off.data_ptr + *s2_b1_on.data_ptr)/2.0 sigs2 = (*s2_b2_off.data_ptr + *s2_b2_on.data_ptr)/2.0 ; Apply smthoff if requested if smthoff gt 1 then begin refs1 = smooth(refs1,smthoff) refs2 = smooth(refs2,smthoff) end ; if tsys is set by user, replace the calculated value if tsys ne 1 then begin mean_Tsys_b1_s1 = tsys*exp(tau/sin(elevation)) mean_Tsys_b1_s2 = tsys*exp(tau/sin(elevation)) mean_Tsys_b2_s1 = tsys*exp(tau/sin(elevation)) mean_Tsys_b2_s2 = tsys*exp(tau/sin(elevation)) endif calb1 = (sigs1 - refs2)/refs2 * mean_Tsys_b1_s1 calb2 = (sigs2 - refs1)/refs1 * mean_Tsys_b2_s2 ; the following should be done with the accumulate functions when they are ; available *!g.s[0].data_ptr = (calb1 + calb2)/2.0 !g.s[0].tsys = (mean_Tsys_b1_s1 + mean_Tsys_b1_s2 + mean_Tsys_b2_s1 + mean_Tsys_b2_s2)/4.0 if units eq 'Jy' or units eq 'Ta*' then begin ; apply opacity and spillover here correction = exp(tau/sin(elevation))*.99 *!g.s[0].data_ptr *= correction if units eq 'Jy' then begin *!g.s[0].data_ptr /= (2.85 * ap_eff) endif endif !g.s[0].units = units ret_tsys = [mean_Tsys_b1_s1,mean_Tsys_b2_s1,mean_Tsys_b1_s2,mean_Tsys_b2_s2] ret_tau = tau/sin(elevation) end