/*hmiopc derive total H- H minus opacity */
#include "cddefines.h"
#include "spline.h"
#include "splint.h"
#include "hmiopc.h"
#	ifdef NCRS
#undef NCRS
#	endif
#define	NCRS	33

double hmiopc(double freq)
{
	double energy, 
	  hmiopc_v, 
	  x, 
	  y;
	static double y2[NCRS];
	static double crs[NCRS]={0.,0.124,0.398,0.708,1.054,1.437,1.805,
	  2.176,2.518,2.842,3.126,3.377,3.580,3.741,3.851,3.913,3.925,
	  3.887,3.805,3.676,3.511,3.306,3.071,2.810,2.523,2.219,1.898,
	  1.567,1.233,.912,.629,.39,.19};
	static double ener[NCRS]={0.,0.001459,0.003296,0.005256,0.007351,
	  0.009595,0.01201,0.01460,0.01741,0.02044,0.02375,0.02735,0.03129,
	  0.03563,0.04043,0.04576,0.05171,0.05841,0.06601,0.07469,0.08470,
	  0.09638,0.1102,0.1268,0.1470,0.1723,0.2049,0.2483,0.3090,0.4001,
	  0.5520,0.8557,1.7669};
	static int lgFirst = TRUE;

#	ifdef DEBUG_FUN
	fputs( "<+>hmiopc()\n", debug_fp );
#	endif

	/* bound free cross section (x10**-17 cm^2) from Doughty et al
	 * 1966, MNRAS 132, 255; good agreement with Wishart MNRAS 187, 59p. */

	/* photoelectron energy, add 0.05552 to get incoming energy (Ryd) */


	if( lgFirst )
	{
		/* set up coefficients for spline */
		spline(ener,crs,NCRS,2e31,2e31,y2);
		lgFirst = FALSE;
	}

	energy = freq - 0.05552;
	if( energy < ener[0] || energy > ener[NCRS-1] )
	{
		hmiopc_v = 0.;
	}
	else
	{
		x = energy;
		splint(ener,crs,y2,NCRS,x,&y);
		hmiopc_v = y*1e-17;
	}


#	ifdef DEBUG_FUN
	fputs( " <->hmiopc()\n", debug_fp );
#	endif
	return( hmiopc_v );
}