;+
; NAME:
;     NLFIT
; PURPOSE:
;     GCALCHEK support routine.  Finds the value of nonlinearity
;     parameter (NL) that gives the least variation in the ratio
;     of 0 DB data to 5 DB data.
; CATEGORY:
;     OVRO APC GCALCHEK SUPPORT ROUTINE
; CALLING SEQUENCE:
;     nlfit,gcdata,gcflag,offsets,res[,/debug]
; INPUTS:
;     gcdata   the GCAL data array, of size (2,7,N,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 for N
;                antennas, and the last index is the frequency
;     gcflag   a parallel array to GCDATA, containing 1 where
;                data in GCDATA are good, and 0 where no good.
;     offsets  an N-element array containing the Total Power offsets,
;                where N is the number of antennas.
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
;     debug    a switch to enable plotting of some intermediate
;                results.
; ROUTINES CALLED:
; OUTPUTS:
;     res      an N-element array containing the residuals for
;                the best fit NL parameter for each antenna.
; COMMENTS:
;     The routine just does a simple loop over NL parameter from
;     -1.0 to 1.0, with step of 0.02.  Could be sped up by using
;     a minimum searching algorithm.  Precision could also nominally
;     be improved but would be misplaced--accuracy is set by
;     systematics and 0.02 is sufficient accuracy.
; SIDE EFFECTS:
;     When the /DEBUG switch is set, a plot is produced, and at
;     STOP is encountered.  Type .C or .CONT to continue
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 08-Mar-1998 by Dale E. Gary
;     09-Mar-1998  DG
;       Finalized, with some alterations. Added /DEBUG switch.
;-

function nlfit,gcdata,gcflag,offsets,res,debug=debug

   ; Determine number of antennas from data array, so that this routine
   ; will be independent of number of antennas.
   nant = (size(gcdata))(3)

   ; Make copies for alteration
   gcd = gcdata
   gcdflg = gcflag

   ; Subtract offsets and flag data where inconsistent
   for i = 0, nant-1 do begin
      gcd(*,*,i,*) = gcd(*,*,i,*) - offsets(i)
   endfor
;   minval = 0.001

   ; Use only values greater than 5 for determining ratios.  This avoids
   ; excessive digitization errors.
   minval = 5
   bad = where(gcd LT minval)  ; Flag all values below minval
   good = where(gcd GE minval)
   IF (bad(0) NE -1) THEN BEGIN
      gcd(bad) = -1
      gcdflg(bad) = 0
   ENDIF
   IF (good(0) NE -1) THEN BEGIN
      gcdflg(good) = 1
   ENDIF

   ; Declare storage for values to be saved during loop
   nlparm = fltarr(nant)
   res = fltarr(nant)

   ; Loop over antennas
   for ant = 0, nant-1 do begin
      ; Make relevant arrays of only 0 DB and 5 DB data
      g1 = reform(gcd(*,0,ant,*))
      g2 = reform(gcd(*,1,ant,*))
      bad = where(g1 eq -1,nbad)
      if (nbad ne 0) then g1(bad) = !values.f_nan
      bad = where(g2 eq -1,nbad)
      if (nbad ne 0) then g2(bad) = !values.f_nan
      good = where(finite(g1),ngood1)
      good = where(finite(g2),ngood2)
      if (ngood1 lt 3 or ngood2 lt 3) then begin
         nlparm[ant] = 0.0
         res[ant] = 100.
      endif else begin
         ; Loop over various nonlinearity values (-1.0 to +1.0, step 0.02)
         nl = -1.02
         if (keyword_set(debug)) then begin
            device,decomposed=0
            window,xsiz=400,ysiz=500
            plot,/nodata,[0,1],[0,1],xran=[0,172],yran=[0.8,1.2],xsty=1,ysty=1
         endif
         nd = fltarr(101)
         nlin = fltarr(101)
         resid = fltarr(101)
         for i = 0, 100 do begin
            nl = nl + 0.02

            ; Calculate curve (and plot it if /debug)
            curve=g1*(1+nl*g1/2047.)/(g2*(1+nl*g2/2047.)*10^0.5)
            if(keyword_set(debug)) then oplot,curve,min_val=0.5

            ; Determine average, residual, and record this non-linearity value
            avg = moment(curve,/nan)
            nd(i) = avg[0]
            resid(i) = sqrt(avg[1])
            nlin(i) = nl
         endfor

         ; Find the non-linearity and residual corresponding to minimum residual
         junk = min(resid,imin)
         nlparm(ant) = nlin(imin)
         res(ant) = resid(imin)
      endelse
      ; If /debug, overplot the result
      if (keyword_set(debug)) then begin
         print,'Best fit: ',1./nd(imin),nlin(imin),resid(imin)
         nl = nlparm(ant)
         clr = 255*(256*256. + 257.)
         curve=g1*(1+nl*g1/2047.)/(nd(imin)*g2*(1+nl*g2/2047.)*10^0.5)
         oplot,curve,color=0,min_val=0.5
         g3 = reform(gcd(*,2,ant,*))
         g4 = reform(gcd(*,3,ant,*))
         curve2=g3*(1+nl*g3/2047.)/(nd(imin)*g4*(1+nl*g4/2047.)*10^0.5)
         clr = 255*(256*256. + 1.)
         oplot,curve2,color=255,min_val=0.5
         stop
      endif
   endfor
return,nlparm
end