/*rayleh compute Rayleigh scattering cross section for Lya */
#include "cddefines.h"
#include "powi.h"
#include "physconst.h"
#include "dampon.h"
#include "rayleh.h"

double rayleh(double ener)
{
	double rayleh_v;

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

	/* do hydrogen Rayleigh scattering cross sections;
	 * fits to Phys Rev 163, 147; and Mihalas radiative damping
	 *
	 * >>chng 96 aug 15, changed logic to do more terms for each part of
	 * rayleigh scattering
	 * if( ener.lt.0.05 ) then
	 *  rayleh = 8.41e-25 * ener**4 * DampOnFac
	 * */
	if( ener < 0.05 )
	{
		rayleh_v = (8.41e-25*powi(ener,4) + 3.37e-24*powi(ener,6))*
		  dampon.DampOnFac;

	}
	else if( ener < 0.646 )
	{
		rayleh_v = (8.41e-25*powi(ener,4) + 3.37e-24*powi(ener,6) + 
		  4.71e-22*powi(ener,14))*dampon.DampOnFac;

	}
	else if( ener >= 0.646 && ener < 1.0 )
	{
		rayleh_v = fabs(0.74959-ener);
		rayleh_v = 1.788e5/POW2(FR1RYD*MAX2(0.001,rayleh_v));
		/*  typical energy between Ly-a and Ly-beta */
		rayleh_v = MAX2(rayleh_v,1e-24)*dampon.DampOnFac;
	}
	else
	{
		rayleh_v = 0.;
	}


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