#ifndef RTS_BEAM
#define RTS_BEAM

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

namespace rts{

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

private:
	
	P na[2];		//numerical aperature of the focusing optics	
	vec<P> f;		//focal point	
	vec<P> k;		//direction vector	
	vec<P> E0;		//polarization direction
	P omega;		//frequency

	function<P, P> apod;	//apodization function
	unsigned int apod_res;	//resolution of complex apodization filters

	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(beam_type _apod = Uniform)
	{
		na[0] = (P)0.0;
		na[1] = (P)1.0;
		f = vec<P>( (P)0.0, (P)0.0, (P)0.0 );
		k = vec<P>( (P)0.0, (P)0.0, (P)1.0 );
		E0 = vec<P>( (P)1.0, (P)0.0, (P)0.0 );
		omega = (P)2 * (P)3.14159;
		apod_res = 256;						//set the default resolution for apodization filters
		set_apod(_apod);						//set the apodization function type
		
	}

	///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;
	}


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

			samples = rtsPointer to sample vectors specified as normalized cartesian coordinates
			N       = number of samples
			kSph	= incident light direction in spherical coordinates
			NAin    = internal obscuration NA
			NAout   = outer cassegrain NA
		*/

		srand(seed);		//seed the random number generator

		///compute the rotation operator to transform (0, 0, 1) to k
		P cos_angle = k.dot(rts::vec<P>(0, 0, 1));
		rts::matrix<P, 3> rotation;
		if(cos_angle != 1.0)
		{
			rts::vec<P> r_axis = rts::vec<P>(0, 0, 1).cross(k).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

		planewave<P> beam_center(k, E0, omega);	//create a plane wave representing the beam center

		//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
			samples.push_back(beam_center.refract(k_prime) * apod(phi/PHI[1]));
		}

		return samples;
	}

};

}

#endif