PRO wave_grid, WavMin, Len, PixSize, OverSampl, TG_Period, X_Rowland, $
              MaxOrd, WavCtr, WavBinSize
;+
; NAME:
;         wave_grid
;
;
; PURPOSE:
;         Create arrays of central wavelengths and bin widths in which
;         the sampling follows the optimal binning for AXAF
;         transmission grating spectra to accurately reflect higher
;         resolution in higher orders.
;
;
; CATEGORY:
;         AXAF transmission grating model spectra
;
;
; CALLING SEQUENCE:
;         wave_grid, wavmin, Len, PixSize, OverSampl, TG_Period, X_Rowland, $
;              MaxOrd, WavCtr, WavBinSize
;      
;
; 
; INPUTS:
; 
;        WavMin (R)   minimum wavelength for output array, in Angstroms.
;        Len    (R)   Approximate length of detector array, from
;                     zero-order to outer limit, in mm. 
;        PixSize  (R) Pixel size of the array in mm (not
;                     necessarily a resolution element.
;        OverSampl (R) Amount of oversampling of pixels of size
;                      PixSize (e.g., if dithering gives half-pixel
;                      resolution, then specify OverSampl of 2); value
;                      is factor to divide PixSize by.
;        TG_Period (R) Period of the diffraction grating, in Angstroms.
;        X_Rowland (R) Rowland diameter of the grating in mm.
;                      (canonically 8635mm, but larger at XRCF)
;        MaxOrd  (I)   Maximum order to consider.
;
;
;
; OPTIONAL INPUTS:
;                   NONE.
;
;	
; KEYWORD PARAMETERS:
;                   NONE.
;
;
; OUTPUTS:
;        WavCtr     (R)  Array of bin central wavelengths (in Angstroms).
;        WavBinSize (R) Array of bin full widths, corresponding to the
;                       WavCtr bins (in Angstroms).
;
; OPTIONAL OUTPUTS:
;         NONE
;
;
; COMMON BLOCKS:
;          NONE.
;
;
; SIDE EFFECTS:
;          NONE.
;
;
; RESTRICTIONS:
;          
;
;
; PROCEDURE:
;        In order to accurately model the higher diffraction orders of
;        a predicted spectrum, the input model spectrum must have
;        sufficient resolution to be sampled by the detector.  A grid
;        which matches the first order spectrum will be undersampled
;        in higher orders, since higher orders are increasingly
;        dispersed by a factor of m, the diffraction order.  If the
;        grid were set to match the highest order to be considered,
;        then the input spectrum would have a very large number of
;        grid points.  This is unnecessary, since much of the spectrum
;        is diffracted distances greater than the length of the
;        detector.  This procedure constructs a grid which has
;        increasingly finer resolution at progressively shorter
;        wavelengths.  The grid spacing is discontinuous at boundaries
;        where the next higher order will fall off the detector.
;
;        Gridding is as follows:
;           From Len to Len/2, grid at the first order resolution
;           (essentially the detector pixel size, possibly oversampled
;           a bit).
;
;           From Len/2 to Len/3, grid at half the first order.
;
;           From Len/3 to Len/4, grid at one-third the first order.
;
;           ... and so forth.
;
;        The general form is:
;            n_m = the number of grid points for order m,
;            n = the number of grid points required to span from
;                wavelength=0 (zero-order centroid) to Len at the
;                first order resolution. 
;            n_m = [ 1/m - 1/(m+1) ] * m  * n  =   n / (m+1)
;
;        Thus the total number of grid-points required is:
;            N = sum_{m=1..MaxOrd} { [1/m - 1/(m+1)] * m  * n}
;              = n * sum_{m=1..MaxOrd} {1 / (m+1)}
;              = n * { Psi(2+MaxOrd) - 1 + gamma},
;
;         where Psi(x) is the digamma function, and gamma is Euler's
;         constant.
;
;         This is a diverging but slowly increasing function.  For
;         MaxOrd of 10, 20, and 30, N is 2.0*n, 2.6*n, and 3.0*n,
;         respectively.  (If the spectrum were gridded at the highest
;         resolution required, N would scale directly with MaxOrd.)
;
;         Given some reasonable MaxOrd (such as 30 for LETGS), the
;         minimum wavelength grid point calculated as outlined above
;         will still be above the minimum needed in model spectra.
;         Thus, we have added a minimum wavelength as a parameter, and
;         the grid will be filled from the minimum grid point to the
;         specified minimum at the gridding of the MaxOrd segment.
;
;
; EXAMPLE:
;         For a typical HRC-S/LETG grid:
;
;           WavMin = 1.2398 ; 10 keV
;           Len = 160.0;   mm, approx length of one side of HRC-S
;           PixSize = 0.0064 ; mm, set by electronic readout
;           OverSampl = 1.0 ;  HRC readout already oversamples the
;                              FWHM of the charge cloud (0.025 mm)
;           TG_Period = 9921 ; Angstroms; LETG baseline period.
;           X_Rowland = 8635 ; mm, Rowland diameter, flight config.
;           MaxOrd = 30 ;      should be plenty; strong lines in 27th
;                              order can contribute a few percent of
;                              counts at some positions.
;
;           wave_grid, WavMin, Len, PixSize, OverSampl, TG_Period, X_Rowland, $
;              MaxOrd, WavCtr, WavBinSize
;                   
;           help,wavctr,wavbinsize; see how big the arrays are
;           plot,wavctr           ; examine the wavelength scale.
;           plot,wavbinsize       ; examine the grid size.
;
;
; MODIFICATION HISTORY:
;           Version 1.0  2 december 1995, d.p. huenemoerder
;
;-

   
   IF n_params() NE 9 THEN BEGIN
      print,''
      print,'wave_grid: incorrect number of parameters.'
      print,''
      print,'USAGE: wave_grid, WavMin, Len, PixSize, OverSampl, ' + $
       'TG_Period, X_Rowland, MaxOrd, WavCtr, WavBinSize'
      print,''
      print,'    WavMin (R)   minimum wavelength for output array, in Angstroms.'
      print,'    Len    (R)   Approximate length of detector array, from'
      print,'                 zero-order to outer limit, in mm. '
      print,'    PixSize  (R) Pixel size of the array in mm (not'
      print,'                 necessarily a resolution element.'
      print,'    OverSampl (R) Amount of oversampling of pixels of size'
      print,'                  PixSize (e.g., if dithering gives half-pixel'
      print,'                  resolution, then specify OverSampl of 2); value'
      print,'                  is factor to divide PixSize by.'
      print,'    TG_Period (R) Period of the diffraction grating, in Angstroms.'
      print,'    X_Rowland (R) Rowland diameter of the grating in mm.'
      print,'                  (canonically 8635mm, but larger at XRCF)'
      print,'    MaxOrd  (I)   Maximum order to consider.'
      print,''
      print,'    OUTPUTS:'
      print,'     WavCtr     (R)  Array of bin central wavelengths (in Angstroms).'
      print,'     WavBinSize (R) Array of bin full widths, corresponding to the'
      print,'                    WavCtr bins (in Angstroms).'
      print,''

      return
   ENDIF
   

;; determine the number of grid points for full first order:
;
   n_first_order = fix(Len / PixSize * OverSampl)
   

;; determine maximum wavelength - wavelength at Len:
;
   wave_max = TG_Period * Len / X_Rowland


; Construct the first-order part of the grid, from Len/2 (maximum
; wavelength/2) to Len (maximum wavelength)
   
   
   this_m = 1
   
   wave_min =  wave_max / (this_m + 1)
   delta_wave = wave_max / (this_m * (this_m + 1))

   n_this_m = n_first_order / (this_m + 1)
   grid_center = findgen(n_this_m)/n_this_m * delta_wave + wave_min
   grid_bins = replicate(delta_wave/n_this_m, n_this_m)


;; form and concatenate grids up to highest order specified
;
   FOR this_m = 2, MaxOrd DO BEGIN
      
      wave_min =  wave_max / (this_m + 1)
      delta_wave = wave_max / (this_m * (this_m + 1))
      n_this_m = n_first_order / (this_m + 1)

      grid_center = [findgen(n_this_m)/n_this_m * delta_wave + wave_min, grid_center]
      grid_bins = [replicate(delta_wave/n_this_m, n_this_m), grid_bins]
      
   ENDFOR
   

;; extend down to the minimum wavelength requested
;
   IF WavMin LT wave_min THEN BEGIN
      
      delta_wave = wave_min-WavMin
      n_left = delta_wave / grid_bins(0)
      
      grid_center = [findgen(n_left)/n_left * delta_wave + WavMin, $
                     grid_center]
      grid_bins = [replicate(grid_bins(0), n_left), grid_bins]
      
   ENDIF
   
   
   WavCtr = grid_center
   WavBinSize = grid_bins
   
END