#ifndef STIM_SPH_HARMONICS
#define STIM_SPH_HARMONICS

#include <GL/gl.h>

#include <stim/gl/error.h>
#include <stim/visualization/colormap.h>
#include <vector>

#define PI 3.14159
#define WIRE_SCALE 1.001
namespace stim{

	class sph_harmonics{

	private:

		double* func;		//stores the raw function data (samples at each point)

		GLuint color_tex;	//texture map that acts as a colormap for the spherical function
		
		unsigned int N;		//resolution of the spherical grid

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


		//evaluates an associated Legendre polynomial (-l <= m <= l)
		double P(int l,int m,double x)
		{
			// evaluate an Associated Legendre Polynomial P(l,m,x) at x
			double pmm = 1.0;
			if(m>0) {
				double somx2 = sqrt((1.0-x)*(1.0+x));
				double fact = 1.0;
				for(int i=1; i<=m; i++) {
					pmm *= (-fact) * somx2;
					fact += 2.0;
				}
			}
			if(l==m) return pmm;
			double pmmp1 = x * (2.0*m+1.0) * pmm;
			if(l==m+1) return pmmp1;
			double pll = 0.0;
			for(int ll=m+2; ll<=l; ++ll) {
				pll = ( (2.0*ll-1.0)*x*pmmp1-(ll+m-1.0)*pmm ) / (ll-m);
				pmm = pmmp1;
				pmmp1 = pll;
			}
			return pll;
		}

		//recursively calculate a factorial given a positive integer n
		unsigned int factorial(unsigned int n) {
			if (n == 0)
			   return 1;
			return n * factorial(n - 1);
		}

		//calculate the SH scaling constant
		double K(int l, int m){

			// renormalisation constant for SH function
			double temp = ((2.0*l+1.0)*factorial(l-m)) / (4.0*PI*factorial(l+m));
			return sqrt(temp);
		}

		//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 a point sample of a Spherical Harmonic basis function
			// l is the band, range [0..N]
			// m in the range [-l..l]
			// theta in the range [0..Pi]
			// phi in the range [0..2*Pi]
			const double sqrt2 = sqrt(2.0);
			if(m==0) return K(l,0)*P(l,m,cos(theta));
			else if(m>0) return sqrt2*K(l,m)*cos(m*phi)*P(l,m,cos(theta));
			else return sqrt2*K(l,-m)*sin(-m*phi)*P(l,-m,cos(theta));
		}

		void gen_function(){

			//initialize the function to zero
			memset(func, 0, sizeof(double) * N * N);

			double theta, phi;
			double result;
			int l, m;

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

				for(unsigned int xi = 0; xi < N; xi++)
					for(unsigned int yi = 0; yi < N; yi++){

						theta = (2 * PI) * ((double)xi / (N-1));
						phi = PI * ((double)yi / (N-1));
						result = C[c] * SH(l, m, phi, theta);		//phi and theta are reversed here (damn physicists)
						func[yi * N + xi] += result;
					}

				m++;			//increment m

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

		void gl_prep_draw(){

			//enable depth testing
				//this has to be used instead of culling because the sphere can have negative values
			glEnable(GL_DEPTH_TEST);
			glDepthMask(GL_TRUE);
			glEnable(GL_TEXTURE_2D);	//enable 2D texture mapping
		}

		//draw a texture mapped sphere representing the function surface
		void gl_draw_sphere() {

			//PI is used to convert from spherical to cartesian coordinates
			//const double PI = 3.14159;

			//bind the 2D texture representing the color map
			glBindTexture(GL_TEXTURE_2D, color_tex);

			//Draw the Sphere
			int i, j;

			for(i = 1; i <= N-1; i++) {
				double phi0 = PI * ((double) (i - 1) / (N-1));
				double phi1 = PI * ((double) i / (N-1));

				glBegin(GL_QUAD_STRIP);
				for(j = 0; j <= N; j++) {

					//calculate the indices into the function array
					int phi0_i = i-1;
					int phi1_i = i;
					int theta_i = j;
					if(theta_i == N)
						theta_i = 0;

					double v0 = func[phi0_i * N + theta_i];
					double v1 = func[phi1_i * N + theta_i];

					v0 = fabs(v0);
					v1 = fabs(v1);


					double theta = 2 * PI * (double) (j - 1) / N;
					double x0 = v0 * cos(theta) * sin(phi0);
					double y0 = v0 * sin(theta) * sin(phi0);
					double z0 = v0 * cos(phi0);

					double x1 = v1 * cos(theta) * sin(phi1);
					double y1 = v1 * sin(theta) * sin(phi1);
					double z1 = v1 * cos(phi1);

					glTexCoord2f(theta / (2 * PI), phi0 / PI);
					glVertex3f(x0, y0, z0);

					glTexCoord2f(theta / (2 * PI), phi1 / PI);
					glVertex3f(x1, y1, z1);
				}
				glEnd();
			}
		}

		//draw a wire frame sphere representing the function surface
		void gl_draw_wireframe() {

					//PI is used to convert from spherical to cartesian coordinates
					//const double PI = 3.14159;

					//bind the 2D texture representing the color map
					glDisable(GL_TEXTURE_2D);
					glColor3f(0.0f, 0.0f, 0.0f);

					//Draw the Sphere
					int i, j;

					for(i = 1; i <= N-1; i++) {
						double phi0 = PI * ((double) (i - 1) / (N-1));
						double phi1 = PI * ((double) i / (N-1));

						glBegin(GL_LINE_STRIP);
						for(j = 0; j <= N; j++) {

							//calculate the indices into the function array
							int phi0_i = i-1;
							int phi1_i = i;
							int theta_i = j;
							if(theta_i == N)
								theta_i = 0;

							double v0 = func[phi0_i * N + theta_i];
							double v1 = func[phi1_i * N + theta_i];

							v0 = fabs(v0);
							v1 = fabs(v1);


							double theta = 2 * PI * (double) (j - 1) / N;
							double x0 = WIRE_SCALE * v0 * cos(theta) * sin(phi0);
							double y0 = WIRE_SCALE * v0 * sin(theta) * sin(phi0);
							double z0 = WIRE_SCALE * v0 * cos(phi0);

							double x1 = WIRE_SCALE * v1 * cos(theta) * sin(phi1);
							double y1 = WIRE_SCALE * v1 * sin(theta) * sin(phi1);
							double z1 = WIRE_SCALE * v1 * cos(phi1);

							glTexCoord2f(theta / (2 * PI), phi0 / PI);
							glVertex3f(x0, y0, z0);

							glTexCoord2f(theta / (2 * PI), phi1 / PI);
							glVertex3f(x1, y1, z1);
						}
						glEnd();
					}
				}

		void init(unsigned int n){

			//set the sphere resolution
			N = n;

			//allocate space for the color map
			unsigned int bytes = N * N * sizeof(unsigned char) * 3;
			unsigned char* color_image;
			color_image = (unsigned char*) malloc(bytes);

			//allocate space for the function
			func = (double*) malloc(N * N * sizeof(double));

			//generate a function (temporary)
			gen_function();

			//generate a colormap from the function
			stim::cpu2cpu<double>(func, color_image, N*N, stim::cmBrewer);

			//prep everything for drawing
			gl_prep_draw();			

			//generate an OpenGL texture map in the current context
			glGenTextures(1, &color_tex);
			//bind the texture
			glBindTexture(GL_TEXTURE_2D, color_tex);

			//copy the color data from the buffer to the GPU
			glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, N, N, 0, GL_RGB, GL_UNSIGNED_BYTE, color_image);

			//initialize all of the texture parameters
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
			glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
			
			//free the buffer
			free(color_image);
		}


	public:

		void glRender(){
			//set all OpenGL parameters required for drawing
			gl_prep_draw();

			//draw the sphere
			gl_draw_sphere();
			//gl_draw_wireframe();

		}

		void glInit(unsigned int n){
			init(n);
		}

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





	};		//end class sph_harmonics




}


#endif