;+
; NAME:
;     TPSOLVE2
; PURPOSE:
;     Routine to take the on- and off-source total power calibrations,
;     and solve for the total power calibration factors.
; CATEGORY:
;     OVRO APC TP CALIBRATION
; CALLING SEQUENCE:
;     tpsolve2,tpon,onrms,onnd,onvar,tpoff,offrms,offnd,offvar,f,tls,tpfac,tprms[,/debug]
; INPUTS:
;     tpon       The on-source  total power data, of size
;                   TPON(2,2,85) => R/L,nant,nharm
;     tpoff      The off-source total power data, of size
;                   TPOFF(2,2,85) => R/L,nant,nharm
;     onrms      The on-source 1-sigma deviation in total power, of
;                   same size as TPON
;     offrms     The off-source 1-sigma deviation in total power, of
;                   same size as TPOFF
;     onnd       The on-source noise diode data, of same size as TPON
;     offnd      The off-source noise diode data, of same size as TPOFF
;     onvar      The on-source 1-sigma deviation in total power, of
;                   same size as TPON
;     offrms     The off-source 1-sigma deviation in total power, of
;                   same size as TPOFF
; OPTIONAL (KEYWORD) INPUT PARAMETERS
; ROUTINES CALLED:
;     calspec, tpinterp
; OUTPUTS:
;     tpfac      The final array of total power calibration factors,
;                   of same size/meaning as input arrays (e.g. TPON)
;     tprms      The total power 1 sigma uncertainties--corresponding to TPFAC
; COMMENTS:
; SIDE EFFECTS:
;     A total power calibration file (.TPC file) is written)
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 30-Aug-2000 by Dale E. Gary   (adapted from TPSOLVE)
;     12-Nov-2000  DG
;       Significant changes--moved MAKE_TPCAL and the plotting to TPCALCHEK2.
;-

pro tpsolve2,tpon,onrms,onnd,onvar,tpoff,offrms,offnd,offvar,$
               f,tls,oldtp,tpfac,tprms,debug=debug

   ; Maximum offsource difference (R-L) + (L-R), which is zero except
   ; for measurement errors and noise.
   MAXDIF = 10

   ; Determine the epoch, in decimal years
   epoch = tls.year + tls.day/365.

   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(onnd) eq 0,nbad)
   if (nbad gt 0) then onnd[bad] = -99.00

   indx = indgen(85)

   ; 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,onnd,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

return
end