#pragma rtGlobals=3		// Use modern global access method and strict wave access.
#include <KBColorizeTraces>


//---------------------------------------------------------------------------------------------

FUNCTION /WAVE TSCal_HistogramLayers(WaveIn, [Disp, OutName])
  Wave WaveIn
  Variable Disp
  String OutName
  
  //NAME: TSCal_HistogramLayers
  //
  //PURPOSE: 
  //  Calculate histograms for each pixel given a multilayer flat-field
  //  calibration image.
  //  Will output a 4D wave where rows & columns correspond to pixel location,
  //     and layer data corresponds to the histogram data for that pixel
  //
  //AUTHOR:
  // Dave Schlossberg
  // Univ. of Wisconsin - Madison
  // 1500 Engineering Drive
  // Madison, WI 53706
  //
  //CATEGORY:
  // Thomson scattering calibration analyses
  //
  //INPUTS:
  // WaveIn		4D wave containing flat-field calibration data to be anlayzed
  //
  //OPTIONAL INPUTS:
  // Disp			Set nonzero to display histogram results
  // OutName		String containing desired name for output wave
  //
  //OUTPUTS:
  // HistoLayers		4D wave containing histograms along layer dimension.
  //						Rows & columns correspond to original row & column locations in original wave
  //
  //EXAMPLE CALLING SEQUENCES:
  //
  //DEPENDENCIES:
  // None.
  //
  //MODIFICATION HISTORY:
  // Created 6/19/16 by Dave Schlossberg
  //
  //
  //
  
  IF( ParamIsDefault(OutName) )
    OutName = UniqueName("HistoLayers", 1, 0)
  ENDIF
  // Make output wave in which to store histogram results
  Duplicate /O WaveIn $OutName
  Wave OutWave = $OutName
  // Insert a new "chunk" to hold the histogram scale information
  InsertPoints /M=3 1, 1, OutWave
  OutWave = 0
  
  // Make temporary waves for histogram of single beams (i.e. p0, q0, all layers)
  Make /O/N=(dimsize(WaveIn, 2))/FREE BeamVals, BeamHisto
  
  // Step through input wave assembling histograms and corresponding scales.
  Variable indy_row, indy_col
  FOR( indy_col=0; indy_col<DimSize(WaveIn, 1); indy_col+=1 )
    FOR( indy_row=0; indy_row<DimSize(WaveIn, 0); indy_row+=1 )
      BeamVals = WaveIn[indy_row][indy_col][p]
      Histogram /B=1 BeamVals, BeamHisto
      // Store histogram values in corresponding row & col in Outwave, Chunk 0
      OutWave[indy_row][indy_col][][0] = BeamHisto[r]
      // Store scale information in corresponding row & col in Outwave, Chunk 1
      OutWave[indy_row][indy_col][][1] = DimOffset(BeamHisto, 0) + DimDelta(BeamHisto,0)*r
    ENDFOR
  ENDFOR
  
  
  IF( !ParamIsDefault(Disp) && (Disp!=0) )
    TSCal_PlotHistogramLayers(OutWave)
  ENDIF
  
  RETURN OutWave
  
END

//------------------------------------------------------------------------------------------

FUNCTION TSCal_PlotHistogramLayers(HistoLayers, [Rows, Cols])
  WAVE HistoLayers, Rows, Cols
  
  // If user doesn't specify rows & columns for which to show historgrams
  //  simply show them all
  IF( ParamIsDefault(Rows) )
    Make /O/N=(dimsize(HistoLayers,0))/FREE Rows = p
  ENDIF
  IF( ParamIsDefault(Cols) )
    Make /O/N=(dimsize(HistoLayers,1))/FREE Cols = p
  ENDIF
  
  // Make a new window
  Display
  
  Variable indy_row, indy_col
  
  FOR( indy_col=0; indy_col<NumPnts(Cols); indy_col+=1 )
    FOR( indy_row=0; indy_row<NumPnts(Rows); indy_row+=1 )
      AppendToGraph HistoLayers[Rows[indy_row]][Cols[indy_col]][][0]/TN=$("Row"+Num2Str(Rows[indy_row])+"Col"+Num2Str(Cols[indy_col])) vs HistoLayers[Rows[indy_row]][Cols[indy_col]][][1]
    ENDFOR
  ENDFOR
  
  
END

//------------------------------------------------------------------------------------------

FUNCTION RunThisForExample()
  
  Make /O/N=(16,8,100) InputWave = enoise(20)
  TSCal_HistogramLayers(InputWave, Disp=1)
  CommonColorsButtonProc("")
  ModifyGraph hidetrace=1, hidetrace[0]=0, mode=5
  
END

//------------------------------------------------------------------------------------------