;Function Return_Best_Slope
;
;PURPOSE:
;   This function will return the slope that represents the best approximation
;   to a line by considering all possible combinations of the form 
;
;	m=(y_i-y_j)/(t_i-t_j)
;   
;   It does not use least squares fitting since it is specifically aimed at 
;   avoiding non-linear behavior in a CMOS ramp.  The nonlinear tail at the 
;   end of a CMOS ramp will contaiminate a LSF fit, so this function only 
;   considers pairwise fitting in order to minimize residuals.
;
;   We expect the error to have the smallest value before saturation, 
;   and in most cases, this will be during the first 1/3 of the ramp.
;
;INPUTS:
;   x - the x (independent) values
;   y - the y (dependent) values
;KEYWORDS
;   FIRSTREAD - The index of the first (x,y) pair that should be considered 
;   LASTREAD  - The index of the last (x,y) pair that should be considered
;   TVFLAG    - 0 don't display anything
;               1 plot the points and the best slope to the screen
;RETURNS:
;   OptimizedSlopes - 4 element array
;       0 - Best Slope
;       1 - Best Intercept
;       2 - Lower Read - j
;       3 - Upper Read - i
;
function Return_Best_Slope, Xg, Yg, FirstRead=FirstRead, LastRead=LastRead, $
	 TvFlag = TvFlag

If N_Elements(FirstRead) eq 0 then FirstRead = 0
If N_Elements(LastRead)  eq 0 then LastRead  = N_Elements(x)
If N_Elements(TvFlag)    eq 0 then TvFlag    = 0

;Only consider points on the way up, before dark current catches up
LastRead = (where(Yg eq max(Yg)))[0]
If LastRead lt 20 then LastRead = (LastRead - 8) > (FirstRead+2)
x = Xg(FirstRead:LastRead)
y = Yg(FirstRead:LastRead)

;Calculate number of slopes that will be found
NumReads = LastRead-FirstRead
Thirds   = NumReads/3+1
NumM     = float(NumReads*(NumReads+1))/2

;Form arrays for them
MArr   = DblArr(NumM)
BArr   = DblArr(NumM)
ErrArr = DblArr(NumM)
k      = 0.

;Keep track of the best slope
LowErr = 10000000000.

If TvFlag then begin
  !p.multi=[0,1,1]
  plot, x, y, color=fsc_color('black'), background=fsc_color('white')
endif 

;Calculate the errors
For i = FirstRead, LastRead-1 do begin
 For j=i+1, LastRead-1 do begin
   Marr(k)   = (m = (y(j)-y(i))/(x(j)-x(i)))
   BArr(k)   = (b = (y(j)-m*x(j)))
   ErrArr(k) = (Err = total(abs(y-(m*x+b))))

   If Err lt LowErr then begin
      LowErr  = Err
      BestM   = M
      BestB   = B
      MinRead = i
      MaxRead = j
   EndIf 
   k+=1.
 EndFor
EndFor

;Error Sectioning - saturated ramps are curved so that 
;more area is under the fit in the beginning reads than 
;after
Err    = y-(BestM*x+BestB)
ErrBef = Err(0:MinRead)
ErrAft = Err(MinRead+1:MaxRead)
ErrAbs = mean(abs([ErrBef,ErrAft]))
ErrDif = total(ErrAft)-total(ErrBef)
ErrRat = ErrDif/ErrAbs

if BestM lt 0 then print, 'One Best M'
 
 If TvFlag then begin
  oplot, x, BestM*x+BestB, color = fsc_color('red')
  print, 'Min Read = ', MinRead
  print, 'Max Read = ', MaxRead 
 EndIf


return, [BestM, BestB, LowErr, MinRead, MaxRead]

end