; The only required parameter is the scan number. Arguments to ; identify the IF number, polarization number, and feed number are ; optional. A zenith opacity can be entered. The ; program calculates Tsys based on the Tcal values, and cal_on/cal_off ; phases. ; ;
; 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. ; ;
; There must be 4 switching states and the SWITCH_STATE value must be 'FSWITCH'. ; ; @param scan {in}{required}{type=integer} scan number ; @keyword ifnum {in}{optional}{type=integer} IF number ; @keyword intnum {in}{optional}{type=integer} integration number ; @keyword plnum {in}{optional}{type=integer} polarization index (0 or 1) ; @keyword fdnum {in}{optional}{type=integer} feed number (default 0) ; @keyword tau {in}{optional}{type=float} tau at zenith ; @keyword ap_eff {in}{optional}{type=float} aperture efficiency ; @keyword units {in}{optional}{type=string} takes the value 'Jy', 'Ta', or 'Ta*' ; @keyword nofold {in}{optional}{type=flag} if set, does no folding. ; When this is set, data container 0 contains the results of ; (sig-ref)/ref while data container 1 contains the results of ; (ref-sig)/sig, calibrated independently and averaged over all ; integrations. Only data container 0 will be shown. ; @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). ; ; @examples ;
; getfs,76 ; accum ; getfs,77 ; accum ; ave ; show ;;
a nofold example, show each half and then fold ;
; getfs,76,/nofold ; oshow,1 ; fold ;; ; @uses dcfold ; ; @version $Id: getfs.pro,v 1.21 2005/05/30 04:03:43 bgarwood Exp $ ;- pro getfs,scan,ifnum=ifnum,intnum=intnum,plnum=plnum,fdnum=fdnum,tau=tau,$ ap_eff=ap_eff,units=units,nofold=nofold,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(units) eq 0 then units = 'Ta' thisnofold=keyword_set(nofold) if keyword_set(noweight) then weight = 1.0 ; else undefined and use default weight ; 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.n_switching_states ne 4 then begin message,string("Cannot handle this scan: number of switching states = ",info.n_switching_states), /info return end if info.n_switching_states ne 4 then begin message,string("This does not appear to be frequency switched data."),/info return 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 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,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,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 if data[0].switch_state ne 'FSWITCH' then begin message,'This is apparently not frequency switched data. switch_state = ' + data[0].switch_state,/info return endif sigwcal = where(data.cal_state eq 1 and data.sig_state eq 1, countSigwcal) sig = where(data.cal_state eq 0 and data.sig_state eq 1, countSig) refwcal = where(data.cal_state eq 1 and data.sig_state eq 0, countRefwcal) ref = where(data.cal_state eq 0 and data.sig_state eq 0, countRef) if (countSigwcal ne expectedCount or countSigwcal ne countSig or countSigwcal ne countRefwcal or countSigwcal ne countRef) 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 ; !g.s[0] gets set in call to getfsint, don't set it here old_frozen = !g.frozen freeze 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 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 return endif if singleInt then begin dofsint,data[sigwcal],data[sig],data[refwcal],data[ref],tau,ap_eff,units,thisnofold,ret_tau,$ result1=result1,result2=result2 endif else begin res1accum = {accum_struct} if thisnofold then res2accum = {accum_struct} for n_int = 0,(info.n_integrations-1) do begin dofsint,data[sigwcal[n_int]],data[sig[n_int]],data[refwcal[n_int]],data[ref[n_int]],tau,ap_eff,$ units,thisnofold,ret_tau,result1=result1,result2=result2 dcaccum,res1accum,result1,weight=weight if thisnofold then dcaccum,res2accum,result2,weight=weight end naccum1 = res1accum.n accumave,res1accum,result1,/quiet if naccum1 ne info.n_integrations then begin message,'unexpected problems in obtaining the average, can not continue',/info ; clean up if old_frozen eq 1 then freeze else unfreeze data_free,result1 if thisnofold then begin data_free,result2 accumclear, res2accum endif data_free, data return endif if thisnofold then begin naccum2 = res2accum.n accumave,res2accum,result2,/quiet if naccum2 ne naccum1 then begin message,'unexpected problems in obtaining the average of one of the results, can not continue',/info ; clean up if old_frozen eq 1 then freeze else unfreeze data_free,result1 data_free,result2 data_free, data return endif endif endelse set_data_container, result1 if thisnofold then set_data_container, result2, index=1 if units eq 'Jy' then $ print,scan,ret_tau,ap_eff,result1.tsys,$ format='("Scan: ",i5," units: Jy tau: ",f5.3," ap_eff: ",f5.3," Tsys: ",f7.2)' $ else if units eq 'Ta*' then $ print,scan,ret_tau,result1.tsys,$ format='("Scan: ",i5," units: Ta* (K) tau: ",f5.3," Tsys: ",f7.2)' $ else print,scan,result1.tsys,format='("Scan: ",i5," units: Ta (K) Tsys: ",f7.2)' data_free, data data_free, result1 if thisnofold then data_free, result2 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 FS scan. ; This is currently an internal function used only by getfs. ; ; @param sigwcal {in}{required}{type=spectrum} An uncalibrated ; spectrum from the signal phase with the cal on. ; @param sig {in}{required}{type=spectrum} An uncalibrated ; spectrum from the signal phase with the cal off. ; @param refwcal {in}{required}{type=spectrum} An uncalibrated ; spectrum from the reference phase with the cal on. ; @param sig {in}{required}{type=spectrum} An uncalibrated ; spectrum from the reference phase with the cal off. ; @param tau {in}{required}{type=float} tau at zenith ; @param ap_eff {in}{required}{type=float} aperture efficiency ; @param units {in}{required}{type=string} units for calibration ('Ta', 'Ta*', or 'Jy') ; @param nofold {in}{required}{type=string} units for calibration ('Ta', 'Ta*', or 'Jy') ; @param ret_tau {out}{optional}{type=float} returned tau; internal ; use ; @keyword result1=result1 {out}{required}{type=spectrum} The primary result. ; If nofold is not set (the default case) then this contains the ; calibrated, folded result. If nofold is set, then this contains the ; calibrated result for 1/2 of the switching phase (sig-ref)/ref. ; @keyword result2=result2 {out}{optional}{type=spectrum} This value ; is only set when nofold is set. In that case, this contains the ; calibrated result for the other 1/2 of the switching phase (ref-sig)/sig. ; ; @private ; ; @version $Id: getfs.pro,v 1.21 2005/05/30 04:03:43 bgarwood Exp $ ;- pro dofsint,sigwcal,sig,refwcal,ref,tau,ap_eff,units,nofold,ret_tau,$ result1=result1, result2=result2 compile_opt idl2 nsig = data_valid(sig) nref = data_valid(ref) nsigwcal = data_valid(sigwcal) nrefwcal = data_valid(refwcal) if (nsig le 0 or nref le 0 or nsigwcal le 0 or nrefwcal le 0) then begin message, 'One or more of the data containers are invalid' endif if (nsig ne nref or nsig ne nsigwcal or nsig ne nrefwcal) then begin message, 'One or more of the data containers does not have the same length' endif nchans = n_elements(*sig.data_ptr) pct10 = nchans/10 pct90 = nchans - pct10 tsysSig = (mean((*sigwcal.data_ptr)[pct10:pct90] + (*sig.data_ptr)[pct10:pct90]) / $ mean((*sigwcal.data_ptr)[pct10:pct90] - (*sig.data_ptr)[pct10:pct90]))/2.0*sig.mean_tcal - $ sig.mean_tcal/2.0 sigResult = (*sigwcal.data_ptr + *sig.data_ptr)/2.0 tsysRef = (mean((*refwcal.data_ptr)[pct10:pct90] + (*ref.data_ptr)[pct10:pct90]) / $ mean((*refwcal.data_ptr)[pct10:pct90] - (*ref.data_ptr)[pct10:pct90]))/2.0*ref.mean_tcal - $ ref.mean_tcal/2.0 refResult = (*refwcal.data_ptr + *ref.data_ptr)/2.0 data_copy, sig, result1 data_copy, ref, result2 *result1.data_ptr = ((sigResult - refResult) / refResult) * tsysRef *result2.data_ptr = ((refResult - sigResult) / sigResult) * tsysSig result1.tsys = tsysRef result2.tsys = tsysSig if not nofold then begin folded=dcfold(result1,result2) data_copy, folded, result1 data_free, folded data_free, result2 endif elevation = sigwcal.elevation * !pi/180.0 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 *result1.data_ptr *= correction if nofold then *result2.data_ptr *= correction if units eq 'Jy' then begin *result1.data_ptr /= (2.85 * ap_eff) if nofold then *result2.data_ptr /= (2.85 * ap_eff) endif endif result1.units = units if nofold then result2.units = units ret_tau = tau/sin(elevation) end