; ; 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 sigscan {in}{required}{type=integer} M&C scan number for the ; signal scan ; @param refscan {in}{required}{type=integer} M&C scan number for the ; reference scan ; @keyword ifnum {in}{optional}{type=integer} IF number (starting with 0) ; @keyword intnum {in}{optional}{type=integer} integration number ; @keyword plnum {in}{optional}{type=integer} polarization number (0 or 1) ; @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: getsigref.pro,v 1.11 2005/05/30 04:03:43 bgarwood Exp $ ;- pro getsigref,sigscan,refscan,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(sigscan) eq 0 or n_elements(refscan) eq 0) then begin message, 'The scan numbers for sig and ref are 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 numbers are valid validscans = get_scan_numbers() sigCount = total(validscans eq sigscan) refCOunt = total(validscans eq refscan) if (sigCount gt 1) then message,"Warning: More than one scan with the sigscan number is in the data file.",/info if (refCount gt 1) then message,"Warning: More than one scan with the refscan number is in the data file.",/info if (sigCount lt 1) then begin message,"sigscan not available.",/info return endif if (refCount lt 1) then begin message,"refscan not available.",/info return endif infoSig = scan_info(sigscan) infoRef = scan_info(refscan) ; check parameters if infoSig.n_channels ne infoRef.n_channels then begin message,string("sig and ref scans have different numbers of channels."),/info return endif if infoSig.n_integrations ne infoRef.n_integrations then begin message,"sig and ref scans have different numbers of integrations",/info return endif if ifnum lt 0 or ifnum ge infoSig.n_ifs then begin sifnum = strcompress(string(ifnum),/remove_all) snif = strcompress(string(infoSig.n_ifs),/remove_all) message,"Illegal IF identifier: " + sifnum + ". The sigscan has " + snif + " IFs, zero-indexed.", /info return endif if ifnum lt 0 or ifnum ge infoRef.n_ifs then begin sifnum = strcompress(string(ifnum),/remove_all) snif = strcompress(string(infoRef.n_ifs),/remove_all) message,"Illegal IF identifier: " + sifnum + ". The refscan has " + snif + " IFs, zero-indexed.", /info return endif if fdnum lt 0 or fdnum gt (infoSig.n_feeds-1) then begin message,"Invalid feed: " + strcompress(string(fdnum),/remove_all) + $ ". The sigscan has " + strcompress(string(infoSig.n_feeds),/remove_all) + " feeds, zero-indexed.",/info return endif if fdnum lt 0 or fdnum gt (infoRef.n_feeds-1) then begin message,"Invalid feed: " + strcompress(string(fdnum),/remove_all) + $ ". The refscan has " + strcompress(string(infoRef.n_feeds),/remove_all) + " feeds, zero-indexed.",/info return endif if plnum lt 0 or plnum gt 1 then begin spol = strcompress(string(plnum),/remove_all) snpol = strcompress(string(infoSig.n_polarizations),/remove_all) message, "Invalid polarization identifier: " + spol + ". The sigscan has " + snpol + " polarizations, zero-indexed.", /info return endif if plnum lt 0 or plnum gt 1 then begin spol = strcompress(string(plnum),/remove_all) snpol = strcompress(string(infoRef.n_polarizations),/remove_all) message, "Invalid polarization identifier: " + spol + ". The refscan 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 infoSig.n_channels/4 then begin message, 'Invalid smthoff value', /info return endif thispol = infoSig.polarizations[plnum] thisfeed = infoSig.feeds[fdnum] ; Retrieve all the data necessary to satisfy this request singleInt = n_elements(intnum) eq 1 ; ignore any SIG states expectedCount = singleInt ? 1 : infoSig.n_integrations if singleInt then begin if intnum ge 0 and intnum le (infoSig.n_integrations-1) then begin data = !g.lineio->get_spectra(count,scan=[sigscan,refscan],ifnum=ifnum,feed=thisfeed,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=[sigscan,refscan],ifnum=ifnum,feed=thisfeed,pol=thispol) endelse if count le 0 then begin message,"No data found, this should never happen, can not continue.",/info return endif sigoff = where(data.scan_number eq sigscan and data.cal_state eq 0, countSigoff) sigon = where(data.scan_number eq sigscan and data.cal_state eq 1, countSigon) refoff = where(data.scan_number eq refscan and data.cal_state eq 0, countRefoff) refon = where(data.scan_number eq refscan and data.cal_state eq 1, countRefon) if (countSigoff ne expectedCount or countSigoff ne countSigon or countSigoff ne countRefoff or countSigoff ne countRefon) then begin message,"Unexpected number of spectra retrieved for some or all of the switching phases in one or both scans, can not continue.",/info data_free, data return endif ; copy first spectra into !g.s[0] as a 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 dosigrefint,data[sigoff],data[sigon],data[refoff],data[refon],tau,tsys,ap_eff,smthoff,units,ret_tsys,ret_tau endif else begin thisaccum = {accum_struct} for i = 0,(infoSig.n_integrations-1) do begin dosigrefint,data[sigoff[i]],data[sigon[i]],data[refoff[i]],data[refon[i]],tau,tsys,ap_eff,smthoff,units,ret_tsys,ret_tau dcaccum,thisaccum,!g.s[0],weight=weight endfor naccum = thisaccum.n accumave,thisaccum,thisave,/quiet if (naccum ne infoSig.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,sigscan,refscan,ret_tau,ap_eff,ret_tsys,$ format='("SigScan: ",i5," RefScan: ",i5," units: Jy tau: ",f5.3," ap_eff: ",f5.3," Tsys: ",f7.2)' $ else if units eq 'Ta*' then $ print,sigscan,refscan,ret_tau,ret_tsys,$ format='("SigScan: ",i5," RefScan: ",i5," units: Ta* (K) tau: ",f5.3," Tsys: ",f7.2)' $ else print,sigscan,refscan,ret_tsys,format='("SigScan: ",i5," RefScan: ",i5," units: Ta (K) Tsys: ",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 signal reference scan pair. ; This is currently an internal function used only by getsigref ; ; @param sig_off {in}{required}{type=spectrum} An uncalibrated ; spectrum from the signal scan with the cal off. ; @param sig_on {in}{required}{type=spectrum} An uncalibrated ; spectrum from signal scan with the cal on. ; @param ref_off {in}{required}{type=spectrum} An uncalibrated ; spectrum from reference scan with the cal off. ; @param ref_on {in}{required}{type=spectrum} An uncalibrated ; spectrum from reference scan 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} tsys applied; internal use ; @param ret_tau {out}{optional}{type=float} tau applied; internal use ; ; @private ; ; @version $Id: getsigref.pro,v 1.11 2005/05/30 04:03:43 bgarwood Exp $ ;- pro dosigrefint,sig_off,sig_on,ref_off,ref_on,tau,tsys,ap_eff,smthoff,units,ret_tsys,ret_tau compile_opt idl2 elevation = ref_off.elevation * !pi/180.0 Tsys_sig = *sig_off.data_ptr/(*sig_on.data_ptr - *sig_off.data_ptr)*sig_off.mean_tcal Tsys_ref = *ref_off.data_ptr/(*ref_on.data_ptr - *ref_off.data_ptr)*ref_off.mean_tcal ; Use the inner 80% of data to calculate mean Tsys pct10 = 0.1*n_elements(Tsys_sig) pct90 = n_elements(Tsys_sig)-pct10 mean_Tsys_sig = mean(Tsys_sig[pct10:pct90]) mean_Tsys_ref = mean(Tsys_ref[pct10:pct90]) avgsig = (*sig_off.data_ptr + *sig_on.data_ptr)/2.0 avgref = (*ref_off.data_ptr + *ref_on.data_ptr)/2.0 ; if tsys is set by user, replace the calculated value if tsys ne 1 then begin mean_Tsys_sig = tsys*exp(tau/sin(elevation)) mean_Tsys_ref = tsys*exp(tau/sin(elevation)) endif if smthoff gt 1 then avgref = smooth(avgref,smthoff) !g.s[0].tsys = mean_Tsys_ref caldata = (avgsig - avgref)/avgref * mean_Tsys_ref ; apply opacity and spillover here *!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_ref ret_tau = tau/sin(elevation) end