/*
Orion Sky Lawlor, olawlor@acm.org, 9/23/2000

Computes so-called "Mie scattering" for a homogenous
sphere of arbitrary size illuminated by coherent
harmonic radiation.
*/
#include <stdio.h>
#include <stdlib.h>
#include <cmath>
#include <complex>

using namespace std;

/******************************************
Computes modified legendre polynomials
*/
class legendre { public:
	//Compute Pi and Tau of order n at z
	static void piTau(int n,double z,double *pi,double *tau);
};

/*Compute Pi and Tau of order n at z
using simple recurrence relation.
*/
void legendre::piTau(int n_stop,double z,double *pi,double *tau)
{
	double p[4],d[4];
	p[0]=0;   d[0]=0;
	p[1]=1;   d[1]=0;
	p[2]=3*z; d[2]=3;
	int n=2;
	if (n>=n_stop) n=n_stop;
	else while (n<n_stop)
	{
		n++;
		p[n&3]=z*(2*n-1)/(n-1)*p[(n-1)&3] - double(n)/(n-1)*p[(n-2)&3];
		d[n&3]=(2*n-1)*p[(n-1)&3] + d[(n-2)&3];
	}
	*pi=p[n&3];
	*tau=z*p[n&3]-pow(sin(acos(z)),2)*d[n&3];
}

/********************************************
Computes the Ricatti-Bessel function Psi_n(z).
*/
class ricatti_bessel { public:
	//Compute the R/B function of order n at z
	// Picks the most accurate function below
	static double psi(int n,double z) {
		if (z<1.0+0.3*n) return inf_series(n,z);
		else return recurse(n,z);
	}
	//Compute derivative with respect to z of
	// the R/B function of order n at z (via finite difference)
	static double psiPrime(int n,double z) {
		double dz=1e-11*z;
		return 0.5*(psi(n,z+dz)-psi(n,z-dz))/dz;
	}
	
	//Compute the R/B function of order n at z
	// using an infinite series.
	//  Blows up (does not converge) for large z (>20)
	//  Converges to 1e-15 in under 10 iterations for z<1.0
	static double inf_series(int n,double z,int *nIt=NULL);
	
	//Compute the R/B function of order n at z
	// using a simple recursion formula.
	//  Blows up for small z (>0.01) and reasonable n
	static double recurse(int n,double z);
private:
	static double fact(int n) {
		register double ret=1;
		while (n>1) {
			ret*=n;
			n--;
		}
		return ret;
	}
};

//Compute the R/B function of order n at z
// using an infinite series.
double ricatti_bessel::inf_series(int n,double z,int *nIt)
{
	long double sum=0;
	int m=0;
	double err=10.0;
	while ((fabs(err)>1e-15*fabs(sum)) && (m<100000)) {
		err=((m&1)?-1:1)*fact(n+m)/ (
			fact(m)*fact(2*n+2*m+1)) * pow(z,2*m);
		sum+=err;
		m++;
	}
	if (nIt!=NULL) *nIt=m;
	return z*pow(2*z,n)*sum;
}

//Compute the R/B function of order n at z
// using a simple recursion formula:
//  Zn+1(z)=(2(n+1/2)/z)*Zn(z)+Zn-1(z)
//  Z0(z)=sin(z), Z1(z)=sinz/z - cos(z)
double ricatti_bessel::recurse(int n_stop,double z)
{
	register long double Znm=sin(z);
	if (n_stop==0) return Znm;
	register long double Zn=Znm/z-cos(z);
	int n=1;
	register long double i_z=2.0/z;
	while (n<n_stop) {
		register long double Znp=(n+0.5)*i_z*Zn-Znm;
		n++;Znm=Zn;Zn=Znp;
	}
	return Zn;
}

/**********************************************
Scattering:
	Now that the preliminaries are done, we actually
calculate the scattering efficiency at various angles.
This solution is a spherical harmonic solution to 
Maxwell's equations for an isotropic homogenous sphere
of refractive index m in vacuum.
*/

class scatter_mie {public:
	typedef complex<double> cpx;
	static const cpx i;//Imaginary unit

	scatter_mie(double l,double r,cpx Nm,int Nmax_n)
		:lambda(l),radius(r),m(Nm),max_n(Nmax_n) 
	{
		k0=2*M_PI/lambda;k1=m*k0;
		An=new cpx[max_n];
		Bn=new cpx[max_n];
		for (int n=1;n<max_n;n++) {
			cpx a=k0*radius,b=k1*radius;
			cpx pa =ricatti_bessel::psi(n,real(a));
			cpx pb =ricatti_bessel::psi(n,real(b));
			cpx paP=ricatti_bessel::psiPrime(n,real(a));
			cpx pbP=ricatti_bessel::psiPrime(n,real(b));
			cpx za =ff_zeta (n,a);
			cpx zaP=ff_zetaP(n,a);
			An[n]=(  pa*pbP-m*pb*paP)/(  za*pbP-m*pb*zaP);
			Bn[n]=(m*pa*pbP-  pb*paP)/(m*za*pbP-  pb*zaP);
		}
	}
	~scatter_mie() {delete[]An;delete[]Bn;}
	
	/*Return the scattering coefficient at the
	angle theta.
	*/
	double ang(double theta) {
		cpx s1(0),s2(0);
		double cosTheta=cos(theta);
		for (int n=1;n<max_n;n++) {
			double pi_n,tau_n;
			legendre::piTau(n,cosTheta,&pi_n,&tau_n);
			s1+=(2*n+1)/double(n*(n+1))*(An[n]* pi_n+Bn[n]*tau_n);
			s2+=(2*n+1)/double(n*(n+1))*(An[n]*tau_n+Bn[n]* pi_n);
		}
		return lambda*lambda*real(s1*conj(s1)+s2*conj(s2))/(radius*radius);
	}
protected:
	//Parameters
	double lambda;//Wavelength of impinging radiation in free space
	double radius;//Sphere size
	cpx m;//Sphere index of refraction
	int max_n;//Largest spherical harmonic to consider

	//Computed values
	double k0; cpx k1;//Wavenumbers in free space and in sphere
	cpx *An,*Bn;
	
	//Computes the far-field zeta of order n at z
	cpx ff_zeta(int n,cpx z) const {
		static const cpx i2n[4]={1,i,-1,-i};
		return i2n[(n+1)&3]*exp(-i*z);
	}
	cpx ff_zetaP(int n,cpx z) const { //Derivative of above w.r.t. z
		static const cpx i2n[4]={1,i,-1,-i};
		return i2n[n&3]*exp(-i*z);
	}
	
};

const scatter_mie::cpx scatter_mie::i(0,1);

#if 1
int main() {
	const int ns=10;
	scatter_mie *s[ns];
	for (int i=0;i<ns;i++)
		s[i]=new scatter_mie(600e-9,(i*1+1)*600e-9,complex<double>(1.5,0),6);
	for (double ang=0;ang<=180;ang+=5) {
		printf("%.0f\t",ang);
		for (int i=0;i<ns;i++)
			printf("%e\t",s[i]->ang(ang*M_PI/180.0));
		printf("\n");
	}
	return 0;
}
#endif
/*<@>
<@> ******** Program output: ********
<@> 0	6.237791e+02	2.845479e+01	5.267256e+00	1.612578e+00	6.480777e-01	3.088940e-01	1.654605e-01	9.648105e-02	6.001013e-02	3.926608e-02	
<@> 5	5.683599e+02	2.609738e+01	4.832564e+00	1.479674e+00	5.947078e-01	2.834698e-01	1.518469e-01	8.854496e-02	5.507489e-02	3.603728e-02	
<@> 10	4.244810e+02	1.992188e+01	3.695010e+00	1.132153e+00	4.552230e-01	2.170424e-01	1.162849e-01	6.781676e-02	4.218580e-02	2.760539e-02	
<@> 15	2.479386e+02	1.218225e+01	2.272140e+00	6.981719e-01	2.812077e-01	1.342215e-01	7.196432e-02	4.199047e-02	2.612978e-02	1.710324e-02	
<@> 20	1.003282e+02	5.420702e+00	1.032226e+00	3.209081e-01	1.301617e-01	6.240011e-02	3.355363e-02	1.961731e-02	1.222470e-02	8.009990e-03	
<@> 25	2.050489e+01	1.343543e+00	2.855218e-01	9.436742e-02	3.963199e-02	1.940506e-02	1.057736e-02	6.241336e-03	3.914593e-03	2.577097e-03	
<@> 30	1.029407e+01	2.192822e-01	7.564882e-02	3.058932e-02	1.413993e-02	7.297575e-03	4.106921e-03	2.474279e-03	1.574198e-03	1.046995e-03	
<@> 35	4.069457e+01	1.051695e+00	2.163619e-01	7.187447e-02	3.032725e-02	1.489430e-02	8.134534e-03	4.806145e-03	3.017293e-03	1.987736e-03	
<@> 40	7.285804e+01	2.337141e+00	4.346277e-01	1.353138e-01	5.502750e-02	2.643242e-02	1.423276e-02	8.329264e-03	5.194177e-03	3.405181e-03	
<@> 45	8.042041e+01	2.920809e+00	5.270758e-01	1.604837e-01	6.439035e-02	3.067196e-02	1.642556e-02	9.576793e-03	5.956370e-03	3.897300e-03	
<@> 50	6.060088e+01	2.479867e+00	4.429470e-01	1.336937e-01	5.335704e-02	2.533272e-02	1.353727e-02	7.881357e-03	4.896809e-03	3.201595e-03	
<@> 55	2.963895e+01	1.443739e+00	2.627244e-01	8.011544e-02	3.218052e-02	1.534241e-02	8.221670e-03	4.796014e-03	2.983989e-03	1.952976e-03	
<@> 60	7.985425e+00	5.145868e-01	1.078091e-01	3.563750e-02	1.499194e-02	7.350919e-03	4.010993e-03	2.368545e-03	1.486385e-03	9.789329e-04	
<@> 65	6.005198e+00	1.536088e-01	5.087688e-02	2.024697e-02	9.292744e-03	4.777527e-03	2.682325e-03	1.613526e-03	1.025499e-03	6.815475e-04	
<@> 70	1.898053e+01	3.482519e-01	8.193944e-02	2.928624e-02	1.282767e-02	6.434134e-03	3.559787e-03	2.121180e-03	1.339520e-03	8.861875e-04	
<@> 75	3.299693e+01	7.398634e-01	1.382018e-01	4.386433e-02	1.807219e-02	8.753490e-03	4.739222e-03	2.783856e-03	1.740654e-03	1.143364e-03	
<@> 80	3.605163e+01	9.475249e-01	1.621485e-01	4.850930e-02	1.930282e-02	9.152431e-03	4.887403e-03	2.844203e-03	1.766673e-03	1.154866e-03	
<@> 85	2.612357e+01	8.299975e-01	1.406378e-01	4.157388e-02	1.641645e-02	7.746043e-03	4.123199e-03	2.394265e-03	1.484861e-03	9.695269e-04	
<@> 90	1.113561e+01	5.217434e-01	1.011593e-01	3.218665e-02	1.325754e-02	6.418883e-03	3.474270e-03	2.040477e-03	1.275617e-03	8.377877e-04	
<@> 95	1.925858e+00	2.674491e-01	7.788664e-02	2.909520e-02	1.294945e-02	6.545117e-03	3.636819e-03	2.172933e-03	1.374604e-03	9.105074e-04	
<@> 100	4.023885e+00	2.172706e-01	8.110428e-02	3.242866e-02	1.484037e-02	7.607776e-03	4.262078e-03	2.559719e-03	1.624988e-03	1.079035e-03	
<@> 105	1.417979e+01	3.422439e-01	9.515045e-02	3.541833e-02	1.574270e-02	7.950080e-03	4.414705e-03	2.636428e-03	1.667313e-03	1.104144e-03	
<@> 110	2.359768e+01	5.086494e-01	1.021100e-01	3.350720e-02	1.403172e-02	6.857620e-03	3.732904e-03	2.200421e-03	1.379177e-03	9.074854e-04	
<@> 115	2.487819e+01	6.182619e-01	1.054754e-01	3.137851e-02	1.242687e-02	5.871632e-03	3.127589e-03	1.816780e-03	1.126990e-03	7.359670e-04	
<@> 120	1.733253e+01	6.819196e-01	1.289234e-01	4.011518e-02	1.627905e-02	7.805598e-03	4.197415e-03	2.454069e-03	1.529235e-03	1.001967e-03	
<@> 125	6.996467e+00	7.691441e-01	1.890598e-01	6.570498e-02	2.821508e-02	1.397673e-02	7.670636e-03	4.545910e-03	2.859614e-03	1.886489e-03	
<@> 130	1.799601e+00	8.962894e-01	2.672764e-01	9.896913e-02	4.374452e-02	2.201042e-02	1.219365e-02	7.270530e-03	4.592723e-03	3.038910e-03	
<@> 135	5.774035e+00	9.765433e-01	3.116865e-01	1.179101e-01	5.260535e-02	2.659898e-02	1.477830e-02	8.827861e-03	5.583544e-03	3.697849e-03	
<@> 140	1.652464e+01	9.039396e-01	2.782108e-01	1.046273e-01	4.656953e-02	2.351707e-02	1.305551e-02	7.794481e-03	4.928190e-03	3.262989e-03	
<@> 145	2.745403e+01	7.199225e-01	1.855795e-01	6.640628e-02	2.891133e-02	1.442571e-02	7.949976e-03	4.723516e-03	2.976777e-03	1.966328e-03	
<@> 150	3.267587e+01	7.166872e-01	1.394468e-01	4.413924e-02	1.804041e-02	8.677324e-03	4.672592e-03	2.733485e-03	1.704176e-03	1.116930e-03	
<@> 155	3.081923e+01	1.350688e+00	2.943462e-01	9.705094e-02	4.052575e-02	1.974587e-02	1.072299e-02	6.309648e-03	3.949556e-03	2.596194e-03	
<@> 160	2.513555e+01	2.970834e+00	7.633413e-01	2.682116e-01	1.156500e-01	5.740066e-02	3.153569e-02	1.870100e-02	1.176905e-02	7.766411e-03	
<@> 165	2.032285e+01	5.517690e+00	1.526552e+00	5.503188e-01	2.402099e-01	1.200374e-01	6.622446e-02	3.938007e-02	2.482991e-02	1.640762e-02	
<@> 170	1.888372e+01	8.406443e+00	2.400672e+00	8.747265e-01	3.837203e-01	1.922790e-01	1.062575e-01	6.325473e-02	3.991334e-02	2.638899e-02	
<@> 175	1.986450e+01	1.071336e+01	3.100579e+00	1.134847e+00	4.988710e-01	2.502666e-01	1.383991e-01	8.242593e-02	5.202648e-02	3.440538e-02	
<@> 180	2.063923e+01	1.159556e+01	3.368382e+00	1.234421e+00	5.429608e-01	2.724721e-01	1.507082e-01	8.976820e-02	5.666579e-02	3.747572e-02	
<@> */