#ifndef STIM_GL_SPHARMONICS_H
#define STIM_GL_SPHARMONICS_H

#include <GL/gl.h>

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

namespace stim{

	
template <typename T>
class gl_spharmonics : public spharmonics<T>{

protected:
	using spharmonics<T>::C;
	using spharmonics<T>::SH;
	
	T* 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

	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 each coefficient
		for(unsigned int c = 0; c < C.size(); c++){				

			//iterate through the entire 2D grid representing the function
			for(unsigned int xi = 0; xi < N; xi++){
				for(unsigned int yi = 0; yi < N; yi++){

					//get the spherical coordinates for each grid point
					theta = (2 * PI) * ((double)xi / (N-1));
					phi = PI * ((double)yi / (N-1));

					//sum the contribution of the current spherical harmonic based on the coefficient
					result = C[c] * SH(l, m, theta, phi);

					//store the result in a 2D array (which will later be used as a texture map)
					func[yi * N + xi] += result;
				}
			}

			//keep track of m and l here
			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() {

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

	void gl_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 functional representation that will be used for the texture map and vertices
		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:
	gl_spharmonics<T>() {

	}

	gl_spharmonics<T>(const spharmonics<T> copy) : gl_spharmonics<T>() {
		C = copy.C;
	}

	gl_spharmonics<T>& operator=(const spharmonics<T> rhs) {
		gl_spharmonics<T> result(rhs.C.size());
		result.C = spharmonics<T>::rhs.C;
		return result;
	}

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

		//draw the sphere
		gl_draw_sphere();
	}

	void glInit(unsigned int n = 256){
		gl_init(n);
		gen_function();
	}

	//overload arithmetic operations
	gl_spharmonics<T> operator*(T rhs) const {
		gl_spharmonics<T> result(C.size());                                //create a new sph    erical harmonics object
		for (size_t c = 0; c < C.size(); c++)                   //for each coefficient
			result.C[c] = C[c] * rhs;                                       //calculat    e the factor and store the result in the new spharmonics object
		return result;
	}

	gl_spharmonics<T> operator+(gl_spharmonics<T> rhs) {
		size_t low = std::min(C.size(), rhs.C.size());                          //store th    e number of coefficients in the lowest object
		size_t high = std::max(C.size(), rhs.C.size());                         //store th    e number of coefficients in the result
		bool rhs_lowest = false;                                                                    //true if rhs has the lowest number of coefficients
		if (rhs.C.size() < C.size()) rhs_lowest = true;                 //if rhs has a low    er number of coefficients, set the flag

		gl_spharmonics<T> result(high);                                                                //create a new object
		size_t c;
		for (c = 0; c < low; c++)                                               //perform     the first batch of additions
			result.C[c] = C[c] + rhs.C[c];                                          //    perform the addition

		for (c = low; c < high; c++) {
			if (rhs_lowest)
				result.C[c] = C[c];
			else
				result.C[c] = rhs.C[c];
		}
		return result;
	}

	gl_spharmonics<T> operator-(gl_spharmonics<T> rhs) {
		return (*this) + (rhs * (T)(-1));
	}



};		//end gl_spharmonics


};		//end namespace stim




#endif