;function reffilter.pro
;
;PURPOSE: 
;	To perform a fourier transform on the reference pixels on the H4RG
;	and detect what frequencies are present in the noise.
;
;INPUTS: 
;	InArr - The input array of pixel values (KeyStr.Naxis2xKeyStr.Naxis2 for H4RG)
;
;KEYWORDS:
;	TVFLAG 
;		 0 - Don't Plot anything
;		 1 - Plot to the X Window
;		 2 - Plot a postscript file that was opened by calling script
;		 3 - Plot a png file tha can be converted for a paper 
;
;	AvgEntireFrame
;		 0 - Don't do a fourier transform on entire frame
;		 1 - Do a fourier transform on entire frame
;
;*********************************************************************

;PLOTTING FUNCTIONS ***************************************************
;************************************************PLOT FOR PAPER
function paper_plot, Xvals, Yvals, ReadNum, TitleStr, $
         Units, FileName, FilStr

  Common KeyParams, KeyStr
  Common ColorStr, Color, ColorNum
  erase
  !p.noerase = 0
  plot, XVals(2048:2248), YVals(KeyStr.Naxis2/2:2248), $
        yrange = [0,max(abs(YVals(KeyStr.Naxis2/2+1:KeyStr.Naxis2/2+200)))], $
        xtitle = 'Frequency ' +'('+Units+')',$
        ytitle = 'Power', color = fsc_color('black'), $
	background = fsc_color('white')

  PngImg = tvrd()
  tvlct,reds,greens,blues,/get
  write_png,KeyStr.PlotDir+FileName, $
    PngImg,reds,greens,blues

  return, 1

end

;***********************************************PLOT FOR POSTSCRIPT
function p_in_p_plot, Xvals, Yvals, ReadNum, TitleStr, $
	 Units

  Common KeyParams, KeyStr
  Common ColorStr, Color, ColorNum
  erase
  plot, Xvals,Yvals,$
    yrange=[0,max(abs(Yvals(KeyStr.Naxis2/2+1:KeyStr.Naxis2-1)))],$
    xrange=[0,max(Xvals)], $
    title = TitleStr, $
    xtitle = 'Frequency ' +'('+Units+')',$
    ytitle = 'Power', color = fsc_color('black'), background=fsc_color('white')
    !p.noerase = 1
    plot, XVals(2048:2148), YVals(KeyStr.Naxis2/2:2148), $
          yrange = [0,max(abs(YVals(KeyStr.Naxis2/2+1:KeyStr.Naxis2/2+100)))],$
          position = [0.5, 0.5, 0.9, 0.9], /normal, $
	  color = fsc_color('black'), background = fsc_color('white')

   If KeyStr.TvFlag eq 1 then stop
   return,1

end

;Median Plots ***************************************************
pro plot_median, Columns, Ch0Med

Common KeyParams, KeyStr
Common ColorStr, Color, ColorNum

Med = Median(Columns,dimension = 1)
MedNoOff = Med - Mean(Columns) 
NoZeros = where(Med ne 0)

plot, Med, yrange=[min(Med(NoZeros)),max(Med(NoZeros))], $
         color=fsc_color('black'), background = fsc_color('white')
oplot, Ch0Med,color=fsc_color('black')
If KeyStr.TvFlag eq 1 then stop
plot, MedNoOff, yrange=[min(MedNoOff(NoZeros)),max(MedNoOff(NoZeros))], $
         color=fsc_color('black'), background = fsc_color('white')
oplot, Ch0Med, color=fsc_color('black')
If KeyStr.TvFlag eq 1 then stop

end

;MAIN FUNCTION **************************************************************
function reffilter, InArr, $
	 AvgEntireFrame = AvgEntireFrame, ReadNum = ReadNum, FilStr = FilStr,$
	 ExcludeCh0 = ExcludeCh0, ExamineEachChannel = ExamineEachChannel

Common KeyParams, KeyStr
Common ColorStr, Color, ColorNum

If Not keyword_set(AvgEntireFrame) then AvgEntireFrame = 0
If Not keyword_set(ReadNum) then ReadNum = 0
If Not keyword_set(FilStr) then FilStr = ''
If Not keyword_set(ExcludeCh0) then ExcludeCh0 = 1
If N_Elements(ExamineEachChannel) eq 0 then ExamineEachChannel = 0

;**********************************************************
;Get the temporal frequencies for the Ref Columns at right and left
ColumnTime = 5.4528 			; Seconds taken to read an entire column
CFreqMax = KeyStr.Naxis2/2/ColumnTime		; Column Frequency in Hertz
CFreqSpacing = CFreqMax/(KeyStr.Naxis2/2)    	; Separation of Frequency in Hertz
CFreqs = CFreqSpacing*(-long(KeyStr.Naxis2)/2 + lindgen(KeyStr.Naxis2))

;Get the spatial frequencies for the Ref Rows at top and bottom
RowDim = KeyStr.Naxis2				; The Dimension of the array
RFreqMax = 1./2 			; The maximum frequency in (1/pix)
RFreqSpacing = 1./(RowDim) 		; The frequency spacing in (1/pix)
RFreqs = RFreqSpacing*(-long(KeyStr.Naxis2)/2 + lindgen(KeyStr.Naxis2))

;************************************************************************
;Entire Frame************************************************************
If AvgEntireFrame then begin

   FFTAlongCol = FFT(InArr, dimension = 2)
   FFTAlongColAvg = Avg(abs(FFTAlongCol),0)
   FFTAlongColAvg = shift(temporary(FFTAlongColAvg),KeyStr.Naxis2/2)
   If KeyStr.TvFlag ne 0 then begin
      TitleStr = 'Average Power along all '+strtrim(KeyStr.Naxis2,2) + $
                 ' columns for file ' + KeyStr.RawKey + ' Read # ' + $
              strtrim(ReadNum, 2) + ' ' + FilStr
      Pl = p_in_p_plot(CFreqs, FFTAlongColAvg, ReadNum, TitleStr, 'Hz')         
   EndIf
EndIf

;Individual Channels (0:127) **************************************************
If ExamineEachChannel then begin
  ;Channel 0 
  Ch0 = InArr(0:KeyStr.Ch0W-1, *)
  Ch0Med = Median(Ch0,dimension = 1) - Mean(Ch0)
  plot_median, Ch0, Ch0Med

  ;Left RefPixels
  LeftRef = InArr(0:KeyStr.Ref-1,*)
  Plot_median, LeftRef, Ch0Med

  ;Right RefPixels
  RightRef = InArr(KeyStr.Naxis1-KeyStr.Ref-1:KeyStr.Naxis1-1,*)
  plot_median, RightRef,Ch0Med

  Masks = Return_Pixel_Masks()
  InArrTemp = InArr
  InArrTemp(Masks.BadPixels)=0

  ;Plot the other channels
  For ChN = 1, 31 do begin 
    Ch = InArr(KeyStr.Ch0W*ChN:KeyStr.Ch0W*(ChN+1)-1,*)
    plot_median, Ch, Ch0Med
    ChNoBad = InArrTemp(KeyStr.Ch0W*ChN:KeyStr.Ch0W*(ChN+1)-1,*)   
    plot_median, ChNoBad,Ch0Med
  EndFor

EndIf

;Ch0 Analysis **********************************************************
FFTChannel0 = shift(FFT(InArr(0:KeyStr.Ch0W-1,*),Dimension = 2), $
        0 , KeyStr.Naxis2/2)
FFTChannel0Avg = Avg(abs(FFTChannel0),0)
If KeyStr.TvFlag ne 3 then begin
  TitleStr = 'Average Power along Columns of Channel 0 for '+ $
              KeyStr.RawKey + ' Read # ' + strtrim(ReadNum, 2) + ' ' + FilStr
  Pl = p_in_p_plot(CFreqs, FFTChannel0Avg, ReadNum, TitleStr, 'Hz')
  !p.noerase = 0
  erase
EndIf Else Begin
  TitleStr = 'Average Power along Columns for Channel 0'
  FileName = KeyStr.RawKey+'_'+FilStr+'Ch0PowerSpectrum.png'
  Pl = paper_plot(CFreqs, FFTChannel0Avg, ReadNum, $
      TitleStr, 'Hz', FileName, FilStr)
  !p.noerase =0
  erase
EndElse

;Ch0 With Median Subtracted *********************************************
Ch0 = InArr(0:KeyStr.Ch0W-1, *)
Ch0Med = Median(Ch0,dimension = 1) - Mean(Ch0)
Ch0MinMed = UIntArr(KeyStr.ch0W,KeyStr.Naxis2)
for row = 0, KeyStr.Naxis2-1 do begin
  Ch0MinMed(*,Row) = Ch0(*,Row) - Ch0Med(Row)
EndFor

FFTCh0MinMed = shift(FFT(Ch0MinMed,Dimension = 2), $
        0 , KeyStr.Naxis2/2)
FFTCh0MinMedAvg = Avg(abs(FFTCh0MinMed),0)
If KeyStr.TvFlag ne 3 then begin
  TitleStr = 'Average Power along Columns of Channel 0 Minus Median for '+ $
              KeyStr.RawKey + ' Read # ' + strtrim(ReadNum, 2) + ' ' + FilStr
  Pl = p_in_p_plot(CFreqs, FFTCh0MinMedAvg, ReadNum, TitleStr, 'Hz')
  !p.noerase = 0
  erase
  stop
EndIf Else Begin
  TitleStr = 'Average Power along Columns for Channel 0'
  FileName = KeyStr.RawKey+'_'+FilStr+'Ch0PowerSpectrum.png'
  Pl = paper_plot(CFreqs, FFTCh0MinMedAvg, ReadNum, $
      TitleStr, 'Hz', FileName, FilStr)
  !p.noerase =0
  erase
EndElse


;Left Ref Pixels *********************************************************
LeftRefPix = FltArr(4, KeyStr.Naxis2)
LeftRefPix(*,*) = InArr(0:3, *)

FFTLeftRefPix = fft(LeftRefPix, dimension = 2)
FFTLeftRefPix = shift(temporary(FFTLeftRefPix), 0, KeyStr.Naxis2/2)
FFTLeftRefPixAvg = Avg(abs(FFTLeftRefPix),0)

;Don't consider 0 value
FFTLeftRefPixAvgNoDC = FFTLeftRefPixAvg
FFTLeftRefPixAvgNoDC(KeyStr.Naxis2/2) = 0

;Find the frequencies with too much power
HiFreqInds = where(FFTLeftRefPixAvgNODC gt 1)
If HiFreqInds(0) ne -1 then begin
  Scale = 1./(FFTLeftRefPixAvgNoDC(HiFreqInds))
  Filter = Fltarr(1,KeyStr.Naxis2)
  Filter(*)=1
  Filter(HiFreqInds) = Scale*Filter(HiFreqInds)
EndIf Else Begin
  Scale = 1
  HiFreqInds = Findgen(KeyStr.Naxis2)
  Filter = Fltarr(1,KeyStr.Naxis2)
  Filter(*)=1
EndElse

If KeyStr.TvFlag ne 0 then begin
   TitleStr = 'Average Power along Left Ref Pixels for '+ $
              KeyStr.RawKey + ' Read # ' + strtrim(ReadNum, 2) + ' ' + FilStr
   Pl = p_in_p_plot(CFreqs, FFTLeftRefPixAvg, ReadNum, TitleStr, 'Hz')
   !p.noerase = 0
   erase
end

;Scale the frequencies with alot of power to kill the fixed pattern noise
for Col=0, 3 do begin
   FFTLeftRefPix(Col, HiFreqInds) = Scale*FFTLeftRefPix(Col,HiFreqInds)
endfor

FFTLeftRefPixFil = shift(temporary(FFTLeftRefPix), 0 , -KeyStr.Naxis2/2)
LeftRefPixFil = fft(FFTLeftRefPixFil, dimension =2, /inverse)
FFTLeftRefPixAvg(HiFreqInds) = Scale*FFTLeftRefPixAvg(HiFreqInds)

;Right Ref Pixels ***************************************************
RightRefPix = FltArr(4,KeyStr.Naxis2)
RightRefPix(*,*) = InArr(KeyStr.Naxis2-4:KeyStr.Naxis2-1,*)

FFTRightRefPix = fft(RightRefPix, dimension = 2)
FFTRightRefPix = shift(temporary(FFTRightRefPix), 0, KeyStr.Naxis2/2)
FFTRightRefPixAvg = Avg(abs(FFTRightRefPix),0)

If KeyStr.TvFlag ne 0 then begin
   TitleStr = 'Average Power along Right Ref Pixels for ' + $
	KeyStr.RawKey + ' Read # ' + strtrim(ReadNum, 2) + ' ' + FilStr
   Pl = p_in_p_plot(CFreqs, FFTRightRefPixAvg, ReadNum, TitleStr, 'Hz')
   !p.noerase = 0
   erase

end

;Further Examination later
return, 1

;Bottom Ref Pixels ***************************************************
If ExcludeCh0 then begin
  InArrNoCh0 = InArr(128:KeyStr.Naxis2-1, *)
  KeyStr.Naxis2NoCh0 = KeyStr.Naxis2 - 128
EndIf Else Begin
  InArrNoCh0 = InArr
  KeyStr.Naxis2NoCh0 = KeyStr.Naxis2
EndElse

BotRefPix = FltArr(KeyStr.Naxis2NoCh0, 4)
BotRefPix(*,*) = InArrNoCh0(*, 0:3)

FFTBotRefPix = fft(BotRefPix, dimension = 1)
FFTBotRefPix = shift(temporary(FFTBotRefPix), KeyStr.Naxis2NoCh0/2, 0)
FFTBotRefPixAvg = avg(abs(FFTBotRefPix),1) 

If KeyStr.TvFlag ne 0 then begin
   TitleStr = 'Average Power along Bottom Ref Pixels for '+ $
              KeyStr.RawKey + ' Read # ' + strtrim(ReadNum, 2) + ' ' + FilStr
   Pl = p_in_p_plot(RFreqs, FFTBotRefPixAvg, ReadNum, TitleStr, '1/Pixels')
   !p.noerase = 0
   erase

EndIf

;Top Ref Pixels ******************************************************
TopRefPix = FltArr(KeyStr.Naxis2NoCh0, 4)
TopRefPix(*,*) = InArrNoCh0(*, KeyStr.Naxis1-4: KeyStr.Naxis1-1)

FFTTopRefPix = fft(TopRefPix, dimension = 1)
FFTTopRefPix = shift(temporary(FFTTopRefPix), KeyStr.Naxis2NoCh0/2, 0)
FFTTopRefPixAvg = avg(abs(FFTTopRefPix),1)

If KeyStr.TvFlag ne 0 then begin
   TitleStr = 'Average Power along Top Ref Pixels for '+ $
              KeyStr.RawKey + ' Read # ' + strtrim(ReadNum, 2) + ' ' + FilStr
   Pl = p_in_p_plot(RFreqs, FFTTopRefPixAvg, ReadNum, TitleStr, '1/Pixels')
   !p.noerase = 0
   erase

EndIf

;**********************************************************************
return, shift(Filter, -KeyStr.Naxis2/2)

end