;+
; NAME:
;	ALPHA_H2O
;
; PURPOSE:
;	Generate synthetic absorption coefficient spectrum of H2O ice.
;
; DESCRIPTION:
;	Subroutine for generating synthetic absorption coefficient
;	spectra of H2O ice for wavelengths between 1 and 2.7 microns
;	and temperatures between 20 and 270 K.
;
;	These data were 1st published in ``The Temperature-Dependent
;	Near-Infrared Absorption Spectrum of Hexagonal H2O Ice'' by
;	W.M. Grundy and B. Schmitt (1998) J. Geophys. Res. E 103,
;	25,809-25,822.
;
;	If you use this subroutine, please cite the JGR paper.
;	Please also read the paper to understand the limitations
;	of the data set and this subroutine.
;
; CATEGORY:
;	Ices
;
; CALLING SEQUENCE:
;	Y = ALPHA_H2O(T,X)
;
; INPUTS:
;	T:	Scalar temperature (in Kelvins).
;	X:	Array of wavelengths (in microns).
;
; OPTIONAL INPUT PARAMETERS:
;	None.
;
; KEYWORD INPUT PARAMETERS:
;	None.
;
; OUTPUTS:
;	Y:	The absorption coefficient (in inverse cm) at each X.
;
; KEYWORD OUTPUT PARAMETERS:
;	None.
;
; OPTIONAL OUTPUT PARAMETERS:
;	None.
;
; COMMON BLOCKS:
;	None.
;
; SIDE EFFECTS:
;	None.
;
; RESTRICTIONS:
;	None.
;
; PROCEDURE:
;	Return the sum of 17 smoothly-varying Gaussians.
;
; MODIFICATION HISTORY:
;	1997 December 2 - created by W.M. Grundy
;-
PRO ICE_GAUSSIAN,X,A,F
    ON_ERROR,2			; Return to caller if an error occurs
    C1 = (X-A[0])^2
    C2 = -2.77258872224		; -4ln(2)
    C3 = -40.			; cut-off for exponential
    if A[1] ne 0.0 then $	; Avoid possible division by zero
	F = A[2]*EXP((C2*C1/(A[1]*A[1])) > C3) else F = 0.0
END

FUNCTION ALPHA_H2O,T,x
    IF n_elements(T) GT 1 THEN BEGIN
	print,'T cannot be an array in this context'
	return,0.0
    ENDIF
    IF T GT 270.0 THEN return,0.0
    IF T LT 20.0 THEN return,0.0
    data = double(bytarr(17,3,5))
    data[0:16,0,0] = [ $
	9483.8818,	9253.9001,	7798.8160,	7615.8508, $
	7200.9562,	6663.2139,	6564.3830,	6380.0270, $
	6319.0626,	6028.9655,	5596.0577,	4952.7117, $
	4955.9127,	4823.6550,	3960.4870,	4241.2252, $
	3381.9357]
    data[0:16,1,0] = [ $
	842.98265,	-2515.6917,	743.49830,	101.10864, $
	357.60369,	128.41358,	334.90895,	113.07332, $
	627.37620,	84.087991,	337.22006,	259.04133, $
	435.93215,	143.63656,	538.18934,	154.57263, $
	256.61140]
    data[0:16,2,0] = [ $
	0.34067697,	-0.94738632,	1.5301105,	0.60498747, $
	0.59147538,	17.114210,	23.949553,	13.952447, $
	21.389032,	30.884946,	3.3610640,	45.281503, $
	48.260185,	28.314718,	26.601682,	8.3508344, $
	2519.5298]
    data[0:16,0,1] = [ $
	1.1939951,	-3.8718575,	0.16506409,	0.38520589, $
	0.21724120,	-0.011057193,	-0.018193390,	-0.27553639, $
	2.9884802,	0.23811586,	-0.84826079,	0.16109639, $
	-0.12277913,	-0.42984049,	0.045487661,	0.039077274, $
	0.29040582]
    data[0:16,1,1] = [ $
	-0.64530980,	23.369389,	0.58789111,	0.10681935, $
	-1.5862691,	0.26451368,	0.46782138,	-0.021004515, $
	1.2247057,	0.13492106,	-0.030634305,	0.25331963, $
	-0.095685670,	-0.19504389,	0.18529675,	0.11303142, $
	0.0029338089]
    data[0:16,2,1] = [ $
	-0.00044105802,	0.0082040029,	-0.0028639194,	0.0019415733, $
	-0.0016811397,	-0.0024274845,	0.048208950,	-0.021326207, $
	-0.013900167,	0.020044197,	-0.0070897602,	0.038146955, $
	0.061494533,	0.027264885,	0.017260260,	-0.010686407, $
	-15.293000]
    data[0:16,0,2] = [ $
	-0.0013815159,	0.010880674,	0.0016978521,	-0.0011801064, $
	0.0018259098,	0.0011561131,	0.00077766615,	0.0029255075, $
	-0.055533051,	0.00028000425,	0.0027371178,	0.0014876875, $
	0.00016533177,	0.0031159540,	-0.00020811705,	-0.00031365548, $
	0.00066707098]
    data[0:16,1,2] = [ $
	0.0,		-0.046524376,	-0.00061644252,	0.00070551001, $
	0.016101357,	0.00023187442,	-0.00013411424,	0.0068366000, $
	0.00056016259,	0.00011142521,	-0.012910437,	-0.0018331047, $
	0.0010062839,	0.0023713817,	0.00089175208,	-0.00082327650, $
	0.0]
    data[0:16,2,2] = [ $
	0.0,		-1.6269149e-05,	1.0158347e-05,	-8.2594070e-05, $
	-2.9964421e-05, -3.6366681e-05,	-0.00018831293,	-0.00029404347, $
	-0.00079668155,	-0.0033033951,	-3.1383799e-06,	-0.00013005703, $
	-0.00024552345, -0.00056810739,	7.3466150e-05,	-4.0785765e-05, $
	0.028044127]
    data[0:16,0,3] = [ $
	0.0,	0.0,	0.0,		0.0,	0.0,	0.0, $
	0.0,	0.0,	0.00027835949,	0.0,	0.0,	0.0, $
	0.0,	0.0,	0.0,		0.0,	0.0]
    data[0:16,1,3] = [ $
	0.0,		0.0,	0.0,		0.0,	-3.6541833e-05, $
	0.0,		0.0,	-2.0374012e-05,	0.0,	0.0, $
	9.1960954e-05,	0.0,	0.0,		0.0,	0.0, $
	0.0,		0.0]
    data[0:16,2,3] = [ $
	0.0,		0.0,		0.0,		4.4933288e-07, $
	3.1187993e-07,	0.0,		0.0,		7.2766211e-07, $
	6.0841622e-06,	1.8298112e-05,	1.4517029e-07,	0.0, $
	0.0,		2.1698425e-06,	0.0,		0.0, $
	0.0]
    data[0:16,0,4] = [ $
	0.0,	0.0,		0.0,	0.0,	0.0,	0.0, $
	0.0,	0.0, -4.4312719e-07,	0.0,	0.0,	0.0, $
	0.0,	0.0,		0.0,	0.0,	0.0]
    data[0:16,1,4] = [ $
	0.0,	0.0,	0.0,	0.0,		0.0,	0.0, $
	0.0,	0.0,	0.0,	0.0, -1.7517519e-07,	0.0, $
	0.0,	0.0,	0.0,	0.0,		0.0]
    data[0:16,2,4] = [ $
	0.0,		0.0,		0.0,		-7.4082871e-10, $
	-6.4752963e-10,	0.0,		0.0,		0.0, $
	-1.2341960e-08,	-2.9321373e-08, -3.4170770e-10,	0.0, $
	0.0,		0.0,		0.0,		0.0, $
	0.0]
    nu = 10000.0/x		; convert to wavenumber
    Y = double(bytarr(n_elements(x)))
    F = double(bytarr(n_elements(x)))
    T = double(T)
    T2 = T*T			; Powers of temperature
    T3 = T2*T
    T4 = T3*T
    FOR i=0,16 DO BEGIN		; sum contributions of 17 Gaussians
	cwav = data[i,0,0] + T*data[i,0,1] + T2*data[i,0,2] + $
	    T3*data[i,0,3] + T4*data[i,0,4]
	fwhm = data[i,1,0] + T*data[i,1,1] + T2*data[i,1,2] + $
	    T3*data[i,1,3] + T4*data[i,1,4]
	peak = data[i,2,0] + T*data[i,2,1] + T2*data[i,2,2] + $
	    T3*data[i,2,3] + T4*data[i,2,4]
	A = [cwav,fwhm,peak > 0]
	ice_gaussian,nu,A,F
	Y = Y + F
    ENDFOR
    nel = n_elements(x) - 1	; Zero outside of wavelength range
    FOR i=0,nel DO BEGIN
	IF x[i] LT 0.97 THEN Y[i] = 0.0
	IF x[i] GT 2.75 THEN Y[i] = 0.0
    ENDFOR
    return,Y
END