#ifndef STIM_SPH_HARMONICS
#define STIM_SPH_HARMONICS

#include <stim/math/vector.h>
#include <boost/math/special_functions/spherical_harmonic.hpp>
#include <vector>

#define PI 3.14159
#define WIRE_SCALE 1.001
namespace stim{

template<class T>
class spharmonics{

protected:

	std::vector<T> C;	//list of SH coefficients

	unsigned int mcN;	//number of Monte-Carlo samples

	//calculate the value of the SH basis function (l, m) at (theta, phi)
		//here, theta = [0, PI], phi = [0, 2*PI]
	double SH(int l, int m, double theta, double phi){
		return boost::math::spherical_harmonic_r(l, m, phi, theta);
	}

	unsigned int coeff_1d(unsigned int l, int m){
		return pow(l + 1, 2) - (l - m) - 1;
	}

	


public:

	void push(double c){
		C.push_back(c);
	}

	void resize(unsigned int n){
		C.resize(n);
	}

	void setc(unsigned int l, int m, T value){
		unsigned int c = coeff_1d(l, m);
		C[c] = value;
	}

	void setc(unsigned int c, T value){
		C[c] = value;
	}

	/// Initialize Monte-Carlo sampling of a function using N spherical harmonics coefficients

	/// @param N is the number of spherical harmonics coefficients used to represent the user function
	void mcBegin(unsigned int coefficients){
		C.resize(coefficients, 0);
		mcN = 0;
	}

	void mcBegin(unsigned int l, int m){
		unsigned int c = pow(l + 1, 2) - (l - m);
		mcBegin(c);
	}

	void mcSample(double theta, double phi, double val){

		int l, m;
		double sh;

		l = m = 0;
		for(unsigned int i = 0; i < C.size(); i++){

			sh = SH(l, m, theta, phi);
			C[i] += sh * val;

			m++;			//increment m

			//if we're in a new tier, increment l and set m = -l
			if(m > l){		
				l++;
				m = -l;
			}
		}	//end for all coefficients

		//increment the number of samples
		mcN++;

	}	//end mcSample()

	void mcEnd(){

		//divide all coefficients by the number of samples
		for(unsigned int i = 0; i < C.size(); i++)
			C[i] /= mcN;
	}

	/// Generates a PDF describing the probability distribution of points on a spherical surface

	/// @param sph_pts is a list of points in spherical coordinates (theta, phi) where theta = [0, 2pi] and phi = [0, pi]
	/// @param l is the maximum degree of the spherical harmonic function
	/// @param m is the maximum order
	void pdf(std::vector<stim::vec<double> > sph_pts, unsigned int l, int m){
			
		mcBegin( l, m );		//begin spherical harmonic sampling

		unsigned int nP = sph_pts.size();

		for(unsigned int p = 0; p < nP; p++){
			mcSample(sph_pts[p][1], sph_pts[p][2], 1.0);
		}

		mcEnd();
	}

	std::string str(){

		std::stringstream ss;

		int l, m;
		l = m = 0;
		for(unsigned int i = 0; i < C.size(); i++){
				
			ss<<C[i]<<'\t';

			m++;			//increment m

			//if we're in a new tier, increment l and set m = -l
			if(m > l){
				l++;
				m = -l;

				ss<<std::endl;
					
			}
		}

		return ss.str();


	}

	/// Returns the value of the function at the coordinate (theta, phi)

	/// @param theta = [0, 2pi]
	/// @param phi = [0, pi]
	double operator()(double theta, double phi){

		double fx = 0;

		int l = 0;
		int m = 0;
		for(unsigned int i = 0; i < C.size(); i++){
			fx += C[i] * SH(l, m, theta, phi);
			m++;
			if(m > l){
				l++;
				m = -l;					
			}

		}

		return fx;
	}

};		//end class sph_harmonics




}


#endif