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

  Common KeyParams, KeyStr
  Common ColorStr, Color, ColorNum, ColorTri, ColorNumTri

  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, ColorTri, ColorNumTri

  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, ColorTri, ColorNumTri

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

;pro Plot_AllNights - When the input to Examine_Ref_Pixels is a list of 
;	files, the median of Ch0 will be plotted for all of the nights on
;	one plot
pro Plot_AllNights, FitsFiles, ReadNum = ReadNum

Common KeyParams, KeyStr
Common ColorStr, Color, ColorNum, ColorTri, ColorNumTri

If N_Elements(ReadNum) eq 0 then ReadNum =0
Readcol, FitsFiles, RawFiles, format='A'

If KeyStr.TvFlag eq 1 then window,xsize=1000,ysize=1000
plot, [0], [0], /nodata, background = fsc_color('white'),$
	color=fsc_color('black'),$
	yrange=[-200,200], xrange=[0, KeyStr.Naxis2]
;!p.noerase=1
!p.multi=[0,2,2]
For FileNum = 0, N_Elements(RawFiles)-1 do begin
    
    Fits_Read_DataCube, RawFiles(FileNum), Ch0, RawHeader, $
	 XStart = 0, XStop = KeyStr.Ch0W-1, YStart=0, $
	 YStop = KeyStr.Naxis2-1, ZStart = ReadNum, ZStop = ReadNum

    Ch0Med = Median(Ch0, dimension = 1)
    Ch0MedNoOffset =  Ch0Med - Mean(Ch0)

    ;KPNO Data are different for different gain.  Individual nights:
    If FileNum eq 0 or FileNum eq 2 or FileNum eq 4 then begin
      plot, Ch0MedNoOffset, color = fsc_color('black'), $
	background = fsc_color('white')
        ColorInd = 1
    EndIf else begin
      ColorInd +=1
      oplot, Ch0MedNoOffset, color = ColorNum.ColorNums(ColorInd)
    EndElse   
EndFor

xyouts,0.55,0.45,'Upper Left Plot', $
	/normal, color=fsc_color('black')
legend,['07Apr24','07Apr25'], psym=[0,0],$
	colors=[fsc_color('black'),ColorNum.ColorNums(2)], $
	position = [0.575,0.42], thick = [4,4],  /normal, $
	textcolors=[fsc_color('black'),fsc_color('black')] 
xyouts,0.55,0.35,'Upper Right Plot', /normal, color=fsc_color('black')
legend,['07Apr26','07Apr27'], psym=[0,0],$
        colors=[fsc_color('black'),ColorNum.ColorNums(2)], $
        position = [0.575,0.32], thick = [4,4], /normal, $
        textcolors=[fsc_color('black'),fsc_color('black')] 
xyouts,0.55,0.25,'Lower Left Plot', /normal, color=fsc_color('black')
legend,['07Apr28','07Apr29','07Apr30'], psym=[0,0,0],$
        colors=[fsc_color('black'),ColorNum.ColorNums(2), $
		ColorNum.ColorNums(3)], $
        position = [0.575,0.22], thick=[4,4,4],/normal, $
        textcolors=[fsc_color('black'),fsc_color('black'),fsc_color('black')] 

If KeyStr.TvFlag eq 3 then begin
  PngImg = tvrd()
  tvlct,reds,greens,blues,/get
  write_png,KeyStr.PlotDir+'Ch0NoiseAllNights.png', $
    PngImg,reds,greens,blues
EndIf

!p.noerase=0
set_plot,'x'

stop

end

;pro Examine_Ref_Pixels
;
;PURPOSE:
;	To take an image and produce quantitative measure-
;	ments of the reference pixels and channel 0 using the fast fourier 
;	transform.
;
;INPUTS:
;	FITSFILENAME:
;	  The File Name for which the analysis should be done
; 
;KEYWORDS:
;       ROOTDIR: 
;         The Directory holding the directories of Raw files
;       DATE:
;         The Date during which they were taken.  Also the subdirectory that 
;         the files are in.
;       REDUCEDDIR:
;         The directory where the processed images and plots will be placed in.
;       BLOCKLENGTH:
;         An offset for files that were written improperly.  The readu 
;         function reads the values with this offset in place.
;       OBJECTNAME: 
;         The string that indicates what files are being sought.  The filenames
;         should begin with this string ('Dark', 'Flat', etc.).
;       NREADS:
;         The number of reads (NAXIS3) that are present in the file.  This
;         is also a string in the filename that will be used as a regular
;         expression in searching for matched files.
;       FIRSTREAD:
;         The first read to be used in the difference 
;             Diff = LastRead - FirstRead
;       LASTREAD: 
;         The last read to be used in the difference 
;       THISFILTER:
;         The filter through which the files were taken. ('I', 'G', or 'Y')
;       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
;	REMOVALMODE
;	  0 - Try to remove noise by applying filter to fourier spectra
;	  1 - Remove noise by subtracting median of Channel 0 from the rest
;	      of the array
;OUTPUTS:
;
;	  TvFlag = 2, 
;	   A postscript file will be written with all of the Fourier Plots
;	  TvFlag = 3
;	   A .png file will be written with the fourier spectrum for channel 0
;CALLING SEQUENCE:
;
;
pro Examine_Ref_Pixels, FitsFiles, Mode = Mode, RootDir = RootDir, Date=Date, $
	ReducedDir = ReducedDir, $
	ObjectName=ObjectName, NReads = NReads, $
	ReadNum = ReadNum, ThisFilter = ThisFilter,$
	OverRide = OverRide, TvFlag = TvFlag, AvgEntireFrame = AvgEntireFrame, $
	ExamineEachChannel = ExamineEachChannel

  Common KeyParams, KeyStr
  Common ColorStr, Color, ColorNum, ColorTri, ColorNumTri

  ;Get all of the defaults for the keywords from the KeywordStruct file
  @KeywordStruct.pro
  @PlotSettings.pro
  File_Mkdir, FourierDir

  If N_Elements(Mode) eq 0 then Mode = 'Channels'
  If stregex(FitsFiles, '.lst', /boolean) then Plot_AllNights, FitsFiles
  If Not keyword_set(AvgEntireFrame) then AvgEntireFrame = 0
  If Not keyword_set(ReadNum) then ReadNum = 0
  If N_Elements(FilStr) eq 0 then FilStr = ''
  If N_Elements(ExcludeCh0) eq 0 then ExcludeCh0 = 1
  If N_Elements(ExamineEachChannel) eq 0 then ExamineEachChannel = 0
 
  ;POSTSCRIPT SET-UP**********************************************************
  If TvFlag eq 2 then begin
    device,filename=FourierDir+$
		    KeyStr.RawKey+'_FourierSpectra.ps'
  EndIf 

  ;**************************************************************************
  ;Get a few kewyords
  FileParams = Return_File_Params(FitsFiles)
 
  ;Read the image that will be worked on
  ReadNum = 0	
  Im = Fits_Read_DataCube_Region(FitsFiles, FCB, $
				 FitsHeader, 0, KeyStr.Naxis1-1, 0, $
				 KeyStr.Naxis2-1, $
				 ReadNum, ReadNum)
  Case Mode of
    'Channels': begin
     ;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))

     ;Ch0 Analysis **********************************************************
     FFTChannel0 = shift(FFT(Im(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
       stop
       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 = Im(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
      stop
      erase
     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

     ;Individual Channels (0:127) *******************************
     ;Channel 0 
     Ch0 = Im(0:KeyStr.Ch0W-1, *)
     Ch0Med = Median(Ch0,dimension = 1) - Mean(Ch0)

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

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

     Masks = Return_Pixel_Masks()
     ImTemp = Im
     ImTemp(Masks.BadPixels)=0

     window,xsize=800,ysize=800
     !p.multi = [0,2,2]
     ;Plot the other channels
      For ChN = 0, 3 do begin
       Ch = Im(KeyStr.Ch0W*ChN:KeyStr.Ch0W*(ChN+1)-1,*)
       ChNoBad = ImTemp(KeyStr.Ch0W*ChN:KeyStr.Ch0W*(ChN+1)-1,*)
       FFTChannel = shift(FFT(Im(KeyStr.Ch0W*ChN:KeyStr.Ch0W*(ChN+1)-1,*),$
	Dimension = 2), 0 , KeyStr.Naxis2/2)
       FFTChannelAvg = Avg(abs(FFTChannel),0)
       print, 'Coefficients of Power for Channel ' + Strtrim(ChN,2)
       print, FFTChannelAvg(2049:2069)
       Plot,CFreqs(2049:2148), FFTChannelAvg(2049:2148)
     EndFor
       Plot,CFreqs(2049:2148), FFTChannelAvg(2049:2148)*FFTChannel0Avg(2049:2148)
      stop

      ;Ch0 Analysis **********************************************************
      FFTChannel = shift(FFT(Im(0:KeyStr.Ch0W-1,*),Dimension = 2), $
        0 , KeyStr.Naxis2/2)
      FFTChannelAvg = Avg(abs(FFTChannel),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
    End
    'Filter': begin
    ;Prepare a cube for the filtered image, Filter contains the coefficients 
    ;for each frequency in F(s) = FFT(f(t))
  
     ImFiltered = fltarr(KeyStr.Naxis1,KeyStr.Naxis2)
     Filter = RefFilter(Im, AvgEntireFrame = AvgEntireFrame, $
	   FilStr = '', ExamineEachChannel = ExamineEachChannel)
  
     ;Filter either with Fourier Coefficients or Median Subtraction
     Case RemovalMode of 
      0 : Begin
      ;Filter the Image by cutting out the power at various frequencies
      for Col = 0, Naxis2 - 1 do begin
        FFTImCol  = fft(Im(Col, *))
        FFTImCol  = Filter*FFTImCol
        ImFiltered(Col,*) = fft(FFTImCol,/inverse)  
      endfor      
      Filtered = reffilter(ImFiltered, $
	AvgEntireFrame = AvgEntireFrame, FilStr = 'FFTFiltered')
      SxAddPar, FitsHeader, 'REFTYPE', 'Fourier Filter'
      FilStr = 'FourierFiltered'
      End
      1 : Begin 
      ;Filter the Image by subtracting the median of channel 0
      ImFiltered = subtract_60hz(Im)
      Filtered = reffilter(ImFiltered, $
        AvgEntireFrame = AvgEntireFrame, FilStr = '60HzRemvoed_')
      SxAddPar, FitsHeader, 'REFTYPE', 'Median Channel 0 Subtracted'
      FilStr = '60HzRemoved'
     End
     EndCase
    End
  EndCase
  ;Close the Postscript device if open
  If TvFlag eq 2 then device, /close

  ;Write the Filtered Image
  Fits_Write, KeyStr.ReducedDir+KeyStr.RawKey+'_'+FilStr+'.fits', $
	      ImFiltered, FitsHeader
  
  ;Default back to the X Window
  set_plot,'x'
stop

end
