;pro PlotSlope
;
;a routine to plot the slopes of pixels for diagnostic or presentation 
;purposes
pro PlotSlope, Col, Row, Raw, Dark, TypeString

       Slope = Raw - Dark
       plot, Slope, yrange=[min(Slope), max(Slope)],$
             title = 'Row : '+Strtrim(Row,2) + ', Col : '+Strtrim(Col,2)
       stop

end

;pro SlopeFit
;
;PURPOSE:
;
;   To fit a slope to each of the pixels in a datacube.  The slope of 
;   the pixel value should represent the signal seen by that pixel, 
;   whether it is darkcurrent or light.
;
;   This routine should properly handle saturated pixels in a multi-read
;   exposure with the exception of pixels that got saturated on the first read.
;NOTES:
;   This should accomodate both raw images and median dark frames if the 
;   filenames are set correctly.  That is, if the file name contains the 
;   word 'Bias' or 'Dark', there will be no attempt to subtract a dark
;   frame.
;
;KEYWORDS:
;   UPDATEMASK: 0 - Don't update the dead pixel mask with the number of times
;		    the pixel differences didn't get above the read noise
;		1 - update the dead pixel mask ''
;
;   TVFLAG:     0 - Default, Don't plot anything
;		1 - Plot to the X window
;	        2 - Plot to postscript
;		3 - Plot to a .png
;
;   TVPPT:	Tv Plot Pixel Type - Designates the type of pixel to be plotted	
;		1 - Plot
;		2 - Plot the first and second read saturated pixels 
;		3 - Plot the cosmic rays that occur somewhere in the ramp
;		4 - Plot the first read cosmics or potentially hot pixels
;CALLING SEQUENCE
;   slopefit, 'MedianDark_Blank_ForEachOf_30Reads.fits', $
;	       rootdir = '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr28/Darks/'
;
;
;REFRENCES:
;   Based on code used in CalnicA for the NICMOS Data Reduction
;****************************************************************************
pro Slopefit_1RG, RawFileName, RootDir = RootDir, Date = Date, $
	 ReducedDir = ReducedDir, ObjectName = ObjectName, $
	 SubtractDark = SubtractDark, FilterStr = FilterStr, $
	 NReads = NReads, ThisNDrops = ThisNDrops, $
	 TvFlag = TvFlag, TvPPT = TvPPT, UpdateMask=UpdateMask, $
	 ThisFilter = ThisFilter, SkipFileNum=SkipFileNum, $ 
	 Rem60Hz = Rem60Hz, Offset = Offset, UseBadMask = UseBadMask

If N_Elements(UpdateMask)  eq 0 then UpdateMask = 0
If N_Elements(SkipFileNum) eq 0 then SkipFileNum = 0
If N_Elements(TvPPT)       eq 0 then TvPPT = 0
Common KeyParams, KeyStr
;Include all of the keywords from the KeywordStruct file

@KeywordStruct.pro
@PlotSettings.pro

;Move to the directory
If stregex(RootDir, 'Reduced', /boolean) then begin
  cd, RootDir
Endif Else begin
  cd, RootDir + Date
EndElse

;If filename is blank, search for matching files in the directory
If N_Elements(RawFileName) eq 0 then begin
   RawFileName = File_Search(KeyStr.ObjectName+'*'+$
		 strtrim(KeyStr.NReads, 2) + '_Reads'+'*.fits')
   RawFileName = Return_File_List(RawFileName, ThisFilter=ThisFilter, $
		 ThisNDrops=ThisNDrops)
   RawFileName = RawFileName(SkipFileNum:N_Elements(RawFileName)-1)
   print, 'About to calculate slopes for: ' 
   print,  RawFileName
   KeyStr.ObjectDir = KeyStr.ReducedDir+KeyStr.ObjectName+KeyStr.PathDelim
   File_MkDir, KeyStr.ObjectDir
EndIf

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

;**************************************************BEGIN LOOP OVER FILES
;Batch mode for given filter or single image
For FileNum = 0, N_Elements(RawFileName)-1 do begin

  print, 'Working on File :' + RawFileName(FileNum)

  ;Get the current file that's being operated on and some keys
  KeyStr.RawName = RawFileName(FileNum)

  ;Get the Parameters necessary to fit this file
  Fits_Read, KeyStr.RawName, NoData, RawHeader, /header_only
  FileParam= Return_File_Params(RawFileName(FileNum))
  SlopeStr = Return_Slope_Params(RawHeader)
  ;Check if File is a dark itself, if it is, don't subtract and 
  If stregex(KeyStr.RawName, 'Median', /boolean) and $
     (stregex(KeyStr.RawName, 'Dark', /boolean) or $
      stregex(KeyStr.RawName, 'Bias', /boolean)) then begin 
     KeyStr.SubtractDark   =0
     KeyStr.DarkSlope      =1
     KeyStr.UseBadMask     =0
     KeyStr.SlopeSigThresh =6
     Masks.BadPixels       =-1
  EndIf

  ;Make a temporary mask that can be chopped to save time
  BadInd        = 0L
  NumNonLin     = 0L
  NReads        = KeyStr.TotalReads
  TimePoints    = SlopeStr.ReadTimes(1:NReads-1)
  LowerLinLim   = KeyStr.LowerLinLim
  UpperLinLim   = KeyStr.UpperLinLim
  FullLinLim    = UpperLinLim - LowerLinLim

  ;Make arrays for the Slopes and what have you
  Slopes     = Fltarr(KeyStr.Naxis1, KeyStr.Naxis2)
  Intercepts = Fltarr(KeyStr.Naxis1, KeyStr.Naxis2)
  Errors     = Fltarr(KeyStr.Naxis1, KeyStr.Naxis2)
  Cosmics    = Fltarr(KeyStr.Naxis1, KeyStr.Naxis2)
  SatMask    = BytArr(KeyStr.Naxis1, KeyStr.Naxis2)

  ;Attempt to remove fixed noise using channel 0
  If KeyStr.Rem60Hz then begin
        FrameFilters = Fltarr(KeyStr.Naxis2, Keystr.Naxis3)
        print, 'Subtracting fixed pattern noise from frames'
        for ReadNum =0, KeyStr.TotalReads-1 do begin
            ;Read in Half of the raw/median image
            Fits_Read_DataCube, KeyStr.RawName , RawFrame , RawHeader,$
		ZStart = ReadNum, ZStop = ReadNum
            FrameFilters(KeyStr.Ref/2:KeyStr.Naxis2-KeyStr.Ref/2-1,ReadNum) = $
		RefSub_Ch0(RawFrame, 'ReturnFilter')
        endfor
  EndIf

  ;Take in Frames one half at a time to avoid memory crash
  For FrHalf = 0, 1 do begin
    print, 'Computing Slopes for frame half : '+strtrim(FrHalf,2)  

    ;Read in Half of the raw/median image
    Fits_Read_DataCube, KeyStr.RawName , RawHalf , RawHeader, $
	XStart= 0, XStop = KeyStr.Naxis1-1, $
	YStart= FrHalf*KeyStr.Naxis2/2, $
	YStop = (FrHalf+1)*KeyStr.Naxis2/2-1, $
	Offset = Offset

    ;Zero the bad pixels now so they can be distinguished in first read
    ;saturation.  Allow for disabling for dark slopes.
    If KeyStr.UseBadMask then begin 
      HalfBadMask=BadMask(0:KeyStr.Naxis1-1, $
			FrHalf*KeyStr.Naxis2/2:(FrHalf+1)*KeyStr.Naxis2/2-1)
      HalfBadPix = where(HalfBadMask eq 1) 
      for ReadNum = 0, KeyStr.Naxis3-1 do begin
        RawHalf(HalfBadPix,ReadNum)=0
      endfor 
      delvarx, HalfBadMask & delvarx, HalfBadPix
    EndIf

    ;Read in Half of the dark image if dark subtraction is requested
    If KeyStr.SubtractDark then $
      Fits_Read_DataCube, KeyStr.DarkName, DarkHalf, DarkHeader, $
        XStart= 0, XStop = KeyStr.Naxis1-1, $
        YStart= FrHalf*KeyStr.Naxis2/2, $
        YStop = (FrHalf+1)*KeyStr.Naxis2/2-1

    ;Loop over the pixels and slope fit them individually and iteratively, 
    ;rejecting cosmic rays  
    For Row = Double(FrHalf*KeyStr.Naxis2/2), $
	      Double((FrHalf+1)*KeyStr.Naxis2/2 - 1) do begin

      ;Get the index for the half frame
      RowH = Row-FrHalf*KeyStr.Naxis2/2
      print, Row

      ;Remove the common mode from each row in each frame
      If KeyStr.Rem60Hz then begin
         For ReadNum =0 , KeyStr.TotalReads - 1 do begin
           RawHalf(*,RowH,ReadNum) = RawHalf(*,RowH,ReadNum) - $
		FrameFilters(Row,ReadNum)
         EndFor
      EndIf

      For Col = 0D, KeyStr.Naxis1-1 do begin 
	 PixelInd = Col+KeyStr.Naxis1*Row 

         If PixelInd eq Masks.BadPixels(BadInd) then begin
	    BadInd = BadInd+1
	 EndIf Else Begin
 	    If KeyStr.SubtractDark then begin 
	       DarkRamp = DarkHalf(Col,RowH,1:KeyStr.TotalReads-1)
	    EndIf Else begin
               DarkRamp = LonArr(KeyStr.TotalReads)
            EndElse

	    ;Take the ramp as the non - bias reads
	    RawRamp    = RawHalf(Col,RowH,1:NReads-1)
	    RawBias    = RawHalf(Col,RowH,0)

            If Not KeyStr.DarkSlope then begin
	      ;Find values in the linear regim
	      GoodMask = where( RawRamp gt LowerLinLim and $
			        RawRamp lt UpperLinLim, NumGood)
	      DarkFirstRead = DarkHalf(Col,RowH,1)
              DarkBias      = DarkHalf(Col,RowH,0)
	      RawFirstRead  = RawHalf(Col,RowH,1)
	      ;Add the saturated reads to the saturated mask
	      SatVals  = where( RawRamp gt UpperLinLim, NumSat) 
              If SatVals(0) ne -1 then SatMask(Col,Row) = NumSat
 	    EndIf Else Begin
              GoodMask = where( (RawRamp +RawBias) gt LowerLinLim and $
                                (RawRamp +RawBias) lt UpperLinLim, NumGood)
	      ;SlopeFit Channel 0 just for the hell of it
	      If Col lt KeyStr.Ch0W then begin
		 NumGood  = N_Elements(RawRamp)
		 GoodMask = IndGen(NumGood) 
	      EndIf
	      DarkFirstRead = 0
              DarkBias      = KeyStr.LowerADCThreshold
	      RawFirstRead  = RawHalf(Col,RowH,1)+RawHalf(Col,RowH,0)
	    EndElse

 	    ;Perform a line fit to the points in the linear region
	    If NumGood le 1 then begin
	      NumNonLin +=1
              ;If pixel is saturated, check neigbhors to see if it's at the
	      ; center of a star.  If so, use read 0 as an estimate of flux
	      ;Make sure that borders of array are not a problem
	      ColLeft = 1 & ColRight = 1 & RowDown = 1 & RowUp = 1
	      If Col  eq 0 or Col eq KeyStr.Ch0W          then ColLeft  = 0 
	      If Col  eq KeyStr.Naxis1 - KeyStr.Ref/2 - 1 then ColRight = 0
              If RowH eq KeyStr.Ref/2 or RowH eq 0	  then RowDown  = 0 
	      If RowH eq KeyStr.Naxis2/2-1		  then RowUp    = 0
	      ;Get the 3x3 (or smaller) region around
              RegCheck = RawHalf(Col-ColLeft:Col+ColRight,$
			         RowH-RowDown:RowH+RowUp,0)
 
	      ;Get the counts of the flagged pixel
	      CenFlux  = RegCheck(ColLeft,RowDown)

	      ;Set it to zero and look at the avg. value of neighbors
	      RegCheck(ColLeft,RowDown) = 0
              SurrPix = where(RegCheck ne 0, NumSurrPix)
              If NumSurrPix ne 0 then begin
                djs_iterstat, RegCheck(where(RegCheck ne 0)), $
                   mean = SurMean, sigma = SurSigma
	      Endif Else begin
                SurMean = 0 & SurSigma = 0
              EndElse 
	      ;If it's an outlier, flag it as a cosmic hit or a bad pixel
	      ;that hasn't been discovered yet
	      If CenFlux gt SurMean+2.5*SurSigma and $
		 Not KeyStr.DarkSlope then begin 
                 Slopes(Col,Row) = 0
		 print, 'Hot Pixel or Cosmic at : ', Col, Row
		 If KeyStr.TvFlag and KeyStr.TvPPT eq 4 then $
		    PlotSlope, Col, Row, RawHalf(Col,Row,*), DarkHalf(Col,Row,*), $
			       'Hot Pixel' 
              EndIf else begin
                 ;If the first read is still in the linear regime and it's
		 ;less than half the way there, use the single cds 
		 If NumGood eq 1 and $
		  RawBias lt (float(FullLinLim)/2+LowerLinLim) $
		  then begin 
	            Slopes(Col,Row)= $
			((RawFirstRead-RawBias)-$
			  DarkFirstRead)/$
			  (SlopeStr.GroupTimes(1)-SlopeStr.GroupTimes(0))
		 EndIf Else begin
		 ;If the first read is already saturated, subract 
		 ;the median dark bias and divide by a frame time
                    Slopes(Col,Row)=(RawBias-DarkBias)/$
			  SlopeStr.FrameTime
                 EndElse
		  print, 'SatPix approx slope = ', Slopes(Col,Row) 
	       EndElse
	    EndIf else begin
	       LineFit, TimePoints, RawRamp, DarkRamp, GoodMask, m, b, $
		sigm, sigb, Error, CosmicEDep, Col, Row, $
		KeyStr.SlopeSigThresh, KeyStr.TvFlag, KeyStr.TvPPT
	       Slopes(Col,Row)     = m
	       Intercepts(Col,Row) = b
	       Errors(Col,Row)     = Error 
	       Cosmics(Col,Row)    = CosmicEDep

	    EndElse
	 EndElse
      EndFor
    EndFor
EndFor
;===========================================================================  

;***************************************************************************
;SLOPEFIT OUTPUT FITS FILE *************************************************
;0 = Slope fit			;2 = Error in Fit
;1 = Y Intercept		;3 = Cosmic Rays
;Form a header with some of the necessary keywords
SlopeArray=[[[Slopes]],[[Intercepts]],[[Errors]],[[Cosmics]]]
MkHdr, SlopeRawHeader, SlopeArray
SlopeRawHeader = Update_Header(SlopeRawHeader)
FxAddPar, SlopeRawHeader, 'FRAME1', 'Slope', 'm in y=mx+b'
FxAddPar, SlopeRawHeader, 'FRAME2', 'Intercept', 'b in y=mx+b'
FxAddPar, SlopeRawHeader, 'FRAME3', 'Error', $
	'Error in fit sqrt(<(data - y=mx+b)^2>)'
FxAddPar, SlopeRawHeader, 'MINTHRESH', KeyStr.LowerLinLim, $
	'Minimum limit for linear slope calculation'
FxAddPar, SlopeRawHeader, 'MAXTHRESH', KeyStr.UpperLinLim, $
	'Maximum limit for linear slope calculation'
FxAddPar, SlopeRawHeader, 'NOISEREM', KeyStr.Rem60HzStr, $
	'Mode of noise removal for raw image'
FxAddPar, SlopeRawHeader, 'MEDDARK',  KeyStr.DarkKey, 'Name of Dark used'
FxAddPar, SlopeRawHeader, 'NOISEREM', KeyStr.Rem60HzStr, $
	'Mode of fixed noise removal'
FxAddPar, SlopeRawHeader, 'PROFILE', 'slopefit.pro', $
	'IDL File used to create slope'

;Form the filename depending on whether is was a dark slope or a lit slope
If Not KeyStr.DarkSlope then begin
   SlopeFileName = KeyStr.ObjectDir + KeyStr.RawKey+KeyStr.ThisFilter+$
        KeyStr.Rem60HzStr+'_Slope.fits'
EndIf else begin
   SlopeFileName = KeyStr.DarkSlopeDir + KeyStr.Date+'_'+$
	KeyStr.RawKey+KeyStr.ThisFilter+$
        KeyStr.Rem60HzStr+'_Slope.fits'
EndElse

Fits_Write, SlopeFileName, SlopeArray, SlopeRawHeader
print, 'Slopes written to : ' , SlopeFileName

;***********************************************FLAT FIELD SLOPE
;Divide by the flat if requested
If not KeyStr.DarkSlope then begin
  If KeyStr.FlatName ne '' then begin
    Fits_Read, KeyStr.FlatName, FlatIm, FlatHeader

    ;Normalize the Flat
    FlatIm = Double(FlatIm(*,*,1))/Max(FlatIm(*,*,1))
    FlatFieldSlope = Double(SlopeArray(*,*,0))/Double(FlatIm)

    If KeyStr.UseBadMask then begin
       FlatFieldSlope(Masks.BadPixels)=0
       SatMask(Masks.BadPixels) = 0
    EndIf
    FxAddPar, SlopeRawHeader, 'FLATFILE', KeyStr.FlatKey, $
	'Flat field used for image'
    FlatSlopeFileName = KeyStr.ObjectDir+ KeyStr.RawKey+KeyStr.ThisFilter+$
      KeyStr.Rem60HzStr+'_SlopeFlatFielded.fits'
    Fits_Write, FlatSlopeFileName, $
	FlatFieldSlope, SlopeRawHeader
    Fits_Open, FlatSlopeFileName, SlopeFCB,/append
    Fits_Write, SlopeFCB, SatMask, extname = 'SatMask'
    Fits_Close, SlopeFCB

    EndIf
Endif

;****************************************************************

EndFor
stop

End