#include <math.h>
#include "units.h"
#include "BlackB.h"

// ctor: Kelvin
BlackB::BlackB(double temp) : T(temp)
{
}

/* Watts per square meter and delta lambda (lambda in Meter)
2*pi*h c^2/lambda^5/[exp{h c/(lambda k T)}-1]
h c^2=5.95521e-17 = 3.74177e-16/(2pi) (Watts/m^2)
Divide by pi to obtain the value per srad.
h*c/k_B about 0.01438768 meter times Kelvin.
Integrated over all lambda: sigma*T^4, where sigma=Stefan-Boltzman's constant
********************************************************************/
double BlackB::densitlam(const double lambda) const
{
	return 2.*M_PI*PLANCK*VACUUMC*VACUUMC/pow(lambda,5.)/(exp(PLANCK*VACUUMC/BOLTZMANNK/(lambda*T))-1.) ;
}

/********* Watts per square meter and delta nu with nu in Hz, therefore units of return value J/m^2
 2*pi*h nu^3/c^2/[exp(h nu/(k T))-1]
 with c=lamda*nu, lambda=c/nu, d lambda= -c/nu^2 d nu back to the
 formula quoted above.
 h /c^2=7.3732e-51
**********************************************************************/
double BlackB::densitnu(const double nu) const
{
	return 2.*M_PI*PLANCK/(VACUUMC*VACUUMC)*pow(nu,3.)/(exp(PLANCK/BOLTZMANNK*nu/T)-1.) ;
} 

// integral over densitnu between these limits (interpreted as wavenumbers if argis_freq=false,
// else as frequency in Hz). Units of return value W/m^2.
// The unit of h nu^3/c^2 is Ws/m^2 and needs to be integrated over 1/s=Hz to get W/m^2
double BlackB::fluxnu(const double nu[2], const bool argis_freq) const
{
	if( argis_freq)
	{
		// integral substitution factor
		const double f = PLANCK/(BOLTZMANNK*T) ;
		return 2.*M_PI*PLANCK/(VACUUMC*VACUUMC)/pow(f,4.)*( debye3(f*nu[1])-debye3(f*nu[0]) ) ;
	}
	else
	{
#if 0
		// 2pi h/c^2 * nu^3/[exp(h nu/k T)-1] dnu = 2*pi*10^8 c^2 nutilde^3/[exp(100 c h nutilde/k/T)-1] d nutilde
		// where the conversion factor is HZ2WN and WN2HZ.
		// = 2pi (k*T)^4/(h^3*c^2) * x^3/[exp(x)-1] dx with x=h*100*VACUUMC*nutilde/(k*T)
		// integral substitution factor
		const double f = 100.*VACUUMC*PLANCK/BOLTZMANNK/T ;
		return 2.e8*M_PI*PLANCK*VACUUMC*VACUUMC/pow(f,4.)*( debye3(f*nu[1])-debye3(f*nu[0]) ) ;
#else
		double frequ[2] ;
		frequ[0] = WN2HZ(nu[0]) ;
		frequ[1] = WN2HZ(nu[1]) ;
		return fluxnu(frequ,true) ;
#endif
	}
} 

// Relative integral over densitnu between these limits (interpreted as wavenumbers or frequencies)
// assuming that the value between nuRef[2] is set to 1.
double BlackB::fluxnuRel(const double nu[2], const double nuRef[2]) const
{
	return fluxnu(nu,true)/fluxnu(nuRef,true) ;
} 

// moment of power 'n' over densitnu between these limits (interpreted as wavenumbers if argis_freq=false)
double BlackB::moment(const double nu[2], const bool argis_freq, const int n) const
{
	if( argis_freq)
	{
		// integral substitution factor
		const double f = PLANCK/(BOLTZMANNK*T) ;
		return 2.*M_PI*PLANCK/(VACUUMC*VACUUMC)/pow(f,3+n)*( debyen(f*nu[1],3+n)-debyen(f*nu[0],3+n) ) ;
	}
	else
	{
		double frequ[2] ;
		frequ[0] = WN2HZ(nu[0]) ;
		frequ[1] = WN2HZ(nu[1]) ;
		return moment(frequ,true,n) ;
	}
} 

		/* 1/(2pi) */
#define M_1_2PI .159154943091895335768883763373
		/* Riemann Zeta of 4 */
#define zeta_4 1.08232323371113819151600369654
#define K 20
double BlackB::debye3(const double x) const
{
	int k ;
	double xk,
	       sum, sumold ;
			/* list of Bernoulli numbers of index 2*k, multiplied by
			(2*pi)^(2k)/(2k+3)/(2k)!, k=0,1,2,3,4,... */
	static double koeff[K]= { 0., 
			.657973626739290574588966066660, -.309235209631753769004572484726,
			.226076235996544253269892873286, -.182559519308717152614306407002,
			.153999165404279705436483993677, -.133366144873774406439818399740,
			.117654264486477494685795122954, -.105264766553621963354919218052,
		       .0952384587898347618895101392038, -.0869566046923507715474881001774,
		       .0800000190760402182186392002920, -.0740740784894954852784799712918,
		       .0689655182690727467837959472180, -.0645161292726021951476423906336,
		       .0606060606625046928739192838104, -.0571428571561617819067228885258,
		       .0540540540572004173448122707508, -.0512820512827975344559508709866,
		       .0487804878050555111974331049360 } ;

	if (x<0.)
		return -1. ;
	sum=0. ;
			/* for values up to 1.80 the list of zeta functions
			and the sum up to k < K are huge enough to gain
			numeric stability in the sum */
	if(x>1.8)
	{
		for(k=1;;k++)
		{
			sumold=sum ;
				/* do not use x(k+1)=xk+x to avoid loss of precision */
			xk=x*k ;
			sum += exp(-xk)*(((xk+3.)*xk +6.)*xk+ 6.)/(k*k*k*k) ;
			if( sum == sumold)
				break ;
		}
				/* 6*zeta_4 is the integral fro x=infinity*/
		return 6.*zeta_4-sum ;
	}
	else
	{
		register double x2pi=x*M_1_2PI ;
		for(k=1;k<K;k++)
		{
			sumold=sum ;
				/* do not precompute x2pi^2 do avoid loss of precision */
			sum += koeff[k]*pow(x2pi,2.*k) ;
			if(sum == sumold)
				break ;
		}
		sum += 1./3. -0.125*x ;
		return sum*pow(x,3.) ;
	}
}
// Debye function of order n, int(t=0..x) t^n dt / [exp(t)-1]
double BlackB::debyen(const double x, const int n) const
{
	int k ;
	double xk,
	       sum=0., sumold ;
			/* list of Bernoulli numbers of index 2*k, multiplied  by
			(2*pi)^k/(2k)!, k=0,1,2,3,4
			Digits :=30 ; for k from 1 to 20 do; evalf( (2*Pi)^(2*k)*bernoulli(2*k)/(2*k)! ) ; od:
			*/
	static double koeff[K]= { 0., 
		3.28986813369645287294483033329, -2.16464646742227638303200739308,
		2.03468612396889827942903585958, -2.00815471239588867875737047702,
		2.00198915025563617067429191780, -2.00049217310661609659727599610,
		2.00012249627011740965851709022, -2.00003056451881730374346514298,
		2.00000763458652999967971292330, -2.00000190792406774559222630408,
		2.00000047690100545546598000730, -2.00000011921637810251895922488,
		2.00000002980310965673008246932, -2.00000000745066804957691410965,
		2.00000000186265486483933636576, -2.00000000046566236673530109841,
		2.00000000011641544175805401779, -2.00000000002910384378208396848,
		2.00000000000727595909475730239 } ;
			/* list of n! zeta(n+1) up to n=7*/
	static double nzetan[] = { 0.,1.64493406684822643647241516665,2.40411380631918857079947632302,
		6.49393940226682914909602217926, 24.8862661234408782319527716750,122.081167438133896765742151575,
		726.011479714984435324654235892, 5060.54987523763947046857360209} ;

		// constrained to the list of nzetan[] given above
	if (x<0. || n > 7 || n <1)
		return -1. ;
			/* for values up to 1.80 the list of zeta functions
			and the sum up to k < K are huge enough to gain
			numeric stability in the sum */
	if(x>1.9)
	{
		// Abramowitz Stegun 27.1.2
		for(k=1;;k++)
		{
			sumold=sum ;
				/* do not use x(k+1)=xk+x to avoid loss of precision */
			xk=x*k ;
			double ksum= 1./xk ;
			double tmp = n*ksum/xk ;	// n/(xk)^2
			for (int s=1 ; s <= n ; s++)
			{
				ksum += tmp ;
				tmp *= (n-s)/xk ;
			}
			sum += exp(-xk)* ksum*pow(x,n+1) ;
			if( sum == sumold)
				break ;
		}
				/* n!*zeta(n) is the integral for x=infinity , 27.1.3 */
		return nzetan[n]-sum ;
	}
	else
	{
		// Abramowitz-Stegun 27.1.1
		register double x2pi=x*M_1_2PI ;
		for(k=1;k<K;k++)
		{
			sumold=sum ;
				/* do not precompute x2pi^2 do avoid loss of precision */
			sum += koeff[k]*pow(x2pi,2*k)/(2*k+n) ;
			if(sum == sumold)
				break ;
		}
		sum += 1./double(n)-x/(2*(1+n)) ;
		return sum*pow(x,n) ;
	}
}

#undef M_1_2PI
#undef zeta_4
#undef K