;+
; Project     : HESI
;                   
; Name        : grid_diffraction
;               
; Purpose     : This function returns the energy as a function of
;	transmission amplitude due to diffraction effects through the
;	HESI grids. 
;               
; Category    : HESI 
;               
; Explanation :  	diffraction_trans = 1. - (reduction in mod. amplitute/100.) 
;		en = 12.39e-7*pi*grid_sep*harmonic_number^2 / 
;      		ARCOSINE(diffraction_trans)*4*((.5*pitch)^2)
;
;               
; Use         :  result = grid_diffraction(fraction, harmonic_number, $
;		pitch, grid_sep)
;    
; Inputs      : 	Fraction - Find this energy at this fraction of the transmission curve.
;		Valid inputs range from 0-1.
;		Input may be a vector.
;		Nharm - Number of harmonics to solve for, starting at 1.
;		Pitch - grid pitch in mm.  Pitch is spacing between slit centers
;		Input may be a vector.
;		Grid_sep - distance between collimator grid planes, mm.
;               
; Opt. Inputs : None
;               
; Outputs     : The result is a 3-d array dimensioned 
;		num_pitch x nharm x num_fraction.
;
; Opt. Outputs: None
;               
; Keywords    : ELIM- Energies below this value are set to 0.0
;		defaults to 0.1 keV
;
; Calls	      :
;
; Common      : None
;               
; Restrictions: 
;               
; Side effects: None.
;               
; Prev. Hist  : Based on grid diffraction formula found in many Hesi Menu 
;	procedures written by Eric Carzon, 1993.
;
; Modified    : Encoded into a function by RAS, 16-apr-1997
;
;-            
;==============================================================================

function grid_diffraction, fraction, nharm, pitch, grid_sep, elim=elim

nfrac=n_elements(fraction)
npitch=n_elements(pitch)
checkvar, elim, 0.1

result = fltarr( npitch, nharm, nfrac)
fr = ((rebin( reform(fraction,1,1,nfrac), npitch, nharm, nfrac)>0.0)<1.0)
hm = rebin(indgen(1,nharm,1)+1, npitch, nharm, nfrac) 
pc = rebin( reform(pitch,npitch,1,1), npitch, nharm, nfrac)

energy= (12.39854e-7*!pi*grid_sep*((hm)^2))/(acos(fr)*4*((.5*pc)^2))
w = where( energy le elim, nw)
if nw gt 0 then energy(w)=0.0

return, energy
end