#ifndef RTS_BEAM
#define RTS_BEAM

#include "../math/vector.h"
#include "../math/function.h"
#include "../optics/planewave.h"
#include <vector>

namespace stim{

template<typename P>
class beam : public planewave<P>
{
public:
	enum beam_type {Uniform, Bartlett, Hamming, Hanning};

private:
	
	P _na[2];		//numerical aperature of the focusing optics	
	vec<P> f;		//focal point	
	function<P, P> apod;	//apodization function
	unsigned int apod_res;	//resolution of apodization filter functions

	void apod_uniform()
	{
		apod = (P)1;
	}
	void apod_bartlett()
	{
		apod = (P)1;
		apod.insert((P)1, (P)0);
	}
	void apod_hanning()
	{
		apod = (P)0;
		P x, y;
		for(unsigned int n=0; n<apod_res; n++)
		{
			x = (P)n/(P)apod_res;
			y = pow( cos( ((P)3.14159 * x) / 2 ), 2);
			apod.insert(x, y);
		}
	}
	void apod_hamming()
	{
		apod = (P)0;
		P x, y;
		for(unsigned int n=0; n<apod_res; n++)
		{
			x = (P)n/(P)apod_res;
			y = (P)27/(P)50 + ( (P)23/(P)50 ) * cos((P)3.14159 * x);
			apod.insert(x, y);
		}
	}

	void set_apod(beam_type type)
	{
		if(type == Uniform)
			apod_uniform();
		if(type == Bartlett)
			apod_bartlett();
		if(type == Hanning)
			apod_hanning();
		if(type == Hamming)
			apod_hamming();
	}

public:

	///constructor: build a default beam (NA=1.0)
	beam(
		vec<P> k = rts::vec<P>(0, 0, rtsTAU), 
		vec<P> _E0 = rts::vec<P>(1, 0, 0), 
		beam_type _apod = Uniform)
		: planewave<P>(k, _E0)
	{
		_na[0] = (P)0.0;
		_na[1] = (P)1.0;
		f = vec<P>( (P)0, (P)0, (P)0 );
		apod_res = 256;						//set the default resolution for apodization filters
		set_apod(_apod);						//set the apodization function type
	}

	beam<P> refract(rts::vec<P> kn) const{

		beam<P> new_beam;
		new_beam._na[0] = _na[0];
		new_beam._na[1] = _na[1];


		rts::planewave<P> pw = planewave<P>::bend(kn);
		//std::cout<<pw.str()<<std::endl;

		new_beam.k = pw.kvec();
		new_beam.E0 = pw.E();

		return new_beam;
	}

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

	/*string str() : 
	{
		stringstream ss;
		ss<<"Beam Center: "<<k<<std::endl;

		return ss.str();
	}*/

	//Monte-Carlo decomposition into plane waves
	std::vector< planewave<P> > mc(unsigned int N = 100000, unsigned int seed = 0) const
	{
		/*Create Monte-Carlo samples of a cassegrain objective by performing uniform sampling
			of a sphere and projecting these samples onto an inscribed sphere.

			seed	=	seed for the random number generator
		*/
		srand(seed);		//seed the random number generator

		vec<P> k_hat = beam::k.norm();

		///compute the rotation operator to transform (0, 0, 1) to k
		P cos_angle = k_hat.dot(rts::vec<P>(0, 0, 1));
		rts::matrix<P, 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)
		{
			rts::vec<P> r_axis = rts::vec<P>(0, 0, 1).cross(k_hat).norm();	//compute the axis of rotation
			P angle = acos(cos_angle);							//compute the angle of rotation
			rts::quaternion<P> 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
		P PHI[2];
		PHI[0] = (P)asin(_na[0]);
		PHI[1] = (P)asin(_na[1]);

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

		std::vector< planewave<P> > samples;	//create a vector of plane waves

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

			//compute and store cartesian coordinates
			rts::vec<P> spherical(1, theta, phi);				//convert from spherical to cartesian coordinates
			rts::vec<P> cart = spherical.sph2cart();
			vec<P> k_prime = rotation * cart;				//create a sample vector

			//store a wave refracted along the given direction
			//std::cout<<"k prime: "<<rotation<<std::endl;
			samples.push_back(planewave<P>::refract(k_prime) * apod(phi/PHI[1]));
		}

		return samples;
	}

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

		return ss.str();
	}



};

}

#endif