;+ ; NAME: ; TPSOLVE ; PURPOSE: ; A temporary routine to take the two total power calibrations, ; one on-source and one off-source, and solve for the total power ; calibration factors. This routine will be replaced when a ; TRAJECTORY has been defined for total power calibrations, so that ; two calibrations are unnecessary. ; CATEGORY: ; OVRO APC TEMPORARY ; CALLING SEQUENCE: ; tpsolve ; INPUTS: ; OPTIONAL (KEYWORD) INPUT PARAMETERS ; ROUTINES CALLED: ; gaussfit ; OUTPUTS: ; COMMENTS: ; SIDE EFFECTS: ; A total power calibration file (.TPC file) is written) ; RESTRICTIONS: ; MODIFICATION HISTORY: ; Written 02-Oct-1998 by Dale E. Gary ; 04-Oct-1998 DG ; Improve results by normalizing to the noise diode (which should ; take out the frequency-to-frequency variations due to the receiver, ; standing waves, etc.), eliminating obviously bad measurements, ; fitting a smooth gaussian, interpolating missing measurements, ; then de-normalizing to obtain final factors. ; 22-Jul-1999 DG ; Changed MAXDIF from 2 to 10, due to the fact that there really is ; a systematic difference between (R-L) and (L-R) as defined here. ; Typically it can reach a value of 6 and still be perfectly good ; data. ; 16-Apr-2000 DG ; Added code to eliminate excessively small or large calibration ; factors (since they can cause problems with under or overflows ; in TPCAL factors). ;- pro tpsolve,debug=debug ; Maximum offsource difference (R-L) + (L-R), which is zero except ; for measurement errors and noise. MAXDIF = 10 file1 = pickfile(path=!defaults.workdir,filter='*.sav',$ title='Select ON SOURCE file.') IF (file1 EQ '') THEN return yrday = strmid(file1,strpos(file1,'.sav')-12,8) yr = strmid(yrday,0,4) day = strmid(yrday,5,3)/365. epoch = yr+day file2 = pickfile(path=!defaults.workdir,filter='*.sav',$ title='Select OFF SOURCE file.') IF (file2 EQ '') THEN return restore,file1 tpon = tpavg ; TPON(2,2,85) R/L,nant,nharm onrms = tprms ; ONRMS(2,2,85) R/L,nant,nharm onnd = ndavg ; ONND(2,2,85) R/L,nant,nharm onvar = ndrms ; ONRMS(2,2,85) R/L,nant,nharm restore,file2 tpoff = tpavg ; TPOFF(2,2,85) R/L,nant,nharm offrms = tprms ; OFFRMS(2,2,85) R/L,nant,nharm offnd = ndavg ; OFFND(2,2,85) R/L,nant,nharm offvar = ndrms ; OFFRMS(2,2,85) R/L,nant,nharm fghz = f/1000. ; Determine the expected flux, in Jy, of Cas A for all possible ; harmonics of 200 MHz flux = calspec('CASA',epoch) ; Eliminate harmonic 89 (85th element in array) which is not part ; of the observing sequence. flux = [flux(0:83),flux(85)] ; Get RMS from individual RMS in quadrature tprms = ((onrms^2 + offrms^2)^0.5)*0.5 ; Sizes are [2,2,85] tp = (tpfac = tpon) ; Loop over polarization for i = 0, 1 do begin ; Determine total power signal, normalizing to the same noise diode value tp(i,*,*) = tpon(i,*,*) - tpoff(i,*,*)*(onnd(i,*,*)/offnd(i,*,*)) ; Loop over antenna for j = 0,1 do begin ; Convert calibrator flux to SFU/polarization, and apply tpfac(i,j,*) = tp(i,j,*)*2e4/flux tprms(i,j,*) = tprms(i,j,*)*2e4/flux endfor endfor ; Do some quality checks ; The off source difference R-L must be equal and opposite to L-R offdiff = reform(tpoff(0,*,*)+tpoff(1,*,*)) ondiff = reform( tpon(0,*,*)+ tpon(1,*,*)) if (keyword_set(debug)) then begin window,2 plot,fghz,ondiff(0,*),psym=3,pos=[.15,.4,.95,.9],ytit='A(g!diR!n + g!diL!n)',$ tit='Ant 1 Raw Signal' oploterr,fghz,ondiff(0,*),tprms(0,0,*)*flux/2e4,psym=3 plot,fghz,offdiff(0,*),psym=1,pos=[.15,.1,.95,.4],/noerase window,3 plot,fghz,ondiff(1,*),psym=3,pos=[.15,.4,.95,.9],ytit='A(g!diR!n + g!diL!n)',$ tit='Ant 2 Raw Signal' oploterr,fghz,ondiff(1,*),tprms(0,1,*)*flux/2e4,psym=3 plot,fghz,offdiff(1,*),psym=1,pos=[.15,.1,.95,.4],/noerase endif ; Flag values that have too great a difference (R-L) - (L-R) for i = 0,1 do begin tp = tpfac(i,*,*) bad = where(abs(offdiff) gt MAXDIF,nbad) if (nbad gt 0) then tp(bad) = 0 tpfac(i,*,*) = tp endfor ; Set all non-finite values to zero (temporarily) bad = where(finite(tpfac) eq 0,nbad) if (nbad gt 0) then tpfac(bad) = 0 ; Now set all 0 or negative values to -99 bad = where(tpfac le 0,nbad) if (nbad gt 0) then begin tpfac(bad) = -99 tprms(bad) = -99 endif ; Now set all non-finite values to -99 in NDAVG bad = where(finite(ndavg) eq 0,nbad) if (nbad gt 0) then ndavg(bad) = -99.00 indx = indgen(85) titles = [['Ant 1, Lin Feed','Ant 1, LCP Feed'],$ ['Ant 2, RCP Feed','Ant 2, Lin Feed']] ; We can recover some missing calibration factors by interpolation between ; good values, using the noise diode to recover the small-scale variations. tpinterp,fghz,tpfac,tprms,ndavg,tpnew tpfac = tpnew ; Find TP factors less than 30 low = where(abs(tpfac) lt 30,nlow) if (nlow ne 0) then begin ; For those values less than 30, if S/N is less than 1, set to 30. bad = where(tpfac[low]/tprms[low] lt 1,nbad) if(nbad ne 0) then begin bad = low[bad] tpfac[bad] = 30. endif endif ; Find TP factors greater than 300, with frequency gt 12.4 GHz high = where(tpfac gt 300,nhigh) ; Keep only indexes greater than 228, which is where frequencies gt 12.4 GHz ; start. keep = where(high gt 228,nkeep) if (nkeep ne 0) then begin high = high[keep] tpfac[high] = 300. endif ; Create a TPCAL segment for inclusion into the raw data files. make_tpcal,tls,tpfac,tprms,ndavg ; Loop over feeds for i = 0, 1 do begin ; Loop over 27-m antennas for j = 0, 1 do begin ; Plot the actual calibration factors and error bars window,i*2+j plot_io,fghz,tpfac(i,j,*),yran=[5,2000],ysty=1,$ tit=titles[i,j],xtit='Freq (GHz)',psym=3 oploterr,fghz,tpfac(i,j,*),tprms(i,j,*) endfor endfor ; ; Here we normalize the calibration factors by the noise diode values, ; ; which should result in a smooth (but noisy) variation with frequency. ; ; We reduce the noise as much as possible by only allowing positive points ; ; whose variation from point to point is within 0.2 of zero to eliminate ; ; outliers. We then fit a simple gaussian to the remaining normalized ; ; values. ; ; ; Normalize the total power factors by the average noise diode value ; tpnorml = reform(tpfac(i,j,*)/ndavg(i,j,*)) ; ; Eliminate negative or high positive values. ; good = where(tpnorml gt 0 and tpnorml lt 2,ngood) ;; tpnrms = tprms(i,j,*)/ndavg(j,*) ; if (ngood gt 0) then begin ; igood = indx(good) ; tpgood = tpnorml(good) ; endif ; ; ; Take derivative and use only points with small derivative, hence ; ; points that agree well with their neighbors, in an attempt to remove ; ; outliers. ; dtpgood = tpgood - shift(tpgood,1) ; good2 = where(dtpgood gt -0.2 and dtpgood lt 0.2) ; ; ; Fit a simple gaussian shape to this second round of good points, ; ; to determine a smooth curve through them. ; fit = gaussfit(igood(good2),tpgood(good2),a,nterms=3) ; ; ; Evaluate the fit for all points, interpolating or extrapolated the bad ; ; points ; z = (indx-a(1))/a(2) ; fit = a(0)*exp(-(z^2)/2) ; oplot,fghz,fit ; ; ; Multiply the fit by the average noide diode to denormalize--these ; ; are the factors to apply to correct the total power. ; tpall = fit*reform(ndavg(i,j,*)) ; tpfac(i,j,*) = tpall ; oplot,fghz,tpfac(i,j,*),psym=3 ; endfor ; endfor ; tpnorml = tpfac(0,0,*)/ndavg(0,*) ; indx = indgen(85) ; good = where(tpnorml gt 0 and tpnorml lt 2) ; tpnrms = tprms(0,0,*)/ndavg(0,*) outfile = !defaults.dbdir+yrday+'.tpc' openw,lun,/get_lun,outfile printf,lun,'Freq ----------------------Ant 1--------------------- ----------------------Ant 2--------------------- Date ' printf,lun,'(GHz) I feed I RMS ND_inc L feed L RMS ND_inc R feed R RMS ND_inc I feed I RMS ND_inc '+yrday printf,lun,'----- -------- ------- ------ -------- ------- ------ -------- ------- ------ -------- ------- ------ ' for i = 0, 84 do begin printf,lun,fghz(i),tpfac(0,0,i),tprms(0,0,i),ndavg(0,0,i),tpfac(1,0,i),tprms(1,0,i),ndavg(1,0,i),$ tpfac(0,1,i),tprms(0,1,i),ndavg(0,1,i),tpfac(1,1,i),tprms(1,1,i),ndavg(1,1,i),$ format='(f4.1,2(2x,f7.2,1x,f7.2,1x,f7.2),2(2x,f7.2,1x,f7.2,1x,f7.2))' endfor free_lun,lun return end