;pro LineFit
;
;An older version of the linefit routine used in the slopefitting
;algorithm.  Not revised yet.
pro LineFitb, x, y1, y2, Mask, m, b, sigm, sigb, Error, TvFlag

;Initialize Accumulators and Results
;Sx   = 0   & Sxx  = 0
;Sy   = 0   & Sxy  = 0
;m    = 0.0 & b    = 0.0
;Sigm = 0.0 & Sigb = 0.0
;Wt   = 1.0

NRej=1
Found=0
Color=0
ThresHold = 300
While NRej eq 1 do begin

  ;Calclate the slopes from the points 
  ValidX = double(X(Mask)) & ValidY1 = double(Y1(Mask))
  ValidY = ValidY1 - double(Y2(Mask))
  S   = double(N_Elements(Mask))
  Sx  = Total(ValidX)
  Sy  = Total(ValidY)
  Sxx = Total(ValidX^2)
  Sxy = Total(ValidX*ValidY)

  ;Plot the differences of the object frames
  RampDiff = y1(1:N_elements(y1)-1) - $
             y1(0:N_elements(y1)-2) - $
             (y2(1:N_elements(y2)-1)- $
              y2(0:N_elements(y2)-2))
  TimeDiff = x(1:N_elements(x)-1) - $
	     x(0:N_elements(x)-2)
  print, 'Diff :', RampDiff/TimeDiff
  print, 'Mean :', Mean(RampDiff/TimeDiff)
  print, 'StdDev:', StdDev(RampDiff/TimeDiff)
  print, 'NumPoints:', S

  ;Assume bright objects will have a large difference for all of 
  ;the consecutive reads in the linear A/D range and don't reject
  ;Get the slope values and the difference of each point from the fit
  Delta  = S*Sxx - Sx*Sx
  m      = (S*Sxy - Sx*Sy)/Delta
  b      = (Sxx*Sy-Sx*Sxy)/Delta
  Diff   = ValidY- (M*ValidX+B)
  Threshold = 3*StdDev(Diff)
  print, 'Diff from line:', Diff
  print, 'StdDev:', StdDev(Diff)
  ;Find Cosmic Rays where there is a negative to positive spike in the 
  ;difference
  
  DiffDiff = Diff(1:S-1) - Diff(0:S-2)
  print, DiffDiff
  print, 'Mean:', Mean(DiffDiff)
  print, 'StdDev:', StdDev(DiffDiff)
 
  PossCos  = where(DiffDiff gt Threshold)
  If PossCos(0) ne -1 then begin
     ;A cosmic ray was found, use slope before it hit
     If PossCos(0) eq 0 then begin
	oplot, validx, m*validx+b, color = 255.^(Color+1)
        NRej = 1
        Mask = Mask(1:N_elements(Mask)-1)
        If N_Elements(Mask) le 2 then begin
           m = 0 & b = 0  & error = 0
           return
        EndIf
     EndIf Else Begin
       Found = 1
       Mask = Mask(0:PossCos(0))
       NRej = 1
     EndElse
  EndIf Else begin
     ;No more cosmic rays were found
     NRej = 0
     Found=0
  EndElse
 
  oplot, validx, m*validx+b, color = fsc_color('red')
  Color+=1
  
  ;Use Error as Gauge
  Error = StdDev((ValidY - (M*ValidX + B)))

EndWhile

end

;****************************************************************************
;pro WindowStartUp
;
;Start up the windows and widgets
pro WindowStartUp

    Common ExState, State, Masks
    Common Images, DarkIm, ObjectIm, FlatIm

    State.WindowBase = widget_base(title='Window', column = 1, uvalue = 'WindowBase')
    State.MainDrawId = widget_draw(State.WindowBase, xsize=512, ysize = 512, frame = 2, $
                             uvalue = 'MainWindowGetPix', /motion_events, /button_events,$
			     keyboard_events=1)
    Widget_Control, State.WindowBase, Set_Uvalue = info, /No_Copy
    Widget_Control, /Realize, State.WindowBase
    Widget_Control, State.MainDrawId, get_value = Tmp_Value
    State.MainWindowID = Tmp_Value
    Widget_Control, State.MainDrawID, Event_Pro = 'Examine_Pixel_List_Event'
    XManager, 'Examine_Pixel_List', State.WindowBase,/no_block

End

;****************************************************************************
;pro Examine_Pixel_List_Event
;
;This routine handles the clicking of the buttons on the image
pro Examine_Pixel_List_Event, Event
 
    Common ExState
    case Event.type of 
         0: begin 
            if Event.PRESS eq 4 then begin
              device,/close
              print, 'File Closed'
            endif else begin 
              print, event.x,event.y
	      State.PixelXSub = floor(event.x * float(State.Width/512.))
              State.PixelYSub = floor(event.y * float(State.Height/512.))
              State.PixelX = State.PixelXSub + State.XStart
              State.PixelY = State.PixelYSub + State.YStart 
	      print, State.PixelX, State.PixelY
	      PlotPixels, State.PixelXSub, State.PixelYSub, $
			State.PixelX, State.PixelY, State.TvFlagLocal 
	    endelse
         end
         else: begin 
         end
    endcase 
end

;****************************************************************************
;pro PlotPixels,
;
;This routine handles the plotting on the second window
pro PlotPixels, PixelXSub, PixelYSub, PixelX, PixelY, TvFlagLocal 
    Common ExState
    Common Images
    Common KeyParams, KeyStr
    Common SlopeParams, SlopeStr, SlopeFitStr
 
    If TvFlagLocal ne 2 then !p.multi=[0,2,2]
    ;Get the true pixel value
    PixTInd = PixelX + PixelY*State.Naxis2
    ;Check to see if it falls in one of the pixel mask categories
    If (where(Masks.DeadInd eq PixTInd))[0] ne -1 then begin
       PixelType = 'Dead'
    EndIf else if (where(Masks.LeakyInd eq PixTInd))[0] ne -1 then begin
       PixelType = 'Leaky'
    Endif else if (where(Masks.LeakyNNInd eq PixTInd))[0] ne -1 then begin
       PixelType = 'Next to Leaky'
    Endif else if (where(Masks.HotInd eq PixTInd))[0] ne -1 then begin
       PixelType = 'Hot'
    Endif else begin
       PixelType = 'Good'
    EndElse

    If KeyStr.TvFlag eq 1 then begin
       resetplt
       !p.charsize=1
       !p.charthick=1
    EndIf Else If KeyStr.TvFlag eq 3 then begin 
      Set_plot,'z'
      device,set_resolution=[1400,1200]
      !p.charsize=2
      !p.charthick=2.5
      !x.thick=1.75
      !y.thick=1.75
      !p.thick=2.25
      !p.noerase=0
      loadct, 0
      !p.font = 1
      device, font = 'Helvetica Bold'
    EndIf
    
    If PixelY gt State.Ref/2 then begin
      ;Subtract Biases
      ObjectRamp=long(ObjectIm(PixelXSub, PixelYSub, 1:State.NReads-1)) 
      DarkRamp  =long(DarkIm(PixelXSub,PixelYSub,1:State.NReads-1)) 

      If TvFlagLocal ne 2 then begin
         ;Plot the Dark Image Ramp========================1
         Plot, SlopeStr.GroupTimes, DarkIm(PixelXSub, PixelYSub, *), $
	  title = 'Dark Image for '+strtrim(PixelX,2)+','+$
	  strtrim(PixelY,2), $
	  yrange=[min(DarkIm(PixelXSub,PixelYSub,*)), $
	 	 max(DarkIm(PixelXSub,PixelYSub,*))], $	
	  color = fsc_color('black'), background = fsc_color('white')

        ;Plot the Object Image Ramp=======================2
        Plot, SlopeStr.GroupTimes, ObjectIm(PixelXSub,PixelYSub,*), $
	    title = 'Object Image for '+strtrim(PixelX,2)+','+strtrim(PixelY,2), $
	    yrange= [min(ObjectIm(PixelXSub,PixelYSub,*)),$
		     max(ObjectIm(PixelXSub,PixelYSub,*))], $
	    color=fsc_color('black')
        xyouts, 0.75,0.6, 'Bias Read : ' + Strtrim(ObjectIm(PixelXSub,PixelYSub,0)), $
	    /normal, $
	    color = fsc_color('black')
        xyouts, 0.75,0.65, 'Last Read : ' + $
	    Strtrim(ObjectIm(PixelXSub,PixelYSub,State.NReads - 1)), /normal, $
	    color=fsc_color('black')
        xyouts, 0.75, 0.70, 'Pixel Type : ', /normal, color = fsc_color('black')
        xyouts, 0.85, 0.70, PixelType, /normal, color = fsc_color('red')

        ;Plot the Object - Dark===========================3
        Plot, FlatIm(PixelXSub, PixelYSub, *), $
	    title = 'Flat', $
	    color = fsc_color('black')
      EndIf
      ;Plot the Object - Dark - Bias ===================4
      Plot, SlopeStr.GroupTimes, ObjectRamp - DarkRamp, $
	    yrange=[min(ObjectRamp-DarkRamp),max(ObjectRamp-DarkRamp)],$
	    title = 'Object - Dark : No Bias '+$
	    strtrim(PixelX,2)+','+strtrim(PixelY,2), $
	    color = fsc_color('black')

      ;Get the linear fit to the object -ramp             
      Mask = where(ObjectRamp gt KeyStr.LowerLinLim and $
                   ObjectRamp lt KeyStr.UpperLinLim)
      If N_Elements(Mask) gt 2 then begin
        LineFit, SlopeStr.GroupTimes, ObjectRamp, DarkRamp, Mask, $
                m, b, sigm, sigb, Error, TvFlag
      EndIf Else begin
                m=0 & b=0 & Error = 0
      EndElse
      ;Output some numbers ====================================
      xyouts, 0.75, .20, 'Slope : ' + $
	    Strtrim(m,2),$
	    /normal, color=fsc_color('black')
      xyouts, 0.75, .15, 'Intercept :'+ $
	    Strtrim(b,2), $
	    /normal, color=fsc_color('black')
      xyouts, 0.75, .10, 'Error :'+ $
	    Strtrim(Error,2), $
	   /normal, color = fsc_color('black')

      
      ;Plot the differences of the object frames
      Diff = ObjectRamp(1:N_elements(ObjectRamp)-1) - $
	     ObjectRamp(0:N_elements(ObjectRamp)-2) - $
	     (DarkRamp(1:N_elements(DarkRamp)-1)- $
	      DarkRamp(0:N_elements(DarkRamp)-2))

      DiffDiff = Diff(1:N_elements(Diff)-1) - $
		 Diff(0:N_elements(Diff)-2)
      ;Plot, Diff, $
      ;      color = fsc_color('black')
      ;Plot, DiffDiff, $
	;    color = fsc_color('black')

     If KeyStr.TvFlag eq 3 then begin 
	   PlotName = 'PixelSlopesFor_'+Strtrim(PixelX,2)+'_'+Strtrim(PixelY,2)+'.png'
           PngImg = tvrd()
            tvlct,reds,greens,blues,/get
            write_png,KeyStr.PlotDir+PlotName, $
                 PngImg,reds,greens,blues
	   KeyStr.TvFlag = 1
           Set_plot,'x'
           resetplt 
      EndIf 

    EndIf

End

;;pro examine_pixel_list
;
;PURPOSE:
;       To take in set of pixels and plot their individual ramps for a 
;       dark and an object file.
;
;INPUTS:
;       PIXELS - A one dimensional list of pixel indices into the 2-d array 
;                (1-d coordinates). The pixel values of those coordinates will
;                be plotted  
;                If pixels is not supplied, the xstart, xstop, ystart, and ystop 
;                keywords will be used to define a region 
;KEYWORDS:
;
;       PIXELSKIP - A number of pixels that will be skipped to reach farther back into
;                   the list
;CALLING SEQUENCE
;   Examine_Pixel_List, objectfile='/nfs/slac/g/ki/ki09/lances/07Apr26/M13_Dither9_20ArcSecondEast_20ArcSecondsSouth_H4RG_SIPIN_30_Reads_Apr27_2007_03_11_35.fits*',/gui
;   Examine_Pixel_List, objectfile='/nfs/slac/g/ki/ki09/lances/07Apr27/Landolt110956_Dither7_0ArcSecondsEast_20ArcSecondsSouth_H4RG_SIPIN_10_Reads_Apr28_2007_02_53_25.fits*',/gui
;   Examine_Pixel_List, objectfile='/nfs/slac/g/ki/ki09/lances/07Apr27/Landolt110956_Dither7_0ArcSecondsEast_20ArcSecondsSouth_H4RG_SIPIN_10_Reads_Apr28_2007_02_40_47.fits',/gui,xstart=1650,ystart=1850

pro Examine_Pixel_List, Pixels, TvFlag = TvFlag, XStart = XStart, $
	XStop = XStop, YStart=YStart, YStop = YStop , $
	PixelSkip = PixelSkip, DarkFile = DarkFile, $
	ObjectFile = ObjectFile, Date = Date, NReads = NReads, $
	GUI = GUI, PlotIm = PlotIm, SlopeFit = SlopeFit, $
	TvFlagLocal = TvFlagLocal

If N_Elements(PixelSkip) eq 0 then PixelSkip = 0
If N_Elements(GUI) eq 0 then GUI = 0
If N_Elements(PlotIm) eq 0 then PlotIm = 0
If N_Elements(SlopeFit) eq 0 then SlopeFit = 0
If N_Elements(TvFlagLocal) eq 0 then TvFlagLocal = 1

;Coords of bright stars in the default image
If N_Elements(XStart) eq 0 then XStart = 1656
If N_Elements(XStop)  eq 0 then XStop  = XStart+50
If N_Elements(YStart) eq 0 then YStart = 1850
If N_Elements(YStop)  eq 0 then YStop  = YStart+50

;Set up some state variables
Common ExState
Common Images
Common KeyParams, KeyStr
Common SlopeParams, SlopeStr, SlopeFitStr
;Include all the keywords from the KeywordStruct.pro file
@KeywordStruct.pro
@PlotSettings.pro

State = { WindowBase: 0L, $
	  MainDrawID: 0L, $
	  MainWindowID: 0L, $
	  Naxis1      : Naxis1, $
	  Naxis2      : Naxis2, $
	  NumPixels   : 0L, $
	  PixelX      : 0L, $
	  PixelY      : 0L, $
	  PixelXSub   : 0L, $
	  PixelYSub   : 0L, $
	  XStart      : XStart, $
	  YStart      : YStart, $
	  XStop       : XStop, $
	  YStop       : YStop, $
	  Height      : 0L, $
	  Width	      : 0L, $
	  Ref         : Ref ,$
          NReads      : NReads , $
	  SlopeFit    : SlopeFit, $
	  TvFlagLocal : TvFlagLocal }
set_plot,'x'
Userwindow = !d.window
wset, Userwindow

;Initialize the widget
WindowStartUp

;A region was supplied to the script or a list of pixels
If (Size(Pixels))[0] eq 0 then begin 
  State.Width = long(XStop-XStart+1) & State.Height = long(YStop-YStart+1)
  State.NumPixels = long(State.Width)*long(State.Height)

  ;Get the relative coordinates of the pixels in the window
  PixelsSub =  (Indgen(State.NumPixels) mod State.Width) + $
	       (Indgen(State.NumPixels)/ State.Width)*State.Width
  PixelsXSub = PixelsSub mod State.Width &  PixelsYSub=PixelsSub/State.Width
  ;Get the absolute value of the pixels on the full frame
  Pixels =    (Indgen(State.NumPixels) mod State.Width)+XStart+$
	      ((Indgen(State.NumPixels)/State.Width)+YStart)*(Naxis2)
  PixelsX=    Pixels mod Naxis1  &      PixelsY = Pixels/Naxis2

EndIf Else begin 
  State.NumPixels = N_Elements(Pixels)
  XStart = 0  & XStop = Naxis1-1
  YStart = 0  & YStop = Naxis2-1
  PixelsXSub = Pixels mod Naxis1 &  PixelsYSub = Pixels/Naxis2
  PixelsX    = PixelsXSub        &  PixelsY    = PixelsYSub

EndElse

;Open the original object if it's a slope fit
If Not keyword_set(ObjectFile) then $
   ObjectFile = '/nfs/slac/g/ki/ki09/lances/07Apr28/'+$
'M13_Dither3_0ArcSecondsEast_20ArcSecondsNorth_H4RG_SIPIN_20_Reads_Apr29_2007_02_47_59.fits'
DisplayFile=''

;Set the display plot to the slope fit and the object to be posted
If stregex(ObjectFile, 'Slope', /boolean) then begin
   DisplayFile = ObjectFile
   ObjectFile = strmid(ObjectFile, 0,stregex(ObjectFile,'_Slope'))+$
		'.fits'
   ObjectFileParts = strsplit(ObjectFile,'/',/extract, count=nparts)
   ObjectFile  = RootDir+Date+PathDelim+ObjectFileParts(NParts-1)
EndIf 

;Read in the Object File
Position = Fits_Open_DataCube(ObjectFile, ObjectFCB, ObjectHeader, 'Read')
FileParams  = Return_File_Params(ObjectFile)
SlopeStr = Return_Slope_Params(ObjectHeader) 
State.NReads=KeyStr.TotalReads
ObjectIm = Fits_Read_DataCube_Region(ObjectFile, ObjectFCB, ObjectHeader, $
                XStart, XStop, YStart, YStop, 0, KeyStr.TotalReads-1)
ObjectClose = Fits_Close_DataCube(ObjectFCB)
Position = Fits_Open_DataCube(ObjectFile, ObjectFCB, ObjectHeader, 'Read')

;Set the display plot to the slope fit and the object to be posted
If stregex(DisplayFile, 'Slope', /boolean) then begin
   Fits_Read, DisplayFile, DisplayIm
   DisplayIm = DisplayIm(XStart:XStop,YStart:YStop)
EndIf Else begin
   DisplayIm = long(ObjectIm(*,*,State.NReads-1)) - $
	       long(ObjectIm(*,*,0))
EndElse

print, 'Opened Object'
;Get some files for the plotting

If Not keyword_set(DarkFile) then DarkFile   = KeyStr.DarkName 
If Not keyword_set(FlatFile) then FlatFile   = KeyStr.FlatName

;Read in the files that will be used for the examination
Position = Fits_Open_DataCube(DarkFile, DarkFCB, DarkHeader, 'Read')
DarkIm = Fits_Read_DataCube_Region(DarkFile, DarkFCB, DarkHeader, $
		XStart, XStop, YStart, YStop, 0,KeyStr.TotalReads-1)
DarkClose = Fits_Close_DataCube(DarkFCB)
print, 'Opened Dark'

;Read in the Flat File
Position = Fits_Open_DataCube(FlatFile, FlatFCB, FlatHeader, 'Read')
FlatIm = Fits_Read_DataCube_Region(FlatFile, FlatFCB, FlatHeader, $
                XStart, XStop, YStart, YStop, 0, 1)
FlatClose = Fits_Close_DataCube(FlatFCB)
print, 'Opened Flat'

wset, state.mainwindowid
DisplayIm = BytScl(DisplayIm, $
        min = Mean(DisplayIm) - 2*StdDev(DisplayIM), $
        max = Mean(DisplayIm) + 2*StdDev(DisplayIm))
tvscl, congrid(DisplayIm, 512,512)
;tvscl, congrid(ObjectIm(*,*,6)/double(FlatIm(*,*,1)),512,512)
window,1,xsize=500, ysize = 500
;Allow for plotting of the image being viewed to a .png
If PlotIm eq 1 then begin
   Set_plot,'z' 
   device,set_resolution=[512,512]
   DisplayIm = Fltarr(State.width, State.Height,3)
   DisplayIm = [[[ObjectIm(*,*,State.NReads-1)]], $
	        [[ObjectIm(*,*,State.NReads-1)]], $
		[[ObjectIm(*,*,State.NReads-1)]]]
   tvimage, Congrid(bytscl(DisplayIm),512,512,3) ,/tv, true = 3
   PlotName = 'M13Region.png'
   PngImg = tvrd()
   tvlct,reds,greens,blues,/get
   write_png,KeyStr.PlotDir+PlotName, $
                 PngImg,reds,greens,blues
   KeyStr.TvFlag = 1
   Set_plot,'x'
   resetplt
EndIf

;Get the bad pixel maps to check what type of pixel we're plotting
Masks = Return_Pixel_Masks()

;If a list was supplied then plot the pixels in that list
If GUI eq 0 then begin
   For PixInd = 0L, State.NumPixels-1 do begin
    ;Allow for skipping of some pixels
    PixTInd = PixInd + PixelSkip
    PixelX    = PixelsX(PixTInd)      & PixelY    = PixelsY(PixTInd)
    PixelXSub = PixelsXSub(PixTInd)   & PixelYSub = PixelsYSub(PixTInd)
    PlotPixels, PixelXSub, PixelYSub, PixelX, PixelY
    stop
   EndFor
EndIf

If State.TvFlagLocal eq 2 then begin
   set_plot,'ps'
   device, filename = KeyStr.PlotDir+KeyStr.ObjectName+'Slopes.ps'
   loadct, 39
   device, /color                   ;allow color on the postscript
   device,ysize=8.5,/inches        ;Height of plot in y
   device,xsize=11.0,/inches        ;Width of plot in x
   device,yoffset=0.50,/inches      ;Y position of lower left corner
   device,xoffset=0.00,/inches
   !p.multi=[0,1,1]
   !p.charsize=1
   !p.thick=3
   !p.charthick=3
   
Endif
End