; Spencer A. Wolfe (2-26-2015) GBT Memo ; Program to calibrate and smooth total power spectral line data by treating the ends of a map as reference positions. ; Inputs a scan range, as well as Tcal for both polarizations for each observing session. ; Calculates Ta = Tcal*{sig - }/{} ; and then corrects for atmosphere --> Ta* = Ta * exp{tau0/sin(el)}/eta ..... Tcal is found from scal.pro using calibrator data ; see http://www.gb.nrao.edu/GBT/DA/gbtidl/integration/contrib/contrib/scal.html ; ASSUMPTION: tau=0.01 ; also truncates spectra to desired frequency range and smooths ~5km/s ; also does an rms test that compares the rms of a spectral region to that of the theoretical value ; if rms > 3*(rms theo), the the data are not added to the final fits file. This takes care of most of the integrations ; with "ripply" baselines. Bad integrations are noted in a .dat file. pro HIREDUCE, scan1, scan2, tcalpol0, tcalpol1 ; call program (after compile), scan range, and Tcal for each polarization OPENW, 1, 'noscan.dat' ; write out a file of all scans not included in the final KEEP file (bad spectra) freeze ; keep the PLOTTER from updating for h = scan1, scan2 do begin gettp, h, intnum=1, ifnum=0, plnum=0 ; call scan to get the procseq # if ((!g.s[0].procseqn mod 2) eq 0) then begin ; check procseq # (to get scans toward lower right ascension) c=0 ; set limits for the reference integrations d=9 print, "Using:", c, d for e = c,d do begin ; loop through limits and polarizations for f = 0,1 do begin get, scan=h, int=e, plnum=f, ifnum=0, cal='T' ; get calon for the reference copy, 0, 1 get, scan=h, int=e, plnum=f, ifnum=0, cal='F' ; get caloff for the reference copy, 0, 2 subtract, 1, 2 ; (calon - caoloff) for the reference accum, f ; accumulate for each polarization endfor endfor ave, 0 ; for pol 0 copy, 0, 5 ave, 1 copy, 0, 6 ; for pol 1 sclear, 0 ; clear accumulators 1 and 2 sclear, 1 for f = 0, 1 do begin ; loop through polarizations select, scan=h, int='0:9', plnum=f, ifnum=0 ; get ref avgstack ; copy, 0, 7+f emptystack endfor for g = 10, 235 do begin ; loop through the signal integrations and polarizations for f = 0, 1 do begin gettp, h, intnum=g, plnum=f, ifnum=0 ; get the signal integrations (sig) copy, 0, 9+f endfor subtract, 9, 7, 11 ; sig - (pol0) subtract, 10, 8, 12 ; (pol1) divide, 11, 5 ; (sig - )/ for pol 0 scale, tcalpol0*exp(0.01/sin(!g.s[0].elevation*3.14159/180.0))/0.99 ; Ta* (pol0) nregion, [2667,2975,4219,5772] nfit, 2 ; fit a second order baseline baseline stats, 4219, 5772, /chan, /quiet, ret=mystats0 if mystats0.rms lt 0.95 then begin ; rms test x=dcextract(!g.s[0], 2667, 5772) ; extract relevent spectral range set_data_container, x boxcar, 16, /decimate ; smooth to 5km/s and remove redundant channels copy, 0, 13 keep, 13 ; keep the spectrum endif else begin print, 'disregarding scan, pl, int:', h, 0, g ; throw out spectra when rms values are too high (bad spectra) printf, 1, h, 0, g endelse data_free, x divide, 12, 6 ; repeat Ta*, smoothing and keep for pol1 scale, tcalpol1*exp(0.01/sin(!g.s[0].elevation*3.14159/180.0))/0.99 nregion, [2667,2975,4219,5772] nfit, 2 baseline stats, 4219, 5772, /chan, /quiet, ret=mystats1 if mystats1.rms lt 0.95 then begin x=dcextract(!g.s[0], 2667, 5772) set_data_container, x boxcar, 16, /decimate copy, 0, 14 keep, 14 endif else begin print, 'disregarding scan, pl, int:', h, 1, g printf, 1, h, 1, g endelse data_free, x endfor endif if ((!g.s[0].procseqn mod 2) eq 1) then begin ; same code, but now for scans toward higher right ascension c=226 d=235 print, "Using:", c, d for e = c,d do begin for f = 0,1 do begin get, scan=h, int=e, plnum=f, ifnum=0, cal='T' copy, 0, 1 get, scan=h, int=e, plnum=f, ifnum=0, cal='F' copy, 0, 2 subtract, 1, 2 accum, f endfor endfor ave, 0 copy, 0, 5 ave, 1 copy, 0, 6 sclear, 0 sclear, 1 for f = 0, 1 do begin select, scan=h, int='226:235', plnum=f, ifnum=0 avgstack copy, 0, 7+f emptystack endfor for g = 0, 225 do begin for f = 0, 1 do begin gettp, h, intnum=g, plnum=f, ifnum=0 copy, 0, 9+f endfor subtract, 9, 7, 11 subtract, 10, 8, 12 divide, 11, 5 scale, tcalpol0*exp(0.01/sin(!g.s[0].elevation*3.14159/180.0))/0.99 nregion, [2667,2975,4219,5772] nfit, 2 baseline stats, 4219, 5772, /chan, /quiet, ret=mystats0 if mystats0.rms lt 0.95 then begin x=dcextract(!g.s[0],2667,5772) set_data_container, x boxcar, 16, /decimate copy, 0, 13 keep, 13 endif else begin print, 'disregarding scan, pl, int:', h, 0, g printf, 1, h, 0, g endelse data_free, x divide, 12, 6 scale, tcalpol1*exp(0.01/sin(!g.s[0].elevation*3.14159/180.0))/0.99 nregion, [2667,2975,4219,5772] nfit, 2 baseline stats, 4219, 5772, /chan, /quiet, ret=mystats1 if mystats1.rms lt 0.95 then begin x=dcextract(!g.s[0],2667,5772) set_data_container, x boxcar, 16, /decimate copy, 0, 14 keep, 14 endif else begin print, 'disregarding scan, pl, int:', h, 1, g printf, 1, h, 0, g endelse data_free, x endfor endif endfor unfreeze ; allow plotter to update again CLOSE, 1 ; close the bad spectra file end