;Single_Window_Modde
;
;PURPOSE:
;	This script will return the difference in coordinates between a 
;	ShackHartman image centroid and the corresponding point on the 
;	reference grid
;
;	Since the center of the grid seems to be drifting 
;
;KEYWORDS:
;	FIRSTREAD: The first read to be used when examining the files
;
;	TVFLAG: 
;	1- Use the xwindow and plot centroids
;	2- Use the postscript
;
;	PSWINDOWS:
;          If TvFlag eq 2, this will display windows in the final postscript
;          that are contained in the readnum array
pro Single_Window_Mode, FitsFileName, TvFlag=TvFlag, FirstRead=FirstRead,$
    PSWindows = PSWindows, Rem60Hz=Rem60Hz, PlotStackRad = PlotStackRad

If N_Elements(Rem60Hz)      eq 0 then Rem60Hz   = 0
If N_Elements(PSWindows)    eq 0 then PSWindows = 0
If N_Elements(TvFlag)       eq 0 then TvFlag    = 0
If N_Elements(FirstRead)    eq 0 then FirstRead = 0
If N_Elements(PlotStackRad) eq 0 then PlotStackRad = 0
If FitsFileName eq '' then FitsFileName = $
   '/nfs/slac/g/ki/ki09/lances/07Apr29/HD103578_Focus_H4RG_SIPIN_OneWindow300_Reads_Apr29_2007_20_57_27.fits*'
@KeywordStruct.pro
@PlotSettings.pro
;Open the FitsFile to the first valid data byte after the header.
;The common file control block FCB holds the structure to the fits file
Position = fits_open_datacube(FitsFileName, FCB, Header, 'Read')

!p.noerase = 1 
If KeyStr.TvFlag eq 1 then window, xsize = 512, ysize =512

Naxis1 = fcb.axis(0)
Naxis2 = fcb.axis(1)
Naxis3 = fcb.axis(2)
Naxis4 = fcb.axis(3)
bitpix = fcb.bitpix(0)
bzero = fcb.bzero(0)
bytes_per_word = (abs(bitpix)/8)
SS    = 250
YRatF  = SS/Naxis1 & XRatF = SS/Naxis2

;Zero the bad pixels in each window 
Masks        = Return_Pixel_Masks()
SubMask      = BytArr(KeyStr.Naxis1,KeyStr.Naxis2)
SubMasks     = BytArr(Naxis1,Naxis2)

;Get the window coordinates
WindowNum=0
XMin   = Long(SxPar(Header, 'WNDW'+Strtrim(WindowNum,2)+'XMI'))
XMax   = Long(SxPar(Header, 'WNDW'+Strtrim(WindowNum,2)+'XMA'))
YMin   = Long(SxPar(Header, 'WNDW'+Strtrim(WindowNum,2)+'YMI'))
YMax   = Long(SxPar(Header, 'WNDW'+Strtrim(WindowNum,2)+'YMA'))
XHWid  = (Naxis1)/2  & YHWid = (Naxis2)/2
XCen   = (Naxis1)/2
YCen   = (Naxis2)/2
Border = 2
SubMask(XMin-1:XMax-1,YMin:YMax)=1
SubMask  = SubMask*Masks.BadMask
SubMasks(*,*) = ~SubMask(XMin-1:XMax-1, YMin:YMax)
CDS           = Fltarr(Naxis1,Naxis2,2)

;Adjust Minimum and Maximum
XMin  = 10  & XMax = 40 & XCen = (XMin+XMax)/2
YMin  = 10  & YMax = 40 & YCen = (YMin+YMax)/2
XMin0 = XMin& XMax0=XMax
YMin0 = YMin& YMax0=YMax
BXWid = XMax-XMin  & BYWid = YMax-YMin 
CentroidBox = Fltarr(BXWid,BYWid)

;Make arrays for the centroid adjusted stack and the undadjusted stack
CStack         = Dblarr(BXWid+1,BYWid+1)
NumCStack      = Dblarr(BXWid+1,BYWid+1)
UStack         = Dblarr(Naxis1,Naxis2)
NumUStack      = Dblarr(Naxis1,Naxis2)
;Keep track of things w/respect to first read
ESArr          = DblArr((Naxis3-FirstRead)/2)
XSCenArr       = DblArr((Naxis3-FirstRead)/2)
YSCenArr       = DblArr((Naxis3-FirstRead)/2)
X2SArr         = DblArr((Naxis3-FirstRead)/2)
Y2SArr         = DblArr((Naxis3-FirstRead)/2)

EIArr          = DblArr((Naxis3-FirstRead)/2)
XICenArr       = DblArr((Naxis3-FirstRead)/2)
YICenArr       = DblArr((Naxis3-FirstRead)/2)
X2IArr         = DblArr((Naxis3-FirstRead)/2)
Y2IArr         = DblArr((Naxis3-FirstRead)/2)

EUArr          = DblArr((Naxis3-FirstRead)/2)
XUCenArr       = DblArr((Naxis3-FirstRead)/2)
YUCenArr       = DblArr((Naxis3-FirstRead)/2)
X2UArr         = DblArr((Naxis3-FirstRead)/2)
Y2UArr         = DblArr((Naxis3-FirstRead)/2)

TimeArr        = DblArr((Naxis3-FirstRead)/2)
NumOff         = 0.

;Window Time
Tp = 10.e-6              ;Pixel time in seconds
Tr = (Naxis1*Tp+18.e-6)  ;Row time in seconds
Tw = Naxis2*Tr+18.e-6   ;Window time in seconds
Tc = 3.*Tw                ;Cycle time in seconds

If PSWindows then begin
   DispReadNums = [821,853,887,929,0]  ; 821, 807, 929, 949, 959
   NumDisp      = N_Elements(DispReadNums)
   DispWindows  = Fltarr(Naxis1,Naxis2,NumDisp)
   DispTime     = Fltarr(NumDisp)
   DispNum      = 0
EndIf

For ReadNum = FirstRead, Naxis3-1 do begin

  Fits_Read_Datacube,FitsFileName, Data, Header, $
           ZStart = ReadNum, $
           ZStop  = ReadNum
  CDS(*,*,ReadNum mod 2) = Data * SubMasks(*,*)

  If (ReadNum mod 2) eq 1 then begin
     ThisWindow = CDS(*,*,1)-CDS(*,*,0)
     If Rem60Hz then begin
          Noise = Median(ThisWindow(60:*,*),dimension = 1)
          For i = 0, Naxis2-1 do ThisWindow(*,i) = ThisWindow(*,i)-Noise(i)
     EndIf

     UStack(*,*) = UStack(*,*)+double(ThisWindow)
     NumUStack(*,*) = NumUStack(*,*)+SubMasks(*,*)
     TimeArr((ReadNum-FirstRead)/2)   = (ReadNum-FirstRead)/2*Tc
     If PSWindows then begin 
        If ReadNum eq DispReadNums(DispNum) then begin
           DispWindows(*,*,DispNum) = ThisWindow
           DispTime(DispNum) = (ReadNum-FirstRead)/2*Tc
           DispNum+=1
        EndIf
     EndIf
     MinWindow  = Min(ThisWindow(where(ThisWindow ne 0)))
     MeanWin    = Mean(ThisWindow(where(ThisWindow ne 0)))
     StdDevWin  = StdDev(ThisWindow(where(ThisWindow ne 0)))
     StarCen = [4,4]
     XMin  = Double(10)  & XMax = Double(40) & XCen = (XMin+XMax)/2
     YMin  = Double(10)  & YMax = Double(40) & YCen = (YMin+YMax)/2
     Star    = ThisWindow(XMin:XMax,YMin:YMax)
     while abs(StarCen(0)) gt 1 or abs(StarCen(1)) gt 1 do begin
       StarMin = Min(Star(where(Star ne 0)))
       StarCen = Return_Centroid(Star-StarMin, /RemoveBad)
       XMin    = (XMin+StarCen(0)) > 2 & XMax = (XMax+StarCen(0)) < (Naxis1-1) 
       YMin    = (YMin+StarCen(1)) > 2 & YMax = (YMax+StarCen(1)) < (Naxis2-1)
       Star    = ThisWindow(XMin:XMax, YMin:YMax)
       If KeyStr.TvFlag eq 1 then begin
         tv, bytscl(congrid(Star,SS,SS), $
             min = MeanWin - 2*StdDevWin, $
             max = MeanWin + 2*StdDevWin)
       EndIf
       print, StarCen
     endwhile

     If KeyStr.TvFlag eq 1 then begin
       XCircF = ((XMin+BXWid/2+3))*XRatF
       YCircF = ((YMin+BYWid/2+3))*YRatF
       XRatL  = SS/(XMax-XMin) & YRatL = SS/(YMax-YMin)
       XCircL = (BXWid/2)*XRatL
       YCircL = (BYWid/2)*YRatL
       tv, bytscl(congrid(Star,SS,SS), $
             min = MeanWin - 2*StdDevWin, $
             max = MeanWin + 2*StdDevWin)
       tv, bytscl(congrid(ThisWindow,SS,SS), $
             min = MeanWin - 2*StdDevWin, $
             max = MeanWin + 2*StdDevWin),$
             SS, 0
       tvcircle, 3, YCircL, YCircL, fsc_color('green'), /fill
       tvcircle, 3, XCircF+SS, YCircF, fsc_color('green'),/fill
     EndIf
     ;Add to the CStack
     if (XMax-Xmin) eq BXWid and (YMax-YMin) eq BYWid then begin
       CStack(*,*)    = CStack(*,*)+Star
       NumCStack(*,*) = NumCStack(*,*)+SubMasks(XMin:XMax,YMin:YMax)
       TempDisSt = CStack/NumCStack
       BadLocs = where(finite(TempDisSt) eq 0, NumBad)
       If NumBad gt 0 then TempDisSt(BadLocs)=0
       MinSt = Min(TempDisSt(where(TempDisSt gt 0)))
       CStackCen = Return_Centroid(TempDisSt-MinSt, /RemoveBad)

       ;Record Values for the uncentered stack
       MinStack = Min(UStack(where(UStack gt 0)))
       UStackCen = Return_Centroid(UStack-MinStack, /RemoveBad)
       EUArr((ReadNum-FirstRead)/2)    = UStackCen(7)
       XUCenArr((ReadNum-FirstRead)/2) = XMin
       YUCenArr((ReadNum-FirstRead)/2) = YMin
       X2UArr((ReadNum-FirstRead)/2)   = UStackCen(2)
       Y2UArr((ReadNum-FirstRead)/2)   = UStackCen(3)
     
       ;Record Values for the stack
       ESArr((ReadNum-FirstRead)/2)    = CStackCen(7)
       XSCenArr((ReadNum-FirstRead)/2) = XMin
       YSCenArr((ReadNum-FirstRead)/2) = YMin
       X2SArr((ReadNum-FirstRead)/2)   = CStackCen(2)
       Y2SArr((ReadNum-FirstRead)/2)   = CStackCen(3)

       ;Record Values for the individual CDS frame
       EIArr((ReadNum-FirstRead)/2)    = StarCen(7)
       XICenArr((ReadNum-FirstRead)/2) = XMin
       YICenArr((ReadNum-FirstRead)/2) = YMin
       X2IArr((ReadNum-FirstRead)/2)   = StarCen(2)
       Y2IArr((ReadNum-FirstRead)/2)   = StarCen(3)
       print, CStackCen
       If KeyStr.TvFlag eq 1 then begin
         tvscl, congrid(TempDisSt,SS,SS), $
      	      0, SS
         wait,.1
       EndIf
     EndIf Else Begin
       NumOff +=1.
     EndElse
  EndIf
endfor

GoodVal = where(ESArr gt 0, NumGood, complement=BadVal) 
If NumGood gt 0 then begin
   ESArr    = ESArr(GoodVal)
   XSCenArr = XSCenArr(GoodVal)
   YSCenArr = YSCenArr(GoodVal)
   X2SArr   = X2SArr(GoodVal)
   Y2SArr   = Y2SArr(GoodVal)

   EIArr    = EIArr(GoodVal)
   XICenArr = XICenArr(GoodVal)
   YICenArr = YICenArr(GoodVal)
   X2IArr   = X2IArr(GoodVal) 
   Y2IArr   = Y2IArr(GoodVal)
   TimeArr  = TimeArr(GoodVal)
EndIf

If PlotStackRad eq 1 then begin
  
  ;SUM OF THE STACK ***********
  UStackCen = Find_Centroid(UStack, XMin0, XMax0, YMin0, YMax0,/global)
  UStackIm  = UStack(Xmin0:XMax0, YMin0:YMax0)
  ;Get the radii and intensities for the centered and uncentered
  Plot_Radial_Star, UStackIm, URadii, UInten, /RjctOut
  Plot_Radial_Star, CStack, CRadii, CInten, /RjctOut
  If KeyStr.TvFlag eq 2 then begin
     device, filename = KeyStr.PlotDir+'HD103578_SumStack.eps', /encapsul
  EndIf
  !p.noerase=0 
  MaxInt = Max([Max(UInten), Max(CInten)])
  plot,sqrt(uradii(1:*)),uinten(1:*),psym=2, background = fsc_color('white'), $
     color  = fsc_color('black'), $
     yrange = [0, MaxInt], xtitle ='Radius (pixels)', ytitle = 'Counts',$
     position = [0.1,.1,0.95,0.95], /normal
  oplot,sqrt(cradii(1:*)),cinten(1:*),psym=2,color=fsc_color('red')
  If KeyStr.TvFlag eq 2 then begin
     device, /close
     device, filename = KeyStr.PlotDir+'HD103578_MeanStack.eps', /encapsul
  EndIf 

  ;MEAN OF THE STACK************
  UStackMean   = UStack/NumUStack
  CStackMean   = CStack/NumCStack
  UStackMeanIm = UStackMean(Xmin0:XMax0, YMin0:YMax0)
  UStackMeanIm(where(finite(UStackMeanIm) eq 0))=0

  ;Get the FWHMs 
  MinUStackMean = Min(UStackMeanIm(where(UStackMeanIm ne 0)))
  UCen  = Return_Centroid(UStackMeanIm-MinUStackMean)
  UFWHM  = 2.3548*sqrt((UCen(2)+UCen(3))/2)
  UFWHM2 = URadii((where(UStackMeanIm ge .5*Max(UStackMeanIm)))[0])
  
  MinCStackMean = Min(CStackMean(where(CStackMean ne 0)))
  CCen  = Return_Centroid(CStackMean-MinCStackMean)
  CFWHM = 2.3548*sqrt((CCen(2)+CCen(3))/2)
  CFWHM2 = CRadii((where(CStackMean ge .5*Max(CStackMean)))[0])

  ;Get the radii and intensities for the centered and uncentered
  Plot_Radial_Star, UStackMeanIm, URadii, UInten, /RjctOut
  Plot_Radial_Star, CStackMean, CRadii, CInten, /RjctOut

  !p.noerase=0
  MaxInt = Max([Max(UInten), Max(CInten)])
  plot,sqrt(uradii(1:*)),uinten(1:*),psym=2, background = fsc_color('white'), $
     color  = fsc_color('black'), $
     xrange = [0, 20], xstyle = 1, $
     yrange = [0, MaxInt], xtitle ='Radius (pixels)', ytitle = 'Counts', $
     position = [0.15,.1,0.95,0.95], /normal
  oplot,sqrt(cradii(1:*)),cinten(1:*),psym=2,color=fsc_color('red')
  items = ['Uncentered Stack', 'Centered Stack']
  pcolors= [fsc_color('black'), fsc_color('red')]
  ptextcolors = [fsc_color('black'), fsc_color('black')]
  psyms = [2,2]
  pthick = [4,4]
  pcharsize = [2,2]
  legend, items, colors = pcolors, textcolors=ptextcolors, psym=psyms, $
          position = [0.525,0.85],/normal, thick = pthick, charsize=3
  If KeyStr.TvFlag eq 2 then begin
     device, /close
     device, filename = KeyStr.PlotDir+'HD103578_CentroidMot.eps', /encapsul
  EndIf
  !p.multi=[0,1,1]
  ArcAx = 0.125*(findgen(23)-11)
  Plot, timearr, xicenarr-mean(xicenarr), ytitle = 'Pixels',$
   xtitle = 'Time (s)', background=fsc_color('white'), color=fsc_color('black'), $
   xrange=[0, 32], xstyle=1, $
   yrange=[-11,11], ystyle=9, $
   position = [0.1,0.1,0.9,0.95]
  Axis, yaxis=1, color=fsc_color('black'), yrange = [min(ArcAx), max(ArcAx)], $
   ytitle = 'Arcseconds'
  OPlot, timearr, yicenarr-mean(yicenarr), color = fsc_color('red')
  items = ['X', 'Y']
  pcolors= [fsc_color('black'), fsc_color('red')]
  ptextcolors = [fsc_color('black'), fsc_color('black')]
  psyms = [0,0]
  pthick = [4,4]
  pcharsize = [2,2]
  legend, items, colors = pcolors, textcolors=ptextcolors, psym=psyms, $
          position = [0.25,0.275],/normal, thick = pthick, $
          charsize = 3
  If KeyStr.TvFlag eq 2 then begin
     device, /close
  EndIf

  stop


EndIf

If KeyStr.TvFlag eq 1 then begin 
  window,1,xsize=1600,ysize=512
  !p.noerase=0
  plot,TimeArr, XICenarr,psym=2 ,color=fsc_color('black'),$
     background=fsc_color('white')
  oplot, TimeArr, XICenarr,psym=0, color=fsc_color('black')
  oplot, TimeArr, YICenarr,psym=2,color=fsc_color('red')
  oplot, TimeArr, YICenarr,psym=0,color=fsc_color('red')


  window,1,xsize=1600,ysize=512
  !p.noerase=0 
  plot,TimeArr, EIArr,psym=2 ,color=fsc_color('black'),$
     background=fsc_color('white')
  oplot, TimeArr, EIArr,psym=0, color=fsc_color('black')

EndIf

If KeyStr.TvFlag eq 2 then begin
  device, filename=KeyStr.PlotDir+'113095Focusing.eps',/encapsul
 
   
  XIFWHM = 2.3548*sqrt(X2IArr)*0.125
  YIFWHM = 2.3548*sqrt(Y2IArr)*0.125
  MinIFWHM = Min([XIFWHM,YIFWHM])
  MaxIFWHM = Max([XIFWHM,YIFWHM])
  !p.thick = 4.0
  plot,TimeArr, XIFWHM,psym=2 ,color=fsc_color('black'),$
     background=fsc_color('white'),$
     yrange=[MinIFWHM, MaxIFWHM], $
     xtitle='Seconds', ytitle = 'Arcseconds'
  !p.thick = 4.0
  oplot, TimeArr, XIFWHM,psym=0, color=fsc_color('black')
  !p.thick = 4.0
  oplot, TimeArr, YIFWHM,psym=2,color=fsc_color('red')
  !p.thick = 4.0
  oplot, TimeArr, YIFWHM,psym=0,color=fsc_color('red')
  !p.noerase=1
  ;Put the Timing info with some pictures of the focusing star
  For WinNum = 0, NumDisp-2 do begin
    ThisWindow =  DispWindows(*,*,WinNum)
    BWIm = Fltarr(Naxis1,Naxis2,3)
    BWIm(*,*,0)=ThisWindow
    BWIm(*,*,1)=ThisWindow
    BWIm(*,*,2)=ThisWindow
    xyouts, 0.1*(2*WinNum+1)+0.075, 0.875, $
            't = '+Strtrim(string(DispTime(WinNum),format='(F10.4)'),2)+' s',$
            color=fsc_color('black'),/normal
    tvimage, bytscl(BWIm), $
	position = [0.1*(2*WinNum+1)+.075, 0.625, $
                    0.1*(2*WinNum+1)+.245 ,0.85],/normal
  EndFor
  ;Make a legend
  items   = ['FWHM in X', 'FWHM in Y']
  syms    = [0, 0]
  pcolors = [fsc_color('black'), fsc_color('red')]
  texts   = [fsc_color('black'), fsc_color('black')]
  legend, items, psym =syms, colors = pcolors, textcolors=texts,$
          position = [0.175, 0.2], /normal
  device,/close
Endif
  stop


end