;pro Examine_Linearity
;
;PURPOSE:
;	To take an image and produce a plot showing the cutoffs for linearity.
;
;INPUTS:
;	FITSFILENAME:
;	  The File Name for which the analysis should be done
; 
;KEYWORDS:
;       ROOTDIR: 
;         The Directory holding the directories of Raw files
;       TVFLAG 
;         0 - Don't Plot anything
;         1 - Plot to the X Window
;         2 - Plot a postscript file that was opened by calling script
;         3 - Plot a png file tha can be converted for a paper 
;	MODE:
;	  0 - Only plot the slopes for bright stars
;	  1 - Plot the ramps for every pixel on the array with the exception
;	      of the pixels in the bad pixel mask.

;OUTPUTS:
;
;	  TvFlag = 2, 
;	   A postscript file will be written with all of the Fourier Plots
;	  TvFlag = 3
;	   A .png file will be written with the fourier spectrum for channel 0
;CALLING SEQUENCE:
;
;Examine_Linearity, '/nfs/slac/g/ki/ki09/lances/07Apr26/M13_Dither0_0ArcSecondsEast_0ArcSecondsSouth_H4RG_SIPIN_30_Reads_Apr27_2007_03_15_28.fits'
;Examine_Linearity,'/nfs/slac/g/ki/ki09/lances/07Apr27/Landolt110956_Dither0_0ArcSecondsEast_0ArcSecondsNorth_H4RG_SIPIN_10_Reads_Apr28_2007_02_31_54.fits'
pro Examine_Linearity, FitsFiles, RootDir = RootDir, Date=Date, $
	ReducedDir = ReducedDir, $
	ObjectName=ObjectName, NReads = NReads, $
	ReadNum = ReadNum, ThisFilter = ThisFilter,$
	OverRide = OverRide, TvFlag = TvFlag, AvgEntireFrame = AvgEntireFrame, $
	ExamineEachChannel = ExamineEachChannel, SubtractDark = SubtractDark, $
	BiasSub = BiasSub, SkipRow = SkipRow, Mode = Mode, Ch0 = Ch0

  If N_Elements(OverRide)     eq 0 then OverRide = 0
  If N_Elements(SkipRow)      eq 0 then SkipRow = 0
  If N_Elements(SubtractDark) eq 0 then SubtractDark = 1
  If N_Elements(Mode)         eq 0 then Mode = 0
  If N_Elements(Ch0)	      eq 0 then Ch0  = 1

  Common KeyParams, KeyStr
  Common ColorStr, Color, ColorNum, ColorTri, ColorNumTri

  ;Get all of the defaults for the keywords from the KeywordStruct file
  @KeywordStruct.pro
  @PlotSettings.pro

  Fits_Read, FitsFiles, NoData, FitsHeader, /header_only       
  FileParams    = Return_File_Params(FitsFiles)
  SlopeStr      = Return_Slope_Params(FitsHeader)

  ;Get the Pixel Masks for use if needed
  Masks   = Return_Pixel_Masks()
  BadMask = Masks.BadMask
  BadInd  = 0D

if KeyStr.TvFlag eq 4 then begin
   window, xsize=768, ysize=768	
   !p.multi = [0,1,3]
   !p.charsize=2
endif else begin
   !p.multi = [0,1,1]
endelse

Case Mode of
  0: Begin
    ;*****************************************************Full Frame
   ;Get the slope file for the corresponding dither
   SlopeFileName = KeyStr.ReducedDir+KeyStr.ObjectName+KeyStr.PathDelim+$
                  KeyStr.RawKey+KeyStr.ThisFilter+'_SlopeFlatFielded.fits'
   Fits_Read, SlopeFileName, SlopeIm, exten_no = 0
   SlopeIm(where(SlopeIm lt 0))=0
   SlopeIm(where(SlopeIm gt 3000))=0

   ;Find the regions that were occupied by bright stars as opposed to 
   ;very high E cosmic rays or messed up pixels
   FitRad   = 20.
   HMin     = 20.           ; minimum counts necessary to see a peak for find.pro
   FWHM     = 10.           ; FWHM cutoff for find.pro
   RoundLim = [-1.,1.]     ; roundness paramtetrs for find.pro
   SharpLim = [0.2,1.]     ; sharpness parameters for find.pro
   PhotAper = [10.]         ; aperture radius used by aper.pro
   RdNoise  = [3.]         ; readnoise used by psf.pro
   PhPerADU = [1.]         ; number of photons per adu
   PSFRad   = 20.          ; radius of the circular area in which psf will be def

   Sky, SlopeIm, SkyMode, SkySig
   Find, SlopeIm, XStar, YStar, Flux, Sharp, Roundness, HMin, $
      FWHM, RoundLim, SharpLim,/silent
   Aper, SlopeIm, XStar, YStar, Mag, ErrAp, Sky, SkyErr, 1., PhotAper, -1, $
      [-1.,3000], setskyval = SkyMode, /silent
   ;Only keep points that had a decent magnitude
   GoodPoints = where(Mag lt 90)
   XStar = XStar(GoodPoints) & YStar = YStar(GoodPoints)
   Group, XStar, YStar, 20, GroupNum
   XStar = XStar(Uniq(GroupNum)) & YStar = YStar(Uniq(GroupNum))
   NumStars = N_Elements(XStar)

   ;Use a box of boxsize for plotting linearity away from the center of a source
   BoxSize = 40
   SatRegSize = 2*BoxSize+1 & RadMax = BoxSize

    For StarNum = 0, NumStars-1 do begin
      print, 'Working on star : '+Strtrim(StarNum,2)

        If XStar(StarNum)+BoxSize gt KeyStr.Naxis1 or $
	   XStar(StarNum)-BoxSize lt 0 or $
	   YStar(StarNum)+BoxSize gt KeyStr.Naxis2 or $
	   YStar(StarNum)-BoxSize lt 0 then goto, NewStar

        XSigs  = Fltarr(NumStars) & YSigs = Fltarr(NumStars)
        
        RelXCoords = Indgen(SatRegSize,SatRegSize) mod SatRegSize
        RelXCoords = RelXCoords - Mean(RelXCoords)
        RelYCoords = reverse(Indgen(SatRegSize,SatRegSize),2) / SatRegSize
        RelYCoords = RelYCoords - Mean(RelYCoords)
        RsqCoords  = RelXCoords^2+RelYCoords^2

        AvgSlopeRatios = Fltarr(RadMax)
	;Go radially outward and capture all pixels within radius
        For Rad = 0, 1 do begin ;RadMax - 1 do begin
          print, 'At a radius of : ' +Strtrim(Rad,2)
          RelPixInRad = where(sqrt(RsqCoords) ge Rad $
                          and sqrt(RsqCoords) lt Rad+1)
          AbsPixInds  = long(RelXCoords(RelPixInRad) + long(XStar(StarNum)) + $
                       (RelYCoords(RelPixInRad) + long(YStar(StarNum)))*$
		       KeyStr.Naxis1)
          NumInAnn    = N_Elements(AbsPixInds)
	  SlopeArr    = DblArr(NumInAnn*(KeyStr.TotalReads-1))
	  ADUArr      = UIntArr(NumInAnn*(KeyStr.TotalReads-1))
          For PixNum = 0, NumInAnn - 1 do begin
            X = long(AbsPixInds(PixNum)) mod KeyStr.Naxis1
            Y = long(AbsPixInds(PixNum)) / KeyStr.Naxis2
            If Masks.BadMask(X, Y) ne 1 then begin
               ;Get the dark for subtraction and fitting purposes
               Fits_Read_DataCube, KeyStr.DarkName, Dark, $
                 XStart = X, XStop = X, YStart = Y, YStop = Y, $
                 ZStart = 0, ZStop = KeyStr.TotalReads-1
	       ;Get the Image pixel ramp
               Fits_Read_DataCube, FitsFiles, Im, $
	         Offset = KeyStr.Offset, $
                 XStart = X, XStop = X, YStart = Y, YStop = Y, $
	         ZStart = 0, ZStop = KeyStr.TotalReads-1
	       ;Try to see if exponential fit works
	       ADUMax  = 40000.
	       ValidX  = (ADUMax-Im(0:KeyStr.NReads-2))
	       ValidY  = (Im(1:KeyStr.NReads-1)-Im(0:KeyStr.NReads-2))/$
		         SlopeStr.FrameTime

	       SN   = double(N_Elements(ValidX))
	       Sx  = Total(ValidX)
	       Sy  = Total(ValidY)
	       Sxx = Total(ValidX^2)
	       Sxy = Total(ValidX*ValidY)

	       Delta  = SN*Sxx - Sx*Sx
	       m      = (SN*Sxy - Sx*Sy)/Delta
	       b      = (Sxx*Sy-Sx*Sxy)/Delta

               ValidXD  = (ADUMax-Dark(1:KeyStr.NReads-2)-Dark(0))
               ValidYD  = (Dark(2:KeyStr.NReads-1)-Dark(1:KeyStr.NReads-2))/$
			  SlopeStr.FrameTime

               SND   = double(N_Elements(ValidXD))
               SxD  = Total(ValidXD)
               SyD  = Total(ValidYD)
               SxxD = Total(ValidXD^2)
               SxyD = Total(ValidXD*ValidYD)

               DeltaD  = SND*SxxD - SxD*SxD
               mD      = (SND*SxyD - SxD*SyD)/DeltaD
               bD      = (SxxD*SyD-SxD*SxyD)/DeltaD

	       ;Remove bias and dark from image
               Deltas  = float(Im(1:KeyStr.TotalReads-1))-float(Im(0))$
                            -Dark(1:KeyStr.TotalReads-1)
               Slopes  = Deltas/(indgen(KeyStr.NReads-1)+1)
               NSlopes = float(Slopes)/float(Slopes(0))
               If (where(finite(NSlopes) eq 0))[0] eq -1 then begin
                 ;Now place them in the appropriate spot in the array
                 for SlopeNum = 0, KeyStr.NReads-2 do begin
                   Index = PixNum*(KeyStr.NReads-1)+SlopeNum
                   If NSlopes(SlopeNum) gt 0 then begin
                     SlopeArr(Index)=NSlopes(SlopeNum)
                     ADUArr(Index)=Im(SlopeNum)
                   EndIf
                 Endfor
               EndIf
            EndIf
          EndFor
	  FullRange   = Max(ADUArr)-Min(ADUArr)
	  AvgSlopeArr = DblArr(FullRange)
	  Valid       = where(ADUArr ne 0)
	  ADUArr      = ADUArr(Valid)
          SlopeArr    = SlopeArr(Valid)
	  for ADUind = 0, N_Elements(SlopeArr)-1 do begin
	  endfor
	  if rad eq 0 then stop
	EndFor
 	NewStar:
      EndFor
  End
  1: Begin
  ;**********************************************Estimate Slope With Finite Diff.
  ;Check to see if the text file with the numbers exists
  TxtFile=File_Search(KeyStr.PlotDir+KeyStr.RawKey+'Method2Linearity.txt')
  If TxtFile ne '' and OverRide eq 0 then begin 
    ReadCol, TxtFile, ADU, NumGSlopes, SumGSlopes, $
                      NumLSlopes1, SumLSlopes1, $
                      NumLSlopes2, SumLSlopes2, $
                      format= 'I,f,f,f,f,f,f'
    FullRange = 65536.
    Fits_Read, Keystr.ObjectDir+KeyStr.RawKey+'LinearLimits.fits', LinLimArr, LinHead
    Fits_Read, KeyStr.ObjectDir+KeyStr.RawKey+'BestSlopes.fits',   SlopeArr, SlHead
    Fits_Read_DataCube, FitsFiles, ImBias, zstart = 0, zstop=0
    Fits_Read_DataCube, FitsFiles, ImFin, zstart=KeyStr.NReads-1, zstop = KeyStr.Nreads-1
    MinADU = Min(ImBias)
    MaxADU = 43300
    BadInd = N_Elements(Masks.BadPixels)+N_Elements(where(SlopeArr(*,*,0) lt 0))
  EndIf Else Begin
    ;Form an array to hold FFD(0), FCD(1:NReads-2), and FBD(NReads-1)
    SlopeArr  = Fltarr(KeyStr.Naxis1,KeyStr.Naxis2, 3)
    LinLimArr = UIntarr(KeyStr.Naxis1,KeyStr.Naxis2, 4)
    ThisSlope = Fltarr(KeyStr.NReads)
    Times     = SlopeStr.ReadTimes
    NReads    = KeyStr.NReads

    NumBad      = 0.     ;Number of negative slopes
    BadInd      = 0.     ;Index for bad pixel mask during scan through array
    NumBSlope   = 0.     ;Num of pixels that have slopes that increase with RN
    NumSSlope   = 0.     ;Num of pixels that have slopes that decrease with RN
    NumBNextLeak= 0.     ;Num of pixels w/inc. slopes next to leaky pixels
    NumSNextLeak= 0.     ;Num of pixels w/dec. sloeps next to leaky pixels
    FullRange   = 65536. ;Full Range of the detector in ADU
    MinADU	= 65536. ;Find minimum during loop	
    MaxADU      = 0.     ;Find maximum during loop
    UpLinInd    = 0.     ;Upper index for best linearity - found by function
    LowLinInd   = 0.	   ;Lower index for best linearity - found by function

    ;Make Arrays to hold slope ratios
    ;Global Slopes (y_i-y_o)/(t_i-t_o)
    SumGSlopes    = DblArr(FullRange)
    NumGSlopes    = DblArr(FullRange)
    ;Local Slopes (y_i+1-y_i-1)/(t_i+1-t_i-1)
    SumLSlopes1   = DblArr(FullRange)
    NumLSlopes1   = DblArr(FullRange)
    SumLSlopes2   = DblArr(FullRange)
    NumLSlopes2   = DblArr(FullRange)

    For Row = 0D, KeyStr.Naxis2-1 do begin
      ;Read the Data in incrementally to prevent a crash
      Fits_Read_DataCube, FitsFiles, Im, FitsHeader, $
          YStart = Row, YStop = Row
      Fits_Read_DataCube, KeyStr.DarkName, Dark, DarkHeader, $
          YStart = Row, YStop = Row
      If Ch0 then begin 
         Fits_Read_DataCube, FitsFiles, ImCh0, FitsHeader, $
           YStart = Row, YStop = Row, XStart = 0, XStop=KeyStr.Ch0W-1
         MedCh0Im   = median(ImCh0,dimension=1)
         Fits_Read_DataCube, KeyStr.DarkName, DarkCh0, DarkHeader, $
           YStart = Row, YStop = Row, XStart = 0, XStop=KeyStr.Ch0W-1
         MedCh0Dark = median(DarkCh0,dimension=1)
      EndIf

      For Col=0D, KeyStr.Naxis1-1 do begin
        PixelInd = Col+KeyStr.Naxis1*Row
        If PixelInd eq Masks.BadPixels(BadInd) then begin
          BadInd +=1.
        EndIf Else Begin
          ;Median files have bias subtracted already; normal ones don't 
          If stregex(KeyStr.RawName, 'Median', /boolean) then begin
            ;Subtract dark current
            Deltas = float(Im(Col,0,1:NReads-1))-$
                     float(Dark(Col,0,1:NReads-1))
            If Ch0 then Deltas = Deltas-(MedCh0Im(1:NReads-1)-$
                                        MedCh0Dark(1:NReads-1))
             
            ;Use the Return_Best_Slope function to get some info
            OptSlo       = Return_Best_Slope(Times(1:NReads-1),Deltas)
            MostLinSlope = OptSlo(0) & MostLinInt = OptSlo(1)
            MostLinErr   = OptSlo(2) 
            MinRead      = OptSlo(3) & MaxRead    = OptSlo(4)
            LowLinInd    = Im(Col,0,MinRead+1)+Im(Col,0,0)
            UpLinInd     = Im(Col,0,MaxRead+1)+Im(Col,0,0)
            SlopeArr(Col, Row, 0)  = MostLinSlope
            SlopeArr(Col, Row, 1)  = MostLinInt
            SlopeArr(Col, Row, 2)  = MostLinErr
            LinLimArr(Col, Row, 0) = LowLinInd 
            LinLimArr(Col, Row, 1) = UpLinInd
            LinLimArr(Col, Row, 2) = MinRead
            LinLimArr(Col, Row, 3) = MaxRead
            If MostLinSlope le 0 then begin
               NumBad+=1.
               goto, SkipPix
            EndIf
            Slopes  = Deltas/(Times(1:NReads-1))
	    NSlopes = float(Slopes)/MostLinSlope
                   
	    ;Create the array with all slopes
	    ThisSlope(0) = float((Im(Col,0,1))-Dark(Col,0,1))/$
				             (Times(1)-Times(0))
            ThisSlope(1) = float((Im(Col,0,2))-Dark(Col,0,2))/$
					     (Times(2)-Times(0))
	    ThisSlope(2:NReads-2) = $
	      float((Im(Col,0,3:NReads-1)-Im(Col,0,1:NReads-3))-$
		      (Dark(Col,0,3:NReads-1)-Dark(Col,0,1:NReads-3)))/$
	              (Times(3:NReads-1)-Times(1:NReads-3))	
	    ThisSlope(NReads-1) = $
              float((Im(Col,0,NReads-1)-Im(Col,0,NReads-2))-$
		      (Dark(Col,0,NReads-1)-Dark(Col,0,NReads-2)))/$
		      (Times(NReads-1)-Times(NReads-2))

	    For SlopeNum =0, NReads-1 do begin
              ;Get the ADU Value for Global and Local
              LIndex  = Im(Col,0,SlopeNum)
              If SlopeNum gt 0 then LIndex = LIndex+Im(Col,0,0)
              If SlopeNum lt NReads-1 then $
	         GIndex  = Im(Col,0,SlopeNum+1)+Im(Col,0,0)

              ;Keep track of the minimum and maximum ADU
	      MinADU = MinADU < LIndex
              MaxADU = MaxADU > LIndex < FullRange    

	      ;Examine Global SLopes
   	      If GIndex gt 0 and GIndex lt FullRange and $
	      SlopeNum lt NReads-1 then begin
		 SumGSlopes(GIndex) = SumGSlopes(GIndex)+$
	         NSlopes(SlopeNum)
		 NumGSlopes(GIndex)+=1.
              EndIf

              ;Examine Local Slopes 
	      If LIndex gt 0 and LIndex lt FullRange then begin	
	        If ThisSlope(0) ne 0 then begin 
           	  SumLSlopes1(LIndex) = SumLSlopes1(LIndex)+$
	          ThisSlope(SlopeNum)/ThisSlope(0)
	          NumLSlopes1(LIndex) +=1.
	        EndIf
	        If MostLinSlope ne 0 then begin
	          SumLSlopes2(LIndex) = SumLSlopes2(LIndex)+$
	          ThisSlope(SlopeNum)/MostLinSlope
	          NumLSlopes2(LIndex) += 1.
	        EndIf
	      EndIf
  	    Endfor

	  EndIf Else Begin
            ;In the case that this is not a median file
          EndElse
	  If KeyStr.TvFlag eq 4 and Row gt 20 and $
                  (where(NSlopes gt 1.25))[0] eq -1 then begin $
	 	  DisSlopes = ThisSlope
		  MinSlopes = min([min(nslopes),min(DisSlopes)])
                  MaxSlopes = max([max(nslopes),max(DisSlopes)])
                  plot,nslopes,yrange=[MinSlopes,MaxSlopes],ystyle=9,$
                   xtitle= 'Read Number', ytitle= 'Slope (ADU)',$
                   position= [0.075,0.55,0.925,0.95], $
                   background=fsc_color('white'), color=fsc_color('black')
                  axis,yaxis=1,yrange=[min(slopes),max(slopes)], $
                   ytitle= 'Slope/Slope!L0!N', color = fsc_color('black')
		  oplot, DisSlopes, color=fsc_color('red')
                  plot,im(Col,0,1:KeyStr.NReads-1)-dark(Col,0,1:KeyStr.NReads-1), $
                   xtitle = 'Read Number', ytitle = 'Image - Dark - Bias', $
                   color = fsc_color('black'), $
                   position = [0.075,0.05,0.5,0.5]
                   Ramp = Im(Col,0,0:KeyStr.NReads-1)
                   Ramp(1:KeyStr.NReads-1) = Ramp(1:KeyStr.NReads-1)+Ramp(0)
                  plot,Ramp, NSlopes, $
                   color=fsc_color('black'),$
                   background=fsc_color('white'),$
                   position = [0.525,0.05,0.95,0.5]
                  stop
	   EndIf
        EndElse
        SkipPix:
      Endfor
      If Row mod 100 eq 0 and Row gt 0 then print, Row
    Endfor 

    ;Write the values to a file
    get_lun, LinFile
    openw, LinFile, KeyStr.PlotDir+KeyStr.RawKey+'Method2Linearity.txt'

    ;Get the average slope at a given ADU value
    for PixInd = 0L, FullRange - 1 do begin
     printf, LinFile, FORMAT = '(%"%I\t%f\t%f\t%f\t%f\t%f\t%f")', $
     PixInd, NumGSlopes(PixInd) , SumGSlopes(PixInd), $
		NumLSlopes1(PixInd), SumLSlopes1(PixInd),$
		NumLSlopes2(PixInd), SumLSlopes2(PixInd)
    endfor

    close, LinFile
    Free_Lun, LinFile
   EndElse

    ;Make a copy to preserve for output to text file
    NumLSlopes1Cp = NumLSlopes1
    NumLSlopes2Cp = NumLSlopes2
    NumGSlopesCp  = NumGSlopes

    ;Eliminate points with insufficient measurements
    NumLSlopes1CP(where(NumLSlopes1 lt 10))=0
    NumLSlopes2CP(where(NumLSlopes2 lt 10))=0
    NumGSlopesCP(where(NumGSlopes lt 10))  =0

    ;Get the averages
    Avg1LRat = SumLSlopes1/NumLSlopes1CP
    Avg2LRat = SumLSlopes2/NumLSlopes2CP
    AvgGRat  = SumGSlopes/NumGSlopesCP

    ;The point where a goes below 95% of orginal slope
    ADRange      = MaxADU-MinADU
    FullADU      = UIndgen(FullRange)
    GNumH95Arr   = where(AvgGRat  lt 0.95 and finite(AvgGRat) eq 1 and $
                         FullADU gt MinADU+ADRange/2,NumG95)
    GlobalHNum95 = GNumH95Arr[0]
    GNumL95Arr   = where(AvgGRat  lt 0.95 and finite(AvgGRat) eq 1 and $
                         FullADU  lt MinADU+ADRange/2,NumG95)
    GlobalLNum95 = GNumL95Arr[N_Elements(GNumL95Arr)-1]
    FNum95Arr    = where(Avg1LRat lt 0.95 and finite(Avg1LRat) eq 1 and $
                         FullADU gt MinADU+ADRange/2,NumF95)
    FirstNum95   = FNum95Arr[0]
    SNum95Arr    = where(Avg2LRat lt 0.95 and finite(Avg2LRat) eq 1 and $
                         FullADU gt MinADU+ADRange/2,NumS95)
    SecNum95     = SNum95Arr[0]
   if KeyStr.TvFlag eq 2 then $
     device, filename = KeyStr.PlotDir+KeyStr.RawKey+'Method2Linearity0.eps', $
	     bits_per_pixel=8, /encapsul, ysize=11.0, xsize=12.0,/inches        
   ADUMaxFrac = 1.0
   ADUMinFrac = 0.0
   Up95Frac   = (GlobalHNum95-MinADU)/float(MaxADU-MinADU)
   Low95Frac  = (GlobalLNum95-MinADU)/float(MaxADU-MinADU)
      !p.charsize=2.0
      !p.charthick=6
      plot, uindgen(FullRange), AvgGRat, xrange=[MinADU, MaxADU],$
       yrange=[0,1.2], xtitle = 'ADU', ytitle = 'm!Ln!N', $
       background = fsc_color('white'), color = fsc_color('black'),$
       position = [0.125,0.15,0.90,0.95]
      xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.1+!X.CRANGE(0),  $
              (!Y.CRANGE(1)-!Y.CRANGE(0))*0.46+!Y.CRANGE(0), $
              'Fraction of pixels  '+$
	      strtrim(string((4096.*4096.-BadInd)/(4096.*4096.),format='(F10.4)'),2),$
              color=fsc_color('black'), charsize=2
      xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.1+!X.CRANGE(0),  $
              (!Y.CRANGE(1)-!Y.CRANGE(0))*0.38+!Y.CRANGE(0), $
              'Maximum ADU Value', $
               color=fsc_color('black'), charsize=2
      xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.50+!X.CRANGE(0),  $
              (!Y.CRANGE(1)-!Y.CRANGE(0))*0.38+!Y.CRANGE(0), $
               strtrim(string(MaxADU,format='(I6.0)'),2), $
               color=fsc_color('black'), charsize=2
      xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.65+!X.CRANGE(0),  $
              (!Y.CRANGE(1)-!Y.CRANGE(0))*0.38+!Y.CRANGE(0), $
              strtrim(string(ADUMaxFrac,format='(D10.4)'),2), $
               color=fsc_color('black'), charsize=2
      xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.1+!X.CRANGE(0),  $
              (!Y.CRANGE(1)-!Y.CRANGE(0))*0.32+!Y.CRANGE(0), $
              '95% Linearity Max', $
               color=fsc_color('black'), charsize=2
      xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.50+!X.CRANGE(0),  $
              (!Y.CRANGE(1)-!Y.CRANGE(0))*0.32+!Y.CRANGE(0), $
               strtrim(string(GlobalHNum95,format='(I6.0)'),2), $
               color=fsc_color('black'), charsize=2
      xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.65+!X.CRANGE(0),  $
              (!Y.CRANGE(1)-!Y.CRANGE(0))*0.32+!Y.CRANGE(0), $
              strtrim(string(Up95Frac,format='(F10.4)'),2), $
               color=fsc_color('black'), charsize=2
      xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.1+!X.CRANGE(0),  $
              (!Y.CRANGE(1)-!Y.CRANGE(0))*0.26+!Y.CRANGE(0), $
              '95% Linearity Min',$
               color=fsc_color('black'), charsize=2
      xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.50+!X.CRANGE(0),  $
              (!Y.CRANGE(1)-!Y.CRANGE(0))*0.26+!Y.CRANGE(0), $
               strtrim(string(GlobalLNum95,format='(I6.0)'),2), $
               color=fsc_color('black'), charsize=2
      xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.65+!X.CRANGE(0),  $
              (!Y.CRANGE(1)-!Y.CRANGE(0))*0.26+!Y.CRANGE(0), $
              strtrim(string(Low95Frac,format='(F10.4)'),2), $
               color=fsc_color('black'), charsize=2
      xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.1+!X.CRANGE(0),  $
              (!Y.CRANGE(1)-!Y.CRANGE(0))*0.20+!Y.CRANGE(0), $
              'Minimum ADU Value', $
               color=fsc_color('black'), charsize=2
      xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.50+!X.CRANGE(0), $
              (!Y.CRANGE(1)-!Y.CRANGE(0))*0.20+!Y.CRANGE(0), $
               strtrim(string(MinADU,format='(I6.0)'),2), $
               color=fsc_color('black'), charsize=2
      xyouts, (!X.CRANGE(1)-!X.CRANGE(0))*0.65+!X.CRANGE(0), $
              (!Y.CRANGE(1)-!Y.CRANGE(0))*0.20+!Y.CRANGE(0), $
              strtrim(string(ADUMinFrac,format='(F10.4)'),2), $
               color=fsc_color('black'), charsize=2
     if KeyStr.TvFlag eq 2 then device,/close
     if KeyStr.TvFlag eq 2 then $
        device, filename = KeyStr.PlotDir+KeyStr.RawKey+'Method2Linearity1.eps', $
             bits_per_pixel=8, /encapsul

      plot, uindgen(FullRange), Avg1LRat, xrange=[MinADU, MaxADU],$
       yrange=[0,1.2], xtitle = 'ADU', ytitle = 'Slope/Slope!L0!N', $
       background = fsc_color('white'), color = fsc_color('black'),$
       position = [0.1,0.1,0.95,0.95]
      xyouts,0.2,0.375, 'Fraction of Pixels = '+$
          strtrim((4096.*4096.-BadInd)/(4096.*4096.),2), $
          color=fsc_color('black'), /normal
      xyouts,0.22,0.325, '95% Linearity at : '+ $
          strtrim(FirstNum95,2), $
          color=fsc_color('black'), /normal

      plot, uindgen(FullRange), Avg1LRat, xrange=[MinADU, MaxADU],$
       yrange=[-3,1.2], xtitle = 'ADU', ytitle = 'Slope/Slope!L0!N', $
       background = fsc_color('white'), color = fsc_color('black'),$
       position = [0.1,0.1,0.95,0.95]
      xyouts,0.2,0.375, 'Fraction of Pixels = '+$
          strtrim((4096.*4096.-BadInd)/(4096.*4096.),2), $
          color=fsc_color('black'), /normal
      xyouts,0.22,0.325, '95% Linearity at : '+ $
          strtrim(SecNum95,2), $
          color=fsc_color('black'), /normal

      plot, uindgen(FullRange), Avg2LRat, xrange=[MinADU, MaxADU],$
       yrange=[0,1.2], xtitle = 'ADU', ytitle = 'Slope/Slope!L1!N', $
       background = fsc_color('white'), color = fsc_color('black'),$
       position = [0.1,0.1,0.95,0.95]
      xyouts,0.2,0.375, 'Fraction of Pixels = '+$
          strtrim((4096.*4096.-BadInd)/(4096.*4096.),2), $
          color=fsc_color('black'), /normal
      xyouts,0.22,0.325, '95% Linearity at :'+ $
          strtrim(SecNum95,2), $
          color=fsc_color('black'), /normal

      ;Default back to the X Windowi
      if KeyStr.TvFlag eq 2 then device, /close
 
      ;Write this to a fits file
      Fits_Write, Keystr.ObjectDir+KeyStr.RawKey+'LinearLimits.fits', LinLimArr
      Fits_Write, KeyStr.ObjectDir+KeyStr.RawKey+'BestSlopes.fits', SlopeArr

    End
  EndCase

  stop
end