; The only required parameter is the scan number. This can be either scan ; in the sequence of two total power scans, and the paired scan is determined ; from the header. Arguments to identify the IF number, ; polarization number and feed number are optional. ; The program 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*' ; @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: getps.pro,v 1.20 2005/05/30 04:03:43 bgarwood Exp $ ;- pro getps,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 that 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 ; check parameters info = scan_info(scan) if info.procedure ne 'OffOn' and info.procedure ne 'OnOff' 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 end ; the rest of the info should be okay end 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 fdnum lt 0 or fdnum gt (info.n_feeds-1) then begin message,"Invalid feed: " + strcompress(string(fdnum),/remove_all) + $ ". This scan has " + strcompress(string(info.n_feeds),/remove_all) + " feeds, zero-indexed.",/info return endif if plnum lt 0 or plnum gt (info.n_polarizations-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 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] ; 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,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,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_off = where(data.scan_number eq scan and data.cal_state eq 0, count1off) s1_on = where(data.scan_number eq scan and data.cal_state eq 1, count1on) s2_off = where(data.scan_number eq (scan+1) and data.cal_state eq 0, count2off) s2_on = where(data.scan_number eq (scan+1) and data.cal_state eq 1, count2on) if (count1off ne expectedCount or count1off ne count1on or count1off ne count2off or count1off ne count2on) 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] as template to hold the result old_frozen = !g.frozen freeze set_data_container,data[0] if n_elements(ap_eff) eq 0 then ap_eff = get_ap_eff(!g.s[0].observed_frequency/1.0e9,!g.s[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(!g.s[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 dopsint,data[s1_off],data[s1_on],data[s2_off],data[s2_on],tau,tsys,ap_eff,smthoff,units,ret_tsys,ret_tau endif else begin thisaccum = {accum_struct} for n_int = 0,(info.n_integrations-1) do begin dopsint,data[s1_off[n_int]],data[s1_on[n_int]],data[s2_off[n_int]],data[s2_on[n_int]], $ 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 if units eq 'Jy' then $ print,scan,ret_tau,ap_eff,ret_tsys[0],ret_tsys[1],$ format='("Scan: ",i5," units: Jy tau: ",f5.3," ap_eff: ",f5.3," Tsys: ",f7.2, 2x, f7.2)' $ else if units eq 'Ta*' then $ print,scan,ret_tau,ret_tsys[0],ret_tsys[1],$ format='("Scan: ",i5," units: Ta* (K) tau: ",f5.3," Tsys: ",f7.2, 2x, f7.2)' $ else print,scan,ret_tsys[0],ret_tsys[1],format='("Scan: ",i5," units: Ta (K) Tsys: ",f7.2, 2x, 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 position switched scan pair. ; This is currently an internal function used only by getps. ; ; @param s1_off {in}{required}{type=spectrum} An uncalibrated ; spectrum from scan 1 with the cal off. ; @param s1_on {in}{required}{type=spectrum} An uncalibrated ; spectrum from scan 1 with the cal on. ; @param s2_off {in}{required}{type=spectrum} An uncalibrated ; spectrum from scan 2 with the cal off. ; @param s2_on {in}{required}{type=spectrum} An uncalibrated ; spectrum from 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: getps.pro,v 1.20 2005/05/30 04:03:43 bgarwood Exp $ ;- pro dopsint,s1_off,s1_on,s2_off,s2_on,tau,tsys,ap_eff,smthoff,units,ret_tsys,ret_tau compile_opt idl2 elevation = s1_off.elevation * !pi/180.0 nchans = n_elements(*s1_off.data_ptr) pct10 = nchans/10 pct90 = nchans - pct10 mean_Tsys_s1 = mean((*s1_off.data_ptr)[pct10:pct90]) / $ mean((*s1_on.data_ptr)[pct10:pct90] - (*s1_off.data_ptr)[pct10:pct90])*s1_off.mean_tcal mean_Tsys_s2 = mean((*s2_off.data_ptr)[pct10:pct90]) / $ mean((*s2_on.data_ptr)[pct10:pct90] - (*s2_off.data_ptr)[pct10:pct90])*s2_off.mean_tcal s1 = (*s1_off.data_ptr + *s1_on.data_ptr)/2.0 s2 = (*s2_off.data_ptr + *s2_on.data_ptr)/2.0 ; if tsys is set by user, replace the calculated value if tsys ne 1 then begin mean_Tsys_s1 = tsys*exp(tau/sin(elevation)) mean_Tsys_s2 = tsys*exp(tau/sin(elevation)) endif if s1_off.procedure eq 'OnOff' then begin if smthoff gt 1 then s2 = smooth(s2,smthoff) !g.s[0].tsys = mean_Tsys_s2 caldata = (s1 - s2)/s2 * mean_Tsys_s2 endif else begin if smthoff gt 1 then s1 = smooth(s1,smthoff) !g.s[0].tsys = mean_Tsys_s1 caldata = (s2 - s1)/s1 * mean_Tsys_s1 endelse *!g.s[0].data_ptr = caldata if units eq 'Jy' or units eq 'Ta*' then begin ; apply opacity and spillover here correction = exp(tau/sin(elevation))*.99 ; this much is common to Ta* and Jy *!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_s1,mean_Tsys_s2] ret_tau = tau/sin(elevation) end