#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 stim::spharmonics<T>::SH;
	stim::spharmonics<T> surface;	//harmonic that stores the surface information
	stim::spharmonics<T> color;	//harmonic that stores the color information
	T* func_surface;		//stores the raw function data (samples at each point)
	T* func_color;			//stores the raw color 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

	///generates the 
	void gen_surface(){
		//allocate memory and initialize the function to zero
		func_surface = (T*) malloc(N * N * sizeof(T));
		memset(func_surface, 0, sizeof(T) * N * N);

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


		l = m = 0;
		//for each coefficient
		for(unsigned int c = 0; c < surface.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 = surface.C[c] * SH(l, m, theta, phi);

					//store the result in a 2D array (which will later be used as a texture map)
					func_surface[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;
			}
		}
	}

	///generates the color function
	void gen_color(){

		//initialize the function to zero
		func_color = (T*) malloc(N * N * sizeof(T));
		memset(func_color, 0, sizeof(T) * N * N);

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


		l = m = 0;
		//for each coefficient
		for(unsigned int c = 0; c < color.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 = color.C[c] * SH(l, m, theta, phi);

					//store the result in a 2D array (which will later be used as a texture map)
					func_color[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_surface[phi0_i * N + theta_i];
				double v1 = func_surface[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();
		}

		glBindTexture(GL_TEXTURE_2D, 0);
	}

	void gen_texture()
	{
		//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);

		//generate a colormap from the function
		stim::cpu2cpu<double>(func_color, 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);
	}

	void gl_init(unsigned int n){

		//set the sphere resolution
		N = n;

		//set up the surface data.
		gen_surface();

		//set up the color data and generate the corresponding texture.
		gen_color();
		gen_texture();

	}

public:
	gl_spharmonics<T>(unsigned int n = 256) {
		gl_init(n);
	}

	gl_spharmonics<T>( stim::spharmonics<T> in_function, unsigned int n = 256)
	{
		surface = in_function;
		color = in_function;	
		gl_init(n);
	}

	gl_spharmonics<T>( stim::spharmonics<T> in_function, stim::spharmonics<T> in_color, unsigned int n = 256)
	{
		surface = in_function;
		color = in_color;	
		gl_init(n);
	}

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


};		//end namespace stim




#endif