/************************************************************************\

				_ALPHA_H2O_

	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.

\************************************************************************/
#ifndef _ALPHA_H2O_
#define _ALPHA_H2O_	"1998 October"
#include	<math.h>
static double dataH2O[5][3][17] =
	{ { {
	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	
	} , {
	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
	} , {
	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
	} } , { {
	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
	} , {
	-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
	} , {
	-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
	} } , { {
	-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
	} , {
	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
	} , {
	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
	} } , { {
	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
	} , {
	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
	} , {
	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
	} } , { {
	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
	} , {
	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
	} , {
	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
	} } } ;

/************************************************************************\
	Gaussian function for a given frequency nu (in inverse cm),
	central frequency cwav, full width at half maximum fwhm, and
	peak amplitude peak.
\************************************************************************/
double	fn_gaussian(double nu, double cwav, double fwhm, double peak)
{
    double	c1,
		c2=-2.77258872224,
		c3,
		result=0.0;
    c1 = nu - cwav;
    c1 = c1*c1;
    c3 = c1*c2/(fwhm*fwhm);
    if(c3 > -40.) result = peak*exp(c3);
    return(result);
}	/****************************************************************/
/************************************************************************\
	Return an array of absorption coefficients at the wavelengths
	provided in the array *lambda, in units of microns.  The number
	of wavelengths is given by n, and the temperature in Kelvins is
	given by T.  Values of T outside the range 20 <= T <= 270 are
	rejected and wavelengths outside the range 1 to 2.7 microns
	return zeros.
	Memory space for *y needs to have already been allocated.
\************************************************************************/
void	*alphaH2O(double *lambda, double *y, int n, double T)
{
    int		i,j;
    double	T2,T3,T4,
		nu,cwav,fwhm,peak;
    if(T > 270.0 || T < 20) {
	printf("ERROR: T outside of valid range.\n");
	exit(1);
    }
    T2 = T*T;
    T3 = T2*T;
    T4 = T3*T;
    for(i=0; i<n; i++) {
	nu = 10000.0/lambda[i];			/* convert to frequency */
	y[i] = 0.0;
	if(lambda[i] > 0.97 && lambda[i] < 2.75) for(j=0; j<17; j++) {
	    cwav = dataH2O[0][0][j] + T*dataH2O[1][0][j] +
		   T2*dataH2O[2][0][j] + T3*dataH2O[3][0][j] +
			T4*dataH2O[4][0][j];
	    fwhm = dataH2O[0][1][j] + T*dataH2O[1][1][j] +
		   T2*dataH2O[2][1][j] + T3*dataH2O[3][1][j] +
			T4*dataH2O[4][1][j];
	    peak = dataH2O[0][2][j] + T*dataH2O[1][2][j] +
		   T2*dataH2O[2][2][j] + T3*dataH2O[3][2][j] +
			T4*dataH2O[4][2][j];
	    if(peak > 0.0)
		y[i] += fn_gaussian(nu,cwav,fwhm,peak);
	}
    }
}	/****************************************************************/
#endif	/* _ALPHA_H2O_ */