;plot_neg_persistence
;
;PURPOSE:
;  To plot areas of the HxRG detector that have been afflicted by 
;  negative persistence. This script accepts a .fits file that 
;  contains the negative persistence effect as well as a .fits file
;  that contains the saturation area that caused it 
;
;KEYWORDS:
;  MODE:
;	0 - Plot both persistence ramp with negative persistence (NP), positive 
;	    persistence (PP), and sky (Sk) pixels
;	1 - Plot 9 pixels out from the central persistence core
;
;  XPPPIX, YPPPIX:
;	  - the x and y pixel values that contain positive persistence
;  XNPPIX, YNPPIX;
;	  - the x and y pixel values that contain negative persistence
;
;EXAMPLES:
;m36=============
;plot_neg_persistence_1rg, '/nfs/slac/g/ki/ki09/lances/H1RG/KPNO_11_07/ASIC/LEACH/2007Nov13_2/M36/M36_Yfilter_200_1_0_-20_20.fits', '', xpppix = 532, ypppix = 685, xnppix = 540, ynppix = 683, mode=0
;plot_neg_persistence_1rg, '/nfs/slac/g/ki/ki09/lances/H1RG/KPNO_11_07/ASIC/LEACH/2007Nov13_2/M36/M36_Yfilter_200_1_0_-20_20.fits', '', xpppix = 532, ypppix = 685, xnppix = 523, ynppix = 684, mode=0,tvflag=1
;
;SUPER LONG RAMP ======================
;plot_neg_persistence_1rg, '/nfs/slac/g/ki/ki09/lances/H1RG/KPNO_11_07/ASIC/07Nov19/SAO5417_I_NegPerDark_H1RG_SIPIN_2500_Reads_Nov19_2007_22_45_50.fits','', mode = 0, date='07Nov19'
;plot_neg_persistence_1rg, '/nfs/slac/g/ki/ki09/lances/H1RG/KPNO_11_07/ASIC/07Nov19/SAO5417_I_NegPerDark_H1RG_SIPIN_2500_Reads_Nov19_2007_22_45_50.fits','',mode=0
;plot_neg_persistence_1rg, '/nfs/slac/g/ki/ki09/lances/H1RG/KPNO_11_07/ASIC/07Nov19/SAO5417_I_NegPerDark_H1RG_SIPIN_2500_Reads_Nov19_2007_22_45_50.fits','',mode=0, xnppix = 343, ynppix=266, xpppix=336, ypppix = 261,xskpix=336,yskpix=260
;plot_neg_persistence_1rg, '/nfs/slac/g/ki/ki09/lances/H1RG/KPNO_11_07/ASIC/07Nov19/SAO5417_I_NegPerDark_H1RG_SIPIN_2500_Reads_Nov19_2007_22_45_50.fits','',mode=1, xnppix = 343, ynppix=266, xpppix=336, ypppix = 261,xskpix=336,yskpix=260
;plot_neg_persistence_1rg, '/nfs/slac/g/ki/ki09/lances/H1RG/KPNO_11_07/ASIC/07Nov19/SAO5417_I_NegPerDark_H1RG_SIPIN_2500_Reads_Nov19_2007_22_45_50.fits','', mode = 0, date='07Nov19', xnppix = 4, ynppix = 645, xpppix = 8, ypppix = 644 ,xskpix=4,yskpix=644
;NGC419=================================
;plot_neg_persistence_1rg, '/nfs/slac/g/ki/ki09/lances/H1RG/KPNO_11_07/ASIC/07Nov19/NGC7419_I_Dither2_20ArcSecondsEast_20ArcSecondsNorth_H1RG_SIPIN_100_Reads_Nov19_2007_21_02_37.fits', '', xpppix=399, ypppix = 469, xnppix= 402, ynppix = 466, xskpix = 396, yskpix=469, mode=0
pro plot_neg_persistence_1rg, NegFileName, SatFileName, $
	XNPPix=XNPPix, YNPPix = YNPPix, $
	XPPPix=XPPPix, YPPPix = YPPPix, $
        XSkPix=XSkPix, YSkPix = YSkPix, $
	TvFlag = TvFlag, Mode = Mode, Date = Date, Leach=Leach

;Include all the keywords from the KeywordStruct.pro file
Common KeyParams, KeyStr
@KeywordStruct_1RG.pro
@PlotSettings_1RG.pro

;Declare the pixel values that will be used for the positive and negative persistence  
If N_Elements(XNPPix)   eq 0  then XNPPix   = 512
If N_Elements(YNPPix)   eq 0  then YNPPix   = 450
If N_Elements(XPPPix)   eq 0  then XPPPix   = 771
If N_Elements(YPPPix)   eq 0  then YPPPix   = 852
If N_Elements(XSkPix)   eq 0  then XSkPix   = 690
If N_Elements(YSkPix)   eq 0  then YSkPix   = 400
If N_Elements(Mode)     eq 0  then Mode     = 1
If N_Elements(TvFlag)   eq 0  then TvFlag   = 0
If N_Elements(Leach)	eq 0  then Leach    = 0
If SatFileName		eq '' then PlotSat  = 0 $
   			      else PlotSat = 1
;Make the Common Block
Common KeyParams, KeyStr

;The Default File is the one with the negative persistence
KeyStr.RawFullPath = KeyStr.RootDir+KeyStr.PathDelim+KeyStr.Date+KeyStr.PathDelim+NegFileName 
KeyStr.RawFullPath = NegFileName
KeyStr.RawName     = NegFileName
Fits_Read, KeyStr.RawFullPath, NoData, NH, /header_only
FileParams = Return_File_Params_1RG(KeyStr.RawFullPath)
If KeyStr.WindowMode then begin
   XNPix = XPix - KeyStr.XWindowStart(0)
   YNPix = YPix - KeyStr.YWindowStart(0)
EndIf

;Get the exposure time in the negative persistence image
If KeyStr.ITime eq 0 then begin
    If Leach ne 0 then begin
      FrameNTime = (KeyStr.Naxis1*10.e-6+40.e-6)*KeyStr.Naxis2
      INTime     = KeyStr.Naxis3*FrameNTime
    EndIf else begin
      FrameNTime = 0.33 
      INTime     = KeyStr.Naxis3*FrameNTime
    Endelse
EndIf Else Begin
    FrameNTime  = KeyStr.ITime/KeyStr.Naxis3
    INTime	= KeyStr.ITime
EndElse

;***************************************************************************************
Case Mode of 
  0: begin ;**************************************PLOT THREE RAMPS, SKY, NEG, POS.

    ;Ramp of pixel exhibiting positive persistence
    Fits_read_datacube, KeyStr.RawFullPath, ImPP, NH, $
	XStart = XPPPix, XStop = XPPPix, $
        YStart = YPPPix, YStop = YPPPix
    ;Ramp of pixel exhibiting negative persistence 
    Fits_read_datacube, KeyStr.RawFullPath, ImNP, NH, $
        XStart = XNPPix, XStop = XNPPix, $
        YStart = YNPPix, YStop = YNPPix
    ;Ramp of pixel integrating skyflux 
    Fits_read_datacube, KeyStr.RawFullPath, ImSk, NH, $
        XStart = XSkPix, XStop = XSkPix, $
        YStart = YSkPix, YStop = YSkPix

    If PlotSat then begin
      ;Get the info for the saturated ramp
      Fits_read_datacube, SatFileName, ImSPP, SH, $
        XStart = XPPPix, XStop = XPPPix, $ 
        YStart = YPPPix, YStop = YPPPix 
      Fits_read_datacube, SatFileName, ImSNP, SH, $ 
        XStart = XNPPix, XStop = XNPPix, $ 
        YStart = YNPPix, YStop = YNPPix 
      ;Ramp of pixel integrating skyflux  
      Fits_read_datacube, SatFileName, ImSSk, SH, $
        XStart = XSkPix, XStop = XSkPix, $ 
        YStart = YSkPix, YStop = YSkPix 
      ;Get the exposure and frame time for the saturation ramp 
      ISTime      = Double(SxPar(SH, 'ITIME'))
      SNaxis3     = Long(SxPar(SH, 'NAXIS3'))
      FrameSTime  = ISTime/SNaxis3
      STime       = IndGen(N_Elements(ImS))*FrameSTime
    EndIf
    
    ;Plot the Sky, PP, and NP pixels
    erase
    MinPlot   = Min([min(ImPP), min(ImNP), min(ImSk)])
    MinBSPlot = Min([min(ImPP-long(ImPP(0))), min(ImNP-long(ImNP(0))), $
		     min(ImSk-long(ImSk(0)))]) 
    MaxPlot   = Max([max(ImPP), max(ImNP), max(ImSk)])
    MaxBSPlot = Max([max(ImPP-long(ImPP(0))), max(ImNP-long(ImNP(0))), $
		     max(ImSk-long(ImSk(0)))])
    Plot, ImPP, Yrange = [MinPlot, MaxPlot], $
          XTitle = 'Read Number', $
          YTitle = 'ADU' , background = fsc_color('white'), $
          color=fsc_color('black'), position = [0.175, 0.10, 0.9250, 0.90], $
          /normal, xstyle = 8
    OPlot, ImNP, color=fsc_color('red')
    OPlot, ImSk, color=fsc_color('blue') 
    !p.noerase=1
    Axis, XAxis = 1, color=fsc_color('black'), xrange=[0, INTime], $
	  XStyle = 1, xtitle = 'Seconds'
    ItemsN = [Strtrim(XPPPix,2)+','+Strtrim(YPPPix,2), $
	     Strtrim(XNPPix,2)+','+Strtrim(YNPPix,2), $
	     Strtrim(XSkPix,2)+','+Strtrim(YSkPix,2)]
    ColorsN = [fsc_color('black'), fsc_color('red'), fsc_color('blue')]
    TextColorsN =[fsc_color('black'),fsc_color('black'),fsc_color('black')]
    PSymsN = [0,0,0]
    legend, ItemsN, PSym=PSymsN, Color=ColorsN, TextColors = TextColorsN 

    If KeyStr.TvFlag eq 1 then begin
       Img  = tvrd(true=3)
       Img24=reverse(Img,2)
       TVLCT, R, G, B, /Get
       Write_Tiff, KeyStr.PlotDir+KeyStr.RawKey+'3Curves.tif', $
         red=Img24(*,*,0), green=Img24(*,*,1), blue=Img24(*,*,2), $
         planarconfig=2, orientation=1, compression=0
    EndIf
 
    If PlotSat then begin
      Plot, ImS, Yrange = [Min(ImS)-100, Max(Ims)+100], $
	  xstyle = 4, ystyle = 4, $
          color = fsc_color('black'), background = fsc_color('white'), $
          position = [0.125, 0.10, 0.90, 0.90], $
          /normal, XRange = [0, INTime]
      Axis, YAxis = 1, color=fsc_color('black'), yrange=[Min(ImS)-100, Max(ImS)+100], $
          ystyle =1, YTitle = 'ADU in First Ramp'
    EndIf 

    !P.Noerase = 0
    If KeyStr.TvFlag eq 3 then begin
           PNGFileName = KeyStr.RawKey+'_'+Strtrim(XPix, 2)+'_'+Strtrim(YPix,2)+'.png'
           PNGImg      = tvrd()
           tvlct, reds, greens, blues, /get
           write_png, KeyStr.PlotDir+PNGFileName, PNGImg, Reds, Greens, Blues
           stop
     EndIf
     stop
  end
  1: begin ;**************************PLOT 9 pixels out from core
    
    ;Ramp of pixel exhibiting positive persistence
    NumPix = 10
    Fits_read_datacube, KeyStr.RawFullPath, Im, NH, $
        XStart = XPPPix, XStop = XPPPix+NumPix-1, $
        YStart = YPPPix, YStop = YPPPix
    NumReads = Long(SxPar(NH, 'NAXIS3'))
    ImMB     = LonArr(NumPix,1,NumReads)
    ImB	     = Long(Im(*,0,0))
    For ReadNum = 0, NumReads-1 do begin
        ImMB(*,0,ReadNum) = Im(*,0,ReadNum)-ImB
    EndFor
    MinImMB    = Min(ImMB)
    MaxImMB    = Max(ImMB)
    Items      = StrArr(NumPix)
    PSyms      = LonArr(NumPix)
    TextColors = LonArr(NumPix)
    PColors    = LonArr(NumPix) 
    LThicks    = LonArr(NumPix)
    LStyles    = LonArr(NumPix) 
    Plot, [0,0], /nodata, background = fsc_color('white'), $
	yrange = [MinImMB, MaxImMB], xrange = [0, NumReads], $
        color  = fsc_color('black'), xtitle ='Read Number', $
	ytitle = 'Bias Subtracted ADU', $
        position = [0.1, 0.1, 0.8, 0.95], /normal, xstyle = 8 
    Axis, XAxis = 1, color=fsc_color('black'), xrange=[0, INTime], $
          XStyle = 1, xtitle = 'Seconds'
    For PixNum = 0, NumPix-1 do begin
      oplot, ImMB(PixNum, 0, *), color = ColorNum.ColorNums(PixNum+1)
      Items(PixNum)      = 'Rad ' + Strtrim(PixNum, 2)
      PColors(PixNum)     = ColorNum.ColorNums(PixNum+1)
      TextColors(PixNum) = fsc_color('black')
      PSyms(PixNum)      = 0  
      LThicks(PixNum)	 = 5
      LStyles(PixNum)    = 0
    EndFor
    legend, Items, Color = PColors, TextColor = TextColors, PSym = PSyms, $
       position = [0.825, 0.9], /normal, Thick = LThicks, LineStyle = LStyles
    If KeyStr.TvFlag eq 3 then begin
      Img  = tvrd(true = 3)
      Img24=reverse(Img,2)
      TVLCT, R, G, B, /Get
      Write_Tiff, KeyStr.PlotDir+KeyStr.RawKey+'RadialNegPersistPlot.tif', $
         red=Img24(*,*,0), green=Img24(*,*,1), blue=Img24(*,*,2), $
         planarconfig=2, orientation=1, compression=0

      PngImg = tvrd()
      tvlct,reds,greens,blues,/get
      write_png,KeyStr.PlotDir+KeyStr.RawKey+'RadialNegPersistPlot.png', $
      PngImg,reds,greens,blues
    EndIf 
    stop
  end
  2: begin
    ;Plot just the persistent ramp
    ;Get info for single pixel of window mode or negative persistence image
    Fits_read_datacube, KeyStr.RawName, ImN, NH, $
        XStart = XPix, XStop = XPix, $
        YStart = YPix, YStop = YPix 
    Plot, ImN, Yrange = [Min(ImN), Max(ImN)], $
          XTitle = 'Read Number', $
          YTitle = 'ADU in Second Ramp' , background = fsc_color('white'), $
          color=fsc_color('black'), position = [0.125, 0.10, 0.90, 0.90], $
          /normal, xstyle = 8, ystyle = 8
    !p.noerase=1
    XYOuts, 0.3, 0.75, 'Signal for Pixel '+Strtrim(Xpix,2)+' , '+Strtrim(YPix,2), $
          /normal, color=fsc_color('black')
    Axis, XAxis = 1, color=fsc_color('black'), xrange=[0, INTime], XStyle = 1, $
          xtitle = 'Seconds'
    stop
  end
  3: begin
     ; Do radial analysis of the negative persistence
     FitsFileName = '~/Desktop/PersistenceImages/SAO5417_I_NegPerSatRamp_H1RG_SIPIN_500_Reads_Nov19_2007_22_09_16.fits'
     FitsFileName = '~/Desktop/PersistenceImages/HD34364NegativePersistenceSaturation_H1RG_SIPIN_OneWindow5000_Reads_Nov18_2007_02_38_09.fits' & XPix = 47 & YPix = 44
    
     BoxSize = 50
     Cols = (Rows = BoxSize)
     Fits_Read_DataCube, FitsFileName, Im, FitsHeader, $
       XStart = XPix-BoxSize/2, XStop = XPix+BoxSize/2,$
       YStart = YPix-BoxSize/2, YStop = YPix+BoxSize/2
     ITIME  = Float(SxPar(FitsHeader, 'ITIME'))     
     Naxis1 = Long(SxPar(FitsHeader, 'NAXIS1'))
     Naxis2 = Long(SxPar(FitsHeader, 'NAXIS2'))
     Depth  = Long(SxPar(FitsHeader, 'NAXIS3'))
     FrameTime = ITIME/Depth
     RelXCoords = Indgen(Cols,Rows) mod Cols
     RelXCoords = RelXCoords - Mean(RelXCoords)
     RelYCoords = reverse(Indgen(Cols, Rows),2) / Cols
     RelYCoords = RelYCoords - Mean(RelYCoords)
     RsqCoords  = RelXCoords^2+RelYCoords^2
     RadMax     = sqrt((Cols/2)^2+(Rows/2)^2)
     ;Plot one set of slopes for each radius     
     For Rad = 0, 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  =   RelXCoords(RelPixInRad) + Cols/2 +$
         ;                (RelYCoords(RelPixInRad)+ Rows/2)*Cols
         AbsPixInds   = RelPixInRad
         For PixNum = 0, N_Elements(AbsPixInds)-1 do begin
           XInd = long(AbsPixInds(PixNum)) mod Cols
           YInd = long(AbsPixInds(PixNum)) / Cols
           MinD = min(Im(XInd,YInd,*))
           MaxD = max(Im(XInd,YInd,*))
	   ImMB = Im(XInd,YInd,*)-Long(Im(XInd,YInd,0))
           MinImMB = min(ImMB)
           MaxImMB = max(ImMB)
           XMin = where(ImMB eq MinImMB)
           print, 'Col : ' + Strtrim(XInd, 2)
           print, 'Row : ' + Strtrim(YInd, 2)
	   print, 'Radius : ' + Strtrim(sqrt(RSqCoords(XInd,YInd)),2)
           Plot,  ImMB, $
	     Yrange = [Min(ImMB), Max(ImMB)], $
             XTitle = 'Read Number', $
             YTitle = 'ADU in Second Ramp' , background = fsc_color('white'), $
             color=fsc_color('black'), position = [0.125, 0.10, 0.90, 0.90], $
             /normal, xstyle = 8, ystyle = 8
          OPlot, [XMin, XMin] , [!Y.CRANGE(0), !Y.CRANGE[1]], linestyle = 3, $
            psym=10, color=fsc_color('red')
          XYOuts, 0.3, 0.75, 'Signal for Pixel '+ $
	    Strtrim(XInd,2)+' , '+Strtrim(YInd,2), $
            /normal, color=fsc_color('black')
	  XYOUTS, 0.3, 0.70, 'At a radius of : ' + $
	     Strtrim(sqrt(RSqCoords(XInd,YInd)),2), $
	    /normal, color=fsc_color('black')
          FrameTime = (Naxis1*10.e-6+40.e-6)*Naxis2
          ITime     = Depth*FrameTime
	  Axis, XAxis = 1, color=fsc_color('black'), xrange=[0, ITime], $
	    XStyle = 1, xtitle = 'Seconds'
           stop
       EndFor 
     EndFor
  end
  4: begin
    For YPix = 40, 89 do begin
      For x = 40, 50 do begin
        MinD = min(ImN(x,YPix,*))
        MaxD = max(ImN(x,YPix,*))
        MinT = where(ImN(x,YPix,*) eq MinD)
        MinArr = [MinArr, MinT]
        print, 'Min in Column '+strtrim(x, 2) + ' : ' + $
         strtrim(MinT,2)
        plot, Lindgen(Depth)*FrameTime, $
         ImN(x, YPix, *), yrange = [MinD, MaxD]
        stop
      EndFor
    Endfor

    MinHist = Histogram(MinArr, locations = LocMin)
    TimeMin = FrameTime * LocMin
    plot, TimeMin, MinHist
  end

endcase

stop

end