;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.
;
;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 ''
;
;CALLING SEQUENCE
;   slopefit, 'MedianDark_Blank_ForEachOf_30Reads.fits', $
;	       rootdir = '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr28/Darks/'
;
;REFRENCES:
;   Based on code written by Figer et. al at STSCI
;****************************************************************************
pro Slopefit, RawFileName, RootDir = RootDir, Date = Date, $
	 ReducedDir = ReducedDir, ObjectName = ObjectName, $
	 SubtractDark = SubtractDark, FilterStr = FilterStr, NReads = NReads, $
	 TvFlag = TvFlag, UpdateMask=UpdateMask, ApplyMask = ApplyMask, $
	 ThisFilter = ThisFilter, SkipFileNum=SkipFileNum

If N_Elements(UpdateMask) eq 0 then UpdateMask = 0
If N_Elements(ApplyMasks)  eq 0 then ApplyMasks = 1
If N_Elements(SkipFileNum) eq 0 then SkipFileNum = 0

;Include all of the keywords from the KeywordStruct file
@KeywordStruct.pro

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

;Enter a '' in order to do all of the files in the directory
If RawFileName eq '' then begin
   RawFileName = File_Search(ObjectName+'*'+strtrim(NReads, 2) +$
                 '_Reads'+'*.fits')
   RawFileName = Return_File_List(RawFileName, ThisFilter=ThisFilter)
   RawFileName = RawFileName(SkipFileNum:N_Elements(RawFileName)-1)
   print, 'About to calculate slopes for: ' 
   print, RawFileName
   ObjectDir = ReducedDir+ObjectName+PathDelim
   File_MkDir, ObjectDir
EndIf

;Treat Dark Frames differently
If stregex(RawFileName(0), 'Dark', /boolean) then begin 
   ObjectName  =  'Dark'
   ObjectDir   =  ReducedDir
   ObjectType  = 'Darks'
EndIf else begin
   ObjectDir = ReducedDir+ObjectName+PathDelim
   File_MkDir, ObjectDir
   ObjectType  = 'Objects'
EndElse

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

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

  ;Get the current file that's being operated on
  RawFile = RawFileName(FileNum)

  ;Get the base of the filename
  RawFileBase = (strsplit(RawFile,'.',/extract))[0]

  ;Get the unique keys in the filename for new filenames
  RawFileKey  = strmid(RawFileBase, stregex(RawFile,'Reads_')+6)

  ;Open the DataCube to be read out later
  Position = Fits_Open_DataCube(RawFile, RawFCB, $
		RawHeader, 'Read')

  ;Get axes and keywords related to the ramp sequence
  Naxis1 =     RAWFcb.axis(0)		;Number of Columns
  Naxis2 =     RAWFcb.axis(1)		;Number of Rows
  TotalReads = RAWFcb.axis(2)		;Number of Reads
  ;Get the slope parameters for this file
  SlopeStr = Return_Slope_Params(RawHeader)

  ;Check if File is a dark itself, if it is, don't subtract and 
  If stregex(RawFile, 'Dark', /boolean) then begin $
     SubtractDark =0
  EndIf

  ;The median files will have lower count rates due to dark subtraction
  If stregex(RawFile, 'Median', /boolean) then begin
        KeyStr.LowerLinLim = -KeyStr.ReadNoiseMax 
        KeyStr.UpperLinLim = 5.8e4
  EndIf

  ;If dark subtraction is requested, find the right median dark  
  If SubtractDark then begin
    DarkFileName = File_Search(DarkDir+'MedianDark*ForEachOf*'+$
	Strtrim(SlopeStr.NumGroups,2)+'Reads'+FilterStr+'.fits')
    If DarkFileName(0) eq '' then $
	DarkFileName = File_Search(DarkDir+'MedianBias*ForEachOf*'+$
        Strtrim(SlopeStr.NumGroups,2)+'Reads'+FilterStr+'.fits')
    If DarkFileName(0) eq '' then stop
    Position2 = Fits_Open_DataCube(DarkFileName, DarkFCB, DarkHeader, 'Read')

  EndIf 

  ;Open up a temp file to hold some difference values
  TempAxes = [Naxis1, Naxis2, TotalReads]
  MkHdr, TempHeader, TempAxes
  FxAddPar, TempHeader, 'BITPIX', 32
  TempDir = ReducedDir+'Temp'+PathDelim
  File_Mkdir, TempDir
  TempFileName = TempDir+RawFileKey+'ObjectMinusDarkTemp.fits'
  Position3 = Fits_Open_DataCube(TempFileName, TempFCB, TempHeader, 'Write')
   
  ;Masks and arrays for interesting things
  SatPix     = bytArr(Naxis1,Naxis2)	  ;A mask that maps saturated pixels
  DeadPix    = UIntarr(Naxis1,Naxis2)     ;A mask that maps dead pixels
  Mask10Sig  = intarr(Naxis1, Naxis2)     ;Values outside the 10 sigma range
  CosmicRays = UIntarr(Naxis1, Naxis2)    ;Values of the cosmic ray hits
  CosmicMask = bytarr(Naxis1, Naxis2)	  ;A mask to eliminate bright objects

  ;Running Sums, Cross-Terms and Values according to 15.2 in Numerical Recipes
  ;Delta = S*Sxx - Sx^2 -- a = (SxxSy - SxSxy)/Delta -- b = (SSxy - SxSy)/Delta
  ;Number of valid points is stored in SN and the indices of these points is 

  Sxy = dblarr(Naxis1, Naxis2)		
  Sx  = dblarr(Naxis1, Naxis2)
  Sy  = dblarr(Naxis1, Naxis2)
  Sxx = dblarr(Naxis1, Naxis2)
  SN  = bytarr(Naxis1, Naxis2, TotalReads)

  ;Hold onto the first read as a bias correction
  RawBiasFrame = long(Fits_Read_DataCube_Region( RawFile, RawFCB,$
		            RawHeader, $
			    0, Naxis1-1, 0, Naxis2-1, 0, 0))
 
  ;Look for Bright Objects or saturated things - use mask for saturation
  BrightInds = where(RawBiasFrame gt KeyStr.UpperLinLim)
  OutRangePix = bytarr(Naxis1,Naxis2)

  ;Get bias frame for dark exposure -- All Zeroes
  If SubtractDark then begin
    DarkBiasFrame = long(Fits_Read_DataCube_Region(DarkFileName, DarkFCB, $
		 	DarkHeader,0, Naxis1-1, 0, Naxis2-1, 0, 0))
  EndIf 

  ;Make an array that will hold the previous frame for loop subtraction
  PreviousFrame = fltarr(Naxis1, Naxis2)

  ;///////////////////////////////////////////////////////BEGIN LOOP OVER READS
  For ReadNum = 0, TotalReads-1 do begin

    print, 'Working on Read ' + Strtrim(ReadNum,2)

    ;Get the time of this read in the exposure
    Time = float(SlopeStr.GroupTimes(ReadNum))
 
    ;Get the current Raw frame
    RawFrame = long(Fits_Read_DataCube_Region( RawFile, RawFCB, $
			RawHeader, $
			0, Naxis1-1, 0, Naxis2-1, ReadNum, ReadNum))

    ;Find the saturated pixels and stop calculating slopes for these
    NewLowPix = where(RawFrame lt KeyStr.LowerLinLim, NumNewLow)
    NewSatPix = where(RawFrame gt KeyStr.UpperLinLim, NumNewSat)
    If NumNewSat ne 0 then OutRangePix(NewSatPix) = 1
    If NumNewLow ne 0 then OutRangePix(NewLowPix) = 1

    ;Subtract the bias
    Frame = RawFrame-RawBiasFrame

    ;Increment the dead pixel mask where the difference between the 
    ;current read and the zeroth read is 0
    DeadPix(where(Frame lt ReadNoiseMax)) += 1

    If SubtractDark then begin
      ;Get the Dark Frame
      DarkFrame = long(Fits_Read_DataCube_Region(DarkFileName, DarkFCB, $
	DarkHeader, 0, Naxis1-1, 0, Naxis2-1, ReadNum, ReadNum))

      ;Subtract bias from dark (is this right?)
      DarkFrame = DarkFrame-DarkBiasFrame
      Frame = long(Frame) - long(DarkFrame)
    EndIf

    ;Append it to the temporary file
    TempWrite = Fits_Append_DataCube(TempFCB, Frame)

    ;Display histogram if keyword is set
    If TvFlag then begin

       window, 0 , xpos = 2*fix(float(XScreenSize)/3), $
                         ypos = 0, $
                         xsize = fix(float(XScreenSize)/3), $
                         ysize = fix(float(YScreenSize)/2)

       FrameHist= histogram(Frame(Ch0W:Naxis1-Ref/2-1,Ref/2:Naxis2-Ref/2-1),$
	locations=LocPix)
       plot,LocPix, FrameHist,psym=10

       window, 1 , xpos = 2*fix(float(XScreenSize)/3), $
                         ypos = fix(float(YScreenSize)/2), $
                         xsize = fix(float(XScreenSize)/3), $
                         ysize = fix(float(YScreenSize)/2)
        FrameIm = bytscl(Frame, $
		  min = min(Frame(Ch0W:Naxis1-Ref/2-1,Ref/2:Naxis2-Ref/2-1)),$
		  max = max(Frame(Ch0W:Naxis1-Ref/2-1,Ref/2:Naxis2-Ref/2-1)))
	tvscl,congrid(FrameIm,512,512)

	window, 2, xpos = 0, ypos = 0, $
	        xsize = fix(float(XScreenSize)/3), $
                ysize = fix(float(YScreenSize)/2)

	NegSlopeHist = histogram(Frame(Masks.HDCInd), locations = LocHDC)
        plot, LocHDC, NegSlopeHist, psym=10
	stop
    EndIf

    ;Get difference between this and last frame
    DiffFrame = Frame - PreviousFrame
    CRSDev   = stddev(DiffFrame[Ch0W:Naxis1-Ref/2-1,Ref/2:Naxis2-Ref/2-1])
    CRMed    = median(DiffFrame[Ch0W:Naxis1-Ref/2-1,Ref/2:Naxis2-Ref/2-1])
    Outliers = where (DiffFrame gt (CRMed + CRSDev*10), NumOutliers)
    If (Outliers[0] ne -1) then begin
        Mask10Sig[Outliers] = Mask10sig[Outliers] + 1
	CosmicMask[Outliers] += 1
        CosmicRays[Outliers] = DiffFrame[Outliers]
    Endif

    ;Stop calculating slopes for the Pixels struck by comsics
    If NumOutliers ne 0 then OutRangePix(Outliers) = 1
    NotSatPix = where(OutRangePix eq 0, NumNotSatPix)

    ;Print some values
    ValidPixVals = Frame(Ch0W:Naxis1-Ref/2-1, Ref/2: Naxis2-Ref/2-1)
    print, 'Number of Pixels in linear regime : ' + $
           strtrim(NumNotSatPix,2)
    print, 'Number of Pixels below zero after bias sub: '+ $
           strtrim(N_elements(where(ValidPixVals lt 0)),2)

    ;Increment the running sums
    Sx(NotSatPix)  = Sx(NotSatPix)  + Time
    Sy(NotSatPix)  = Sy(NotSatPix)  + Frame(NotSatPix)
    Sxy(NotSatPix) = Sxy(NotSatPix) + Time * Frame(NotSatPix)
    Sxx(NotSatPix) = Sxx(NotSatPix) + Time^2
    SNThisFrame    = BytArr(Naxis1, Naxis2)
    SNThisFrame(NotSatPix) = 1
    SN(*,*,ReadNum) = SNThisFrame

   ;Update previous frame
   PreviousFrame = Frame 
   If NumNewLow ne 0 then OutRangePix(NewLowPix)=0

  EndFor
  ;////////////////////////////////////////////////////////////////////////

  ;Close the temporary file that holds the differences 
  TempClose = Fits_Close_DataCube(TempFCB)
  delvarx, TempFCB

  ;Calculate slope according to p. 662 in Numerical Recipes and prepare errors
  Delta = SlopeStr.NumGroups*Sxx-Sx*Sx
  SlopeArray     = (SlopeStr.NumGroups * Sxy - Sx * Sy) / Delta
  InterceptArray = (Sxx*Sy-Sx*Sxy)/ Delta
  ErrorArray	 = DblArr(Naxis1,Naxis2)
 
  ;Get some headers from the file that was just slope fitted
  RA     = float(SxPar(RawHeader, 'RA'))
  DEC    = float(SxPar(RawHeader, 'DEC'))
  FILTER = SxPar(RawHeader, 'FILTER')

  ;==========================================================ERROR CALCULATION
  ;Go through the file once again (necessary to avoid memory crashes) to 
  ;calculate errors in the pixels (RawBiasFrame and DarkBiasFrame still in mem)
  Position3 = Fits_Open_DataCube(TempFileName, TempFCB, TempHeader, 'Read')

  For FrameHalf = 0, 1 do begin
    print, 'Computing Error for frame half : '+strtrim(FrameHalf,2)  
    ;Read in Half the image
    RawMinusDarkHalf = Fits_Read_DataCube_Region( TempFileName, TempFCB, $
			TempHeader, $ 
			0, Naxis1-1, $
			FrameHalf*Naxis2/2, (FrameHalf+1)*Naxis2/2,$
			0, TotalReads-1)
    print, 'Got Half of image'
    ;Get the error for each pixel with respect to the slope that was calculated
    For Row = FrameHalf*Naxis2/2, (FrameHalf+1)*Naxis2/2 - 1 do begin
      RealRow = Row-FrameHalf*Naxis2/2
      For Col = Ch0W, Naxis1-1 do begin 
        If finite(SlopeArray(Col, RealRow)) then begin
          DataPoints = SN(Col,Row,*)
	  ValidDataPoints = Where(DataPoints eq 1,NumValid)
          If NumValid ge 2 then begin
     	    m = SlopeArray(Col,Row)
            b = InterceptArray(Col,Row)
            ErrorPts = RawMinusDarkHalf(Col,RealRow,ValidDataPoints) - $
		     m*SlopeStr.GroupTimes(ValidDataPoints)+b
            ErrorArray(Col,Row) = StdDev(ErrorPts,/nan)
	  EndIf Else Begin
	    ErrorArray(Col,Row) = 0
	  EndElse
        EndIf
      EndFor
    EndFor
  EndFor
  ;===========================================================================  

  ;Close the datacubes
  If SubtractDark then DarkClose=Fits_Close_DataCube(DarkFCB)
  RawClose=Fits_Close_DataCube(RawFCB)
  TempClose = Fits_Close_DataCube(TempFCB)

  ;***********************************************************************
  ;Form a header with some of the necessary keywords
  SlopeArray=[[[SlopeArray]],[[InterceptArray]],[[ErrorArray]]]
  MkHdr, SlopeRawHeader, SlopeArray
  FxAddPar, SlopeRawHeader, 'RA',  RA, 'Right Ascension (Hours)'
  FxAddPar, SlopeRawHeader, 'DEC', DEC, 'Declination (Degrees)'
  FxAddPar, SlopeRawHeader, 'FILTER', SlopeStr.Filter, 'Filter Used'
  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, 'PROFILE', 'slopefit.pro', $
	'IDL File used to create slope'
  Fits_Write, ObjectDir + RawFileBase+'_Slope.fits', SlopeArray, SlopeRawHeader

  ;***********************************************FLAT FIELD SLOPE
  ;Divide by the flat if requested
  If not stregex(ObjectName, 'Dark', /boolean) then begin
    FlatTest = File_Search(FlatDir + 'MedianFlat_'+strtrim(SlopeStr.Filter,2)+$
	'_ForEachOf_2Reads.fits')
    If FlatTest ne '' then begin
      Fits_Read, FlatDir + 'MedianFlat_'+strtrim(SlopeStr.Filter,2)+$
	'_ForEachOf_2Reads.fits', $
	FlatIm, FlatHeader

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

      If ApplyMasks then begin
         Masks = Return_Pixel_Masks()
         FlatFieldSlope(Masks.BadPixels)=0
      EndIf
      FxAddPar, SlopeRawHeader, 'FLATFILE', 'MedianFlat_'+$
	strtrim(SlopeStr.Filter,2)+'_ForEachOf_2Reads.fits', 'Flat field used for image'
      Fits_Write, ObjectDir+ RawFileBase+'_SlopeFlatFielded.fits', $
	FlatFieldSlope, SlopeRawHeader 
    EndIf
  Endif

  ;************************************************COSMIC RAYS TO FILE
  ;Write the Cosmic rays out to a file
  CosmicFile = CosmicsDir + 'Cosmics_'+RawFileKey+'.fits'
  FxAddPar, SlopeRawHeader, 'BITPIX', 16 
  FxAddPar, SlopeRawHeader, 'BZERO', 32768
  
  ;Remove objects that repeatedly had big difference value
  File_Mkdir, CosmicsDir
  BrightInds = where(CosmicMask gt 1)
  If BrightInds(0) ne -1 then CosmicRays(BrightInds) = 0
  Fits_Write, CosmicFile, CosmicRays, SlopeRawHeader

  ;************************************************DEAD PIXELS TO FILE
  ;Write the pixels that should be working, but aren't to a mask
  If UpdateMask then begin
    DeadPix(0:Ch0W-1,*)                  = 0	;Ch0
    DeadPix(Naxis1-Ref/2-1:Naxis1-1,*)   = 0    ;Right Ref
    DeadPix(*, 0:Ref/2-1) 	         = 0 	;Lower Ref
    DeadPix(*, Naxis2-1-Ref/2: Naxis2-1) = 0	;Upper Ref
    ;DeadPix(CosmicMask)                  = 0	;Saturated Pixels
    File_MkDir, DeadDir
    DeadMaskUp = Update_Pixel_Mask('Dead', DeadPix, TotalReads, DeadDir, $
				 ObjectType)
  EndIF
  ;Delete the old file control blocks
  Delvarx, RawFCB
  Delvarx, DarkFCB
  Delvarx, TempFCB

EndFor
;*****************************************************************

stop

End