//
//	a_n , b_n calc for Mie computation
//	c++ an_bn.cpp
//

#include <iostream>
#include <cmath>
#include <complex>

using std::cout;
using std::cerr;
using std::endl;
using std::complex;


int main(int argc, char* argv[])
{
	//	definition of the imaginary number, $i$
	const complex<double> Z_ONE = complex<double>(0.0, 1.0);


	complex<double> z = complex<double>(10.0, 0.0);


	// Riccati-Bessel functions and relating val for recurrence calcuration
	complex<double> psi_n;	//	== S
	complex<double> psi_n1;	//	psi_{n-1}
	complex<double> psi_n2;	//	psi_{n-2}

	complex<double> chi_n;	//	== C
	complex<double> chi_n1;	//	chi_{n-1}
	complex<double> chi_n2;	//	chi_{n-2}

	complex<double> xi_n;
	complex<double> xi_n1;//   = psi_n1 - Z_ONE * chi_n1;
	complex<double> zeta_n;
	complex<double> zeta_n1;// = psi_n1 + Z_ONE * chi_n1;

	//	derivatives
	complex<double> dpsi_n;
	complex<double> dzeta_n;



	//////////////////////////////////////////////
	//
	//	set initial values

	//	at n = -1
	psi_n2 = cos(z);	//	psi_{n-2}
	psi_n1 = sin(z);	//	psi_{n-1}

	chi_n2 = -sin(z);	//	chi_{n-2}
	chi_n1 =  cos(z);	//	chi_{n-1}

	//	at n = 0
	xi_n1   = psi_n1 - Z_ONE * chi_n1;
	zeta_n1 = psi_n1 + Z_ONE * chi_n1;

	const int MAX_n = 10;

	for(int n=1; n<=MAX_n; n++) {

		//	calc n'th values
		psi_n = (2.0*n-1.0)/(z) * psi_n1 - psi_n2;	//	psi_1, psi_2, ...
		chi_n = (2.0*n-1.0)/(z) * chi_n1 - chi_n2;	//	chi_1, chi_2, ...

		xi_n   = psi_n - Z_ONE * chi_n;
		zeta_n = psi_n + Z_ONE * chi_n;

		dpsi_n  = psi_n1  - (double)n / z * psi_n;
		dzeta_n = zeta_n1 - (double)n / z * zeta_n;

		cout << n << "," << z << ",," << psi_n << "," << chi_n << "," << xi_n << "," << zeta_n << "," << dpsi_n << "," << dzeta_n << "\n";

		{
			//	update values for the next loop
			psi_n2 = psi_n1;
			psi_n1  = psi_n;

			chi_n2 = chi_n1;
			chi_n1  = chi_n;
		}
	}

}