/*
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:
	typedef complex<double> cpx;

	//Compute the R/B function of order n at z
	// Picks the most accurate function below
	static cpx psi(int n,cpx z) {
		if (abs(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 cpx psiPrime(int n,cpx z) {
		cpx 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 cpx inf_series(int n,cpx 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 cpx recurse(int n,cpx 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.
ricatti_bessel::cpx ricatti_bessel::inf_series(int n,cpx z,int *nIt)
{
	cpx sum=0;
	int m=0;
	cpx err=10.0;
	while ((abs(err)>1e-15*abs(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;
	cpx p=pow(2.0*z,n);
	return z*p*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)
ricatti_bessel::cpx ricatti_bessel::recurse(int n_stop,cpx z)
{
	register cpx Znm=sin(z);
	if (n_stop==0) return Znm;
	register cpx Zn=Znm/z-cos(z);
	int n=1;
	register cpx i_z=2.0/z;
	while (n<n_stop) {
		register cpx 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,a);
			cpx pb =ricatti_bessel::psi(n,b);
			cpx paP=ricatti_bessel::psiPrime(n,a);
			cpx pbP=ricatti_bessel::psiPrime(n,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,(0.5)*600e-9,complex<double>(15,-0.0),i);
	for (double ang=0;ang<=180;ang+=5) {
		printf("%.0f\t",ang);
		for (int i=0;i<ns;i++)
			printf("%.15e\t",s[i]->ang(ang*M_PI/180.0));
		printf("\n");
	}
	return 0;
}
#endif
/*<@>
<@> ******** Program output: ********
<@> 0	0.000000000000000e+00	0.000000000000000e+00	1.807564153342996e+00	1.088195997838730e+02	3.843341069004833e+02	3.638074673349945e+02	2.819261120832588e+02	2.783183638892113e+02	2.856167915910687e+02	2.860214653875759e+02	
<@> 5	0.000000000000000e+00	0.000000000000000e+00	1.800698913209609e+00	1.069697857139339e+02	3.732914043955274e+02	3.534484298688645e+02	2.769496361173179e+02	2.737180078787511e+02	2.802146325998521e+02	2.805593985223245e+02	
<@> 10	0.000000000000000e+00	0.000000000000000e+00	1.780311791243833e+00	1.015691161178109e+02	3.416862886314333e+02	3.238601348582383e+02	2.621552685817657e+02	2.598862969793978e+02	2.642717360815136e+02	2.644684225010599e+02	
<@> 15	0.000000000000000e+00	0.000000000000000e+00	1.747022244578880e+00	9.304911349060120e+01	2.938106009141519e+02	2.791983204640229e+02	2.380364030519437e+02	2.369054683714753e+02	2.386310150348476e+02	2.386665560771277e+02	
<@> 20	0.000000000000000e+00	0.000000000000000e+00	1.701841767119110e+00	8.208128124515740e+01	2.359519657156754e+02	2.254399999795244e+02	2.056847455751026e+02	2.054232355475743e+02	2.048480536110453e+02	2.047812876891419e+02	
<@> 25	0.000000000000000e+00	0.000000000000000e+00	1.646143155749663e+00	6.951102737589839e+01	1.752863929338550e+02	1.692321415679229e+02	1.670954966982573e+02	1.671827988722522e+02	1.653354049019794e+02	1.652512753299386e+02	
<@> 30	0.000000000000000e+00	0.000000000000000e+00	1.581618798514515e+00	5.627546289797461e+01	1.186644530974336e+02	1.167282801661082e+02	1.253751785550733e+02	1.253419109615514e+02	1.234062652363145e+02	1.233606940038476e+02	
<@> 35	0.000000000000000e+00	0.000000000000000e+00	1.510229252156188e+00	4.331395398006938e+01	7.153365637485606e+01	7.267585647722072e+01	8.464555391966978e+01	8.430953739654953e+01	8.309138378947237e+01	8.308593913449810e+01	
<@> 40	0.000000000000000e+00	0.000000000000000e+00	1.434143671463375e+00	3.148099424501222e+01	3.720059314380941e+01	3.991113887196894e+01	4.949191382156695e+01	4.899107654824001e+01	4.864202302837968e+01	4.864304179079538e+01	
<@> 45	0.000000000000000e+00	0.000000000000000e+00	1.355673900451560e+00	2.147056975732183e+01	1.655404616438950e+01	1.928387605514659e+01	2.398392909154106e+01	2.360314676787323e+01	2.373711023741604e+01	2.371727549473484e+01	
<@> 50	0.000000000000000e+00	0.000000000000000e+00	1.277204227983901e+00	1.375976729717639e+01	8.265140782377525e+00	9.921911409317023e+00	1.054169425595174e+01	1.046921040239883e+01	1.056523681234333e+01	1.052873171887197e+01	
<@> 55	0.000000000000000e+00	0.000000000000000e+00	1.201118942173578e+00	8.577099207851001e+00	9.377887074031991e+00	9.691700186110296e+00	9.094712306705821e+00	9.288175833772611e+00	9.092434214707049e+00	9.064118171767575e+00	
<@> 60	0.000000000000000e+00	0.000000000000000e+00	1.129729884791841e+00	5.898156148472723e+00	1.612688968836378e+01	1.571929959407649e+01	1.695678235512467e+01	1.721004880177803e+01	1.686697334068863e+01	1.686283551956433e+01	
<@> 65	0.000000000000000e+00	0.000000000000000e+00	1.065206206905974e+00	5.468089053438511e+00	2.478499530418948e+01	2.488532096987969e+01	2.958067080598872e+01	2.972296161758415e+01	2.957981121713283e+01	2.959242678953597e+01	
<@> 70	0.000000000000000e+00	0.000000000000000e+00	1.009508460088243e+00	6.847425885894230e+00	3.235761243660335e+01	3.426044550384527e+01	4.196817616537928e+01	4.202171131115328e+01	4.230014699604039e+01	4.231227142276565e+01	
<@> 75	0.000000000000000e+00	0.000000000000000e+00	9.643290258028624e-01	9.475192749208949e+00	3.699583437790795e+01	4.147547537600309e+01	5.019330116532557e+01	5.035255242451255e+01	5.096510274652981e+01	5.097138267814638e+01	
<@> 80	0.000000000000000e+00	0.000000000000000e+00	9.310406929957851e-01	1.274152035974006e+01	3.808097907777230e+01	4.501268255067326e+01	5.242927731789243e+01	5.287846048700741e+01	5.352147854005784e+01	5.353122148553753e+01	
<@> 85	0.000000000000000e+00	0.000000000000000e+00	9.106549463331804e-01	1.606107245451036e+01	3.601711806615567e+01	4.437864395919854e+01	4.910324963739966e+01	4.983186624642897e+01	5.024275645193602e+01	5.026287984153944e+01	
<@> 90	0.000000000000000e+00	0.000000000000000e+00	9.037912324813715e-01	1.893881551052514e+01	3.183443116324055e+01	4.010190426420446e+01	4.221594212516653e+01	4.300066303834684e+01	4.314122510638791e+01	4.316426572578902e+01	
<@> 95	0.000000000000000e+00	0.000000000000000e+00	9.106581382588205e-01	2.102081992691132e+01	2.673944603784426e+01	3.351682387358351e+01	3.423551703830671e+01	3.481671698666756e+01	3.483842862700285e+01	3.485288121359905e+01	
<@> 100	0.000000000000000e+00	0.000000000000000e+00	9.310470525545833e-01	2.212478578944943e+01	2.174253670894116e+01	2.635297945237565e+01	2.711124676568842e+01	2.741927053265729e+01	2.747799825154485e+01	2.748702341070731e+01	
<@> 105	0.000000000000000e+00	0.000000000000000e+00	9.643385045947369e-01	2.224756891833954e+01	1.745313207795659e+01	2.022697120418109e+01	2.180102958885043e+01	2.203128149462122e+01	2.214488500009332e+01	2.216571500543976e+01	
<@> 110	0.000000000000000e+00	0.000000000000000e+00	1.009520985973864e+00	2.154982684844492e+01	1.407237485544184e+01	1.619260152115658e+01	1.838035400542056e+01	1.883386169696641e+01	1.890791234634000e+01	1.894804981546170e+01	
<@> 115	0.000000000000000e+00	0.000000000000000e+00	1.065221684555788e+00	2.032066456249768e+01	1.154969989192846e+01	1.450965092363374e+01	1.649995618479620e+01	1.730577096155735e+01	1.727331968173203e+01	1.731220963411784e+01	
<@> 120	0.000000000000000e+00	0.000000000000000e+00	1.129748196411588e+00	1.892750584295348e+01	9.821850749554578e+00	1.472587839319016e+01	1.584662629809575e+01	1.679713791579722e+01	1.671037892748618e+01	1.671378279429686e+01	
<@> 125	0.000000000000000e+00	0.000000000000000e+00	1.201139948400775e+00	1.775807206484888e+01	9.032922343297763e+00	1.604902051292192e+01	1.636523746325728e+01	1.702172963469574e+01	1.699546019029061e+01	1.695187983380223e+01	
<@> 130	0.000000000000000e+00	0.000000000000000e+00	1.277227768948475e+00	1.716212727639402e+01	9.646354789077368e+00	1.786146305735379e+01	1.821436078403325e+01	1.822766491427936e+01	1.829339922117529e+01	1.823163803852908e+01	
<@> 135	0.000000000000000e+00	0.000000000000000e+00	1.355699796992555e+00	1.740046674688204e+01	1.239935431342819e+01	2.015548604866562e+01	2.159656017558110e+01	2.100294995997648e+01	2.102718800824341e+01	2.099008050489664e+01	
<@> 140	0.000000000000000e+00	0.000000000000000e+00	1.434171726492478e+00	1.860752373523740e+01	1.810434392840648e+01	2.368123681397658e+01	2.662899115139160e+01	2.587490485841706e+01	2.563766216905896e+01	2.563502940384364e+01	
<@> 145	0.000000000000000e+00	0.000000000000000e+00	1.510259252157688e+00	2.077209595255642e+01	2.735429438791441e+01	2.970635961583157e+01	3.331721850290475e+01	3.295234706502233e+01	3.233468654930602e+01	3.232613277522226e+01	
<@> 150	0.000000000000000e+00	0.000000000000000e+00	1.581650515170286e+00	2.373829166140872e+01	4.022361074281338e+01	3.944765712307905e+01	4.157049402329351e+01	4.179406622760674e+01	4.092536716426685e+01	4.085836278039788e+01	
<@> 155	0.000000000000000e+00	0.000000000000000e+00	1.646176347676803e+00	2.722617044987101e+01	5.606339193429269e+01	5.338565071143518e+01	5.115694844453798e+01	5.152812930827475e+01	5.077635086748534e+01	5.065059157599737e+01	
<@> 160	0.000000000000000e+00	0.000000000000000e+00	1.701876181707012e+00	3.086903476522040e+01	7.346718212940404e+01	7.074627862127772e+01	6.157029434760967e+01	6.110832595769178e+01	6.091823787253954e+01	6.081722103384421e+01	
<@> 165	0.000000000000000e+00	0.000000000000000e+00	1.747057619911750e+00	3.426222847746336e+01	9.043724355320866e+01	8.939282032479065e+01	7.189848547819487e+01	6.955738893084300e+01	7.024350345547580e+01	7.029773265658444e+01	
<@> 170	0.000000000000000e+00	0.000000000000000e+00	1.780347858094028e+00	3.701687671409640e+01	1.047250506031431e+02	1.062250263059120e+02	8.084309433107114e+01	7.610707720395605e+01	7.771733658646551e+01	7.801928433255971e+01	
<@> 175	0.000000000000000e+00	0.000000000000000e+00	1.800735397086620e+00	3.881142953724693e+01	1.142692140581923e+02	1.179863002790172e+02	8.697521998373135e+01	8.023409572045540e+01	8.253377752819202e+01	8.306496543929244e+01	
<@> 180	0.000000000000000e+00	0.000000000000000e+00	1.807600776582490e+00	3.943421637309149e+01	1.176238450187905e+02	1.222123925293087e+02	8.916325830610711e+01	8.164164028377591e+01	8.419484802902407e+01	8.481908903552670e+01	
<@> */