;+
; NAME:
;     PLOT_NL_TEST
; PURPOSE:
;     Routine to plot 5 DB attenuation factor for various values of NL
;     (nonlinearity parameter), assuming nominal offset.
; CATEGORY:
;     OVRO APC DIAGNOSTIC TEST
; CALLING SEQUENCE:
;     plot_nl_test,gcdata
; INPUTS:
;     gcdata   the GCAL data array, of size (2,7,5,86), where
;                the first index is ND off/on, second index is the
;                attenuation state (0,5,10,15,20,35,20*DB), the
;                third index is the total power channel, and the
;                last index is the frequency
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
; ROUTINES CALLED:
; OUTPUTS:
; COMMENTS:
; SIDE EFFECTS:
;     A plot is produced
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 08-Mar-1998 by Dale E. Gary
;-

pro plot_nl_test,gcdata,gcflag,ant,os

   if (n_elements(os) eq 0) then os = 0.0

   ; Make relevant arrays of only 0 DB and 5 DB data
   g1 = reform(gcdata(*,0,ant-1,*))
   g2 = reform(gcdata(*,1,ant-1,*))
   g1f = reform(gcflag(*,0,ant-1,*))
   g2f = reform(gcflag(*,1,ant-1,*))

   ; Loop over various nonlinearity values (-0.6 to +0.6, step 0.2)
   nl = -0.8
   window,xsiz=400,ysiz=500
   plot,/nodata,[0,1],[0,1],xran=[0,172],yran=[0.8,1.2],xsty=1,ysty=1
   nd = fltarr(7)
   nlin = fltarr(7)
   resid = fltarr(7)
   for i = 0, 6 do begin
      nl = nl + 0.2

      ; Calculate curve and plot it
      curve=(g1-os)*(1+nl*(g1-os)/2047.)/((g2-os)*(1+nl*(g2-os)/2047.)*10^0.5)
      curflg = g1f
      oplot,curve

      ; Zero any values where either channel are flagged bad
      bad1 = where(g1f eq 0)
      bad2 = where(g2f eq 0)
      if (bad1(0) ne -1) then curve(bad1) = 0
      if (bad2(0) ne -1) then begin
         curve(bad2) = 0
         curflg(bad2) = 0
      endif

      ; The average of the curve is the noise diode factor
      nd(i) = total(curve)/total(curflg)
      dif = curve - nd(i)
      bad = where(curve eq 0)
      if (bad(0) ne -1) then dif(bad) = 0
      ; Determine the rms variation about this average and print it.
      resid(i) = (total(dif^2)/total(curflg))^0.5
      nlin(i) = nl
      print,'ND factor, Nonlinearity and residual: ', nd(i), nlin(i), resid(i)
   endfor

   junk = min(resid,imin)
   print,'Best fit: ', nd(imin), nlin(imin)
return
end