#ifndef RTS_ESPHERE
#define RTS_ESPHERE

#include "../math/complex.h"
#include "../math/bessel.h"
#include "../visualization/colormap.h"
#include "../optics/planewave.h"
#include "../cuda/devices.h"
#include "../optics/efield.cuh"

namespace stim{

/*  This class implements a discrete representation of an electromagnetic field
    in 2D scattered by a sphere. This class implements Mie scattering.
*/
template<typename P>
class esphere : public efield<P>
{
private:
	stim::complex<P> n;		//sphere refractive index
	P a;					//sphere radius

	//parameters dependent on wavelength
	unsigned int Nl;		//number of orders for the calculation
	rts::complex<P>* A;		//internal scattering coefficients
	rts::complex<P>* B;		//external scattering coefficients

	void calcNl(P kmag)
	{
        //return ceil( ((P)6.282 * a) / lambda + 4 * pow( ((P)6.282 * a) / lambda, ((P)1/(P)3)) + 2);
		Nl = ceil( kmag*a + 4 * pow(kmag * a, (P)1/(P)3) + 2);
	}

	void calcAB(P k, unsigned int Nl, rts::complex<P>* A, rts::complex<P>* B)
	{
		/*  These calculations require double precision, so they are computed
			using doubles and converted to P at the end.

			Input:

			k	=	magnitude of the k vector (tau/lambda)
			Nl	=	highest order coefficient ([0 Nl] are computed)
		*/

		//clear the previous coefficients
		rts::complex<P>* cpuA = (rts::complex<P>*)malloc(sizeof(rts::complex<P>) * (Nl+1));
		rts::complex<P>* cpuB = (rts::complex<P>*)malloc(sizeof(rts::complex<P>) * (Nl+1));

		//convert to an std complex value
		complex<double> nc = (rts::complex<double>)n;

		//compute the magnitude of the k vector
		complex<double> kna = nc * k * (double)a;

		//compute the arguments k*a and k*n*a
		complex<double> ka(k * a, 0.0);

		//allocate space for the Bessel functions of the first and second kind (and derivatives)
		unsigned int bytes = sizeof(complex<double>) * (Nl + 1);
		complex<double>* cjv_ka = (complex<double>*)malloc(bytes);
		complex<double>* cyv_ka = (complex<double>*)malloc(bytes);
		complex<double>* cjvp_ka = (complex<double>*)malloc(bytes);
		complex<double>* cyvp_ka = (complex<double>*)malloc(bytes);
		complex<double>* cjv_kna = (complex<double>*)malloc(bytes);
		complex<double>* cyv_kna = (complex<double>*)malloc(bytes);
		complex<double>* cjvp_kna = (complex<double>*)malloc(bytes);
		complex<double>* cyvp_kna = (complex<double>*)malloc(bytes);

		//allocate space for the spherical Hankel functions and derivative
		complex<double>* chv_ka = (complex<double>*)malloc(bytes);
		complex<double>* chvp_ka = (complex<double>*)malloc(bytes);

		//compute the bessel functions using the CPU-based algorithm
		double vm;
		cbessjyva_sph(Nl, ka, vm, cjv_ka, cyv_ka, cjvp_ka, cyvp_ka);
		cbessjyva_sph(Nl, kna, vm, cjv_kna, cyv_kna, cjvp_kna, cyvp_kna);

		//compute A for each order
		complex<double> i(0, 1);
		complex<double> a, b, c, d;
		complex<double> An, Bn;
		for(int l=0; l<=Nl; l++)
		{
			//compute the Hankel functions from j and y
			chv_ka[l] = cjv_ka[l] + i * cyv_ka[l];
			chvp_ka[l] = cjvp_ka[l] + i * cyvp_ka[l];

			//Compute A (internal scattering coefficient)
			//compute the numerator and denominator for A
			a = cjv_ka[l] * chvp_ka[l] - cjvp_ka[l] * chv_ka[l];
			b = cjv_kna[l] * chvp_ka[l] - chv_ka[l] * cjvp_kna[l] * nc;

			//calculate A and add it to the list
			rts::complex<double> An = (2.0 * l + 1.0) * pow(i, l) * (a / b);
			cpuA[l] = (rts::complex<P>)An;

			//Compute B (external scattering coefficient)
			c = cjv_ka[l] * cjvp_kna[l] * nc - cjv_kna[l] * cjvp_ka[l];
			d = cjv_kna[l] * chvp_ka[l] - chv_ka[l] * cjvp_kna[l] * nc;

			//calculate B and add it to the list
			rts::complex<double> Bn = (2.0 * l + 1.0) * pow(i, l) * (c / d);
			cpuB[l] = (rts::complex<P>)Bn;

			std::cout<<"A: "<<cpuA[l]<<"     B: "<<cpuB[l]<<std::endl;
		}

		
		if(A != NULL)	cudaFree(A);	//free any previous coefficients
		if(B != NULL)	cudaFree(B);
		cudaMalloc(&A, sizeof(rts::complex<P>) * (Nl+1));	//allocate memory for new coefficients
		cudaMalloc(&B, sizeof(rts::complex<P>) * (Nl+1));

		//copy the calculations from the CPU to the GPU
		cudaMemcpy(A, cpuA, sizeof(rts::complex<P>) * (Nl+1), cudaMemcpyDeviceToHost);
		cudaMemcpy(B, cpuB, sizeof(rts::complex<P>) * (Nl+1), cudaMemcpyDeviceToHost);
	}

public:

	esphere(unsigned int res0, unsigned int res1, P _a, rts::complex<P> _n, bool _scalar = false) :
		efield<P>(res0, res1, _scalar)
	{
		std::cout<<"Sphere scattered field created."<<std::endl;
		n = _n;		//save the refractive index
		a = _a;		//save the radius
	}

	//assignment operator: build an electric field from a plane wave
    efield<P> & operator= (const planewave<P> & rhs)
	{
		calcNl(rhs.kmag());		//compute the number of scattering coefficients
		std::cout<<"Nl: "<<Nl<<std::endl;
		calcAB(rhs.kmag(), Nl, A, B);	//compute the scattering coefficients

		//determine important parameters for the scattering domain
		unsigned int sR = ceil(sqrt( (P)(pow((P)esphere::R[0],2) + pow((P)esphere::R[1],2))) );
		//unsigned int thetaR = 256;

		/////////////////////continue scattering code here/////////////////////////

		esphere::clear();				//clear any previous field data
		from_planewave(rhs);	//create a field from the planewave
		return *this;			//return the current object
	}

	string str()
	{
		stringstream ss;
		ss<<"Mie Scattered Field"<<std::endl;
		ss<<(*this).efield<P>::str()<<std::endl;
		ss<<"a = "<<a<<std::endl;
		ss<<"n = "<<n;

		return ss.str();
	}

};

}	//end namespace rts

#endif