#ifndef RTS_BEAM
#define RTS_BEAM

#include "../math/vec3.h"
#include "../optics/scalarwave.h"
#include "../math/bessel.h"
#include <vector>

//Boost
//#include <boost/math/special_functions/bessel.hpp>
//#include <boost/math/special_functions/legendre.hpp>

namespace stim{

		/// Function returns the value of the scalar field produced by a beam with the specified parameters

		template<typename T>
		std::vector< stim::vec3<T> > generate_focusing_vectors(size_t N, stim::vec3<T> d, T NA, T NA_in = 0){

			std::vector< stim::vec3<T> > dirs(N);					//allocate an array to store the focusing vectors

			///compute the rotation operator to transform (0, 0, 1) to k
			T cos_angle = d.dot(vec3<T>(0, 0, 1));
			stim::matrix<T, 3> rotation;

			//if the cosine of the angle is -1, the rotation is just a flip across the z axis
			if(cos_angle == -1){
				rotation(2, 2) = -1;
			}
			else if(cos_angle != 1.0)
			{
				vec3<T> r_axis = vec3<T>(0, 0, 1).cross(d).norm();	//compute the axis of rotation
				T angle = acos(cos_angle);							//compute the angle of rotation
				quaternion<T> quat;							//create a quaternion describing the rotation
				quat.CreateRotation(angle, r_axis);
				rotation = quat.toMatrix3();							//compute the rotation matrix
			}

			//find the phi values associated with the cassegrain ring
			T PHI[2];
			PHI[0] = (T)asin(NA);
			PHI[1] = (T)asin(NA_in);

			//calculate the z-axis cylinder coordinates associated with these angles
			T Z[2];
			Z[0] = cos(PHI[0]);
			Z[1] = cos(PHI[1]);
			T range = Z[0] - Z[1];

			//draw a distribution of random phi, z values
			T z, phi, theta;
			//T kmag = stim::TAU / lambda;
			for(int i=0; i<N; i++){								//for each sample
				z = (T)((double)rand() / (double)RAND_MAX) * range + Z[1];			//find a random position on the surface of a cylinder
				theta = (T)(((double)rand() / (double)RAND_MAX) * stim::TAU);
				phi = acos(z);													//project onto the sphere, computing phi in spherical coordinates

				//compute and store cartesian coordinates
				vec3<T> spherical(1, theta, phi);								//convert from spherical to cartesian coordinates
				vec3<T> cart = spherical.sph2cart();
				dirs[i] = rotation * cart;										//create a sample vector
			}
			return dirs;
		}
		
/// Class stim::beam represents a beam of light focused at a point and composed of several plane waves
template<typename T>
class scalarbeam
{
public:
	//enum beam_type {Uniform, Bartlett, Hamming, Hanning};

private:
	
	T NA[2];				//numerical aperature of the focusing optics	
	vec3<T> f;				//focal point
	vec3<T> d;				//propagation direction
	stim::complex<T> A;		//beam amplitude
	T lambda;				//beam wavelength
public:

	///constructor: build a default beam (NA=1.0)
	scalarbeam(T wavelength = 1, stim::complex<T> amplitude = 1, vec3<T> focal_point = vec3<T>(0, 0, 0), vec3<T> direction = vec3<T>(0, 0, 1), T numerical_aperture = 1, T center_obsc = 0){
		lambda = wavelength;
		A = amplitude;
		f = focal_point;
		d = direction.norm();					//make sure that the direction vector is normalized (makes calculations more efficient later on)
		NA[0] = numerical_aperture;
		NA[1] = center_obsc;
	}

	///Numerical Aperature functions
	void setNA(T na)
	{
		NA[0] = (T)0;
		NA[1] = na;
	}
	void setNA(T na0, T na1)
	{
		NA[0] = na0;
		NA[1] = na1;
	}

	//Monte-Carlo decomposition into plane waves
	std::vector< scalarwave<T> > mc(size_t N = 100000) const{

		std::vector< stim::vec3<T> > dirs = generate_focusing_vectors(N, d, NA[0], NA[1]);	//generate a random set of N vectors forming a focus
		std::vector< scalarwave<T> > samples(N);											//create a vector of plane waves
		T kmag = (T)stim::TAU / lambda;								//calculate the wavenumber
		stim::complex<T> apw;										//allocate space for the amplitude at the focal point
		stim::vec3<T> kpw;											//declare the new k-vector based on the focused plane wave direction
		for(size_t i=0; i<N; i++){										//for each sample
			kpw = dirs[i] * kmag;									//calculate the k-vector for the new plane wave
			apw = exp(stim::complex<T>(0, kpw.dot(-f)));				//calculate the amplitude for the new plane wave
			samples[i] = scalarwave<T>(kpw, apw);			//create a plane wave based on the direction
		}

		return samples;
	}

	/// Calculate the field at a given point
	/// @param x is the x-coordinate of the field point
	/// @O is the approximation accuracy
	stim::complex<T> field(T x, T y, T z, size_t O){
		std::vector< scalarwave<T> > W = mc(O);
		T result = 0;											//initialize the result to zero (0)
		for(size_t i = 0; i < O; i++){							//for each plane wave
			result += W[i].pos(x, y, z);
		}
		return result;
	}

	/// Calculate the field at a set of positions
	/*void field(stim::complex<T>* F, T* x, T* y, T* z, size_t N, size_t O){

		memset(F, 0, N * sizeof(stim::complex<T>));
		std::vector< scalarwave<T> > W = mc(O);					//get a random set of plane waves representing the beam
		size_t o, n;
		T px, py, pz;
		for(n = 0; n < N; n++){									//for each point in the field
			(x == NULL) ? px = 0 : px = x[n];								// test for NULL values
			(y == NULL) ? py = 0 : py = y[n];
			(z == NULL) ? pz = 0 : pz = z[n];
			for(o = 0; o < O; o++){								//for each plane wave
				F[n] += W[o].pos(px, py, pz);
			}
		}
	}*/

	std::string str()
	{
		std::stringstream ss;
		ss<<"Beam:"<<std::endl;
		//ss<<"	Central Plane Wave: "<<beam::E0<<" e^i ( "<<beam::k<<" . r )"<<std::endl;
		ss<<"	Beam Direction: "<<d<<std::endl;
		if(NA[0] == 0)
			ss<<"	NA: "<<NA[1];
		else
			ss<<"	NA: "<<NA[0]<<" -- "<<NA[1];

		return ss.str();
	}



};			//end beam

template<typename T>
void cpu_scalar_psf(stim::complex<T>* F, size_t N, T* x, T* y, T* z, T lambda, T A, stim::vec3<T> f, T NA, T NA_in, int Nl){

	memset(F, 0, N * sizeof(stim::complex<T>));
	T k = stim::TAU / lambda;
	T jl, Pl, C, kr, cos_phi;
	T cos_alpha_1 = cos(asin(NA_in));
	T cos_alpha_2 = cos(asin(NA));
	stim::vec3<T> p, ps;

	/*double vm;
	size_t table_bytes = (Nl + 1) * sizeof(double);
	double* jv = (double*) malloc( table_bytes );
	double* yv = (double*) malloc( table_bytes );
	double* djv= (double*) malloc( table_bytes );
	double* dyv= (double*) malloc( table_bytes );
	*/

	T vm;
	size_t table_bytes = (Nl + 1) * sizeof(T);
	T* jv = (T*) malloc( table_bytes );
	T* yv = (T*) malloc( table_bytes );
	T* djv= (T*) malloc( table_bytes );
	T* dyv= (T*) malloc( table_bytes );

	for(size_t n = 0; n < N; n++){
		(x == NULL) ? p[0] = 0 : p[0] = x[n];								// test for NULL values and set positions
		(y == NULL) ? p[1] = 0 : p[1] = y[n];
		(z == NULL) ? p[2] = 0 : p[2] = z[n];

		ps = p.cart2sph();												//convert this point to spherical coordinates
		kr = k * ps[0];
		cos_phi = std::cos(ps[2]);
		stim::bessjyv_sph<T>(Nl, kr, vm, jv, yv, djv, dyv);

		for(int l = 0; l <= Nl; l++){
			//jl = boost::math::sph_bessel<T>(l, kr);
			//jl = stim::bessjyv(l, kr
			jl = (T)jv[l];
			Pl = 1;//boost::math::legendre_p<T>(l, cos_phi);
			C = 1;//boost::math::legendre_p<T>(l+1, cos_alpha_1) - boost::math::legendre_p<T>(l + 1, cos_alpha_2) - boost::math::legendre_p<T>(l - 1, cos_alpha_1) + boost::math::legendre_p<T>(l - 1, cos_alpha_2);
			F[n] += pow(complex<T>(0, 1), l) * jl * Pl * C;
		}
		F[n] *= A * stim::TAU;
	}
}

}			//end namespace stim

#endif