#include <complex.h>
#include <math.h>
#include <stdlib.h>

#define CGLM_ALL_UNALIGNED
#include <cglm/vec3.h>
#include <cglm/vec4.h>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#ifndef CMPLX
#define CMPLX(a,b) (a+I*b)
#endif

void riemman(float *d_surface, int * coords, int grid_size)
{
	complex double eq;
	float u = 2 * ((float)coords[0] / grid_size) - 1;
	float v = 2 * ((float)coords[1] / grid_size) - 1;

	eq = csqrt(CMPLX(u,v));

	d_surface[0] = u;
	d_surface[1] = v;
	d_surface[2] = creal(eq);
	d_surface[3] = cimag(eq);
}


void cube( float *d_surface, int * coord, int grid_size )
{
	unsigned char i;

	for(int i=0; i<4; i++ )
		d_surface[i]=(float)coord[i]/grid_size;
}

void mobius(float *d_surface, int * coord, int grid_size)
{
	const float width = 0.5;
	float u = (2 * M_PI) * ((float)coord[0] / grid_size);
	float v = (2 * width) * ((float)coord[1] / grid_size) - width;

	d_surface[0] = cos(u) + v * cos(u / 2) * cos(u);
	d_surface[1] = sin(u) + v * cos(u / 2) * sin(u);
	d_surface[2] = v * sin(u / 2);
}

void torus(float *d_surface, int * coord, int grid_size)
{
	float u = (2 * M_PI) * ((float)coord[0] / grid_size);
	float v = (2 * M_PI) * ((float)coord[1] / grid_size);

	d_surface[0] = (1 + 0.5 * cos(v)) * cos(u);
	d_surface[1] = (1 + 0.5 * cos(v)) * sin(u);
	d_surface[2] = 0.5 * sin(v);
}

void klein(float *d_surface, int * coord, int grid_size)
{
	float u = (2 * M_PI) * ((float)coord[0] / grid_size);
	float v = (2 * M_PI) * ((float)coord[1]/ grid_size);

	d_surface[0] = (0.5 * cos(v) + 0.5) * cos(u);
	d_surface[1] = (0.5 * cos(v) + 0.5) * sin(u);
	d_surface[2] = sin(v) * cos(u / 2);
	d_surface[3] = sin(v) * sin(u / 2);
}

typedef void (*function_t)(float *, int *, int);

float *generate_data_surface(int grid_size, unsigned char *s)
{
	unsigned int i, j, k, o, p, l, n, m;
	long size, q=0;
	function_t f;
	float *d_surface;

	const int dim =2;
	int cara[dim];
	char bits[dim+1];
	bits[dim]=0;

	f =klein ;
	*s = 4;

	size = grid_size * grid_size * 6 * (*s) * 24;
	d_surface = malloc((size + 1) * sizeof(float));
	d_surface[0] = size;
	
	for(o = 0; o < dim; o ++)
	{  
		for (p = 0; p < o; p++)
		{
			for (k = 0; k < (1 << (dim-2)); k++)
			{
				unsigned char skip=0;
				for(n = 0; n < dim-2; n++)
				{
					if( n==(o-1) || n==p )
						skip++;

					cara[n+skip] = (k & (1<<n))?grid_size:0;
				}

				for(i = 0; i < grid_size; i++)
				{
					for (j = 0; j < grid_size; j++)
					{
						cara[o] = i;
						cara[p] = j;
						f(&d_surface[q + 1], cara, grid_size);
						q += *s;

						cara[o] = i + 1;
						cara[p] = j;
						f(&d_surface[q + 1], cara, grid_size);
						q += *s;

						cara[o] = i + 1;
						cara [p] = j + 1;
						f(&d_surface[q + 1], cara, grid_size);
						q += *s;

						cara[o] = i;
						cara [p] = j; 
						f(&d_surface[q + 1], cara, grid_size);
						q += *s;

						cara[o] = i;
						cara [p] = j + 1;
						f(&d_surface[q + 1], cara, grid_size);
						q += *s;

						cara[o] = i + 1;
						cara [p] = j + 1;
						f(&d_surface[q + 1], cara, grid_size);
						q += *s;
					}
				}
			}	
		}
	}

	return d_surface;
}



static void __calculate_normal(
	float *p1, float *p2, float *p3, float *normal, unsigned char n)
{
	float alpha;
	vec4 v1, v2, v3;
	vec4 u1, u2, u3;

	switch (n)
	{
	case 3:

		glm_vec3_sub(p2, p1, v1);
		glm_vec3_sub(p3, p1, v2);

		glm_vec3_cross(v1, v2, normal);
		glm_vec3_normalize(normal);
		return;

	case 4:

		/*
			In Grant-Shmidth we need 3 linearly independian vector that forms a
		   basis, so we can have a ortonormal version of that basis, since, we
		   must have v1 = p3 - p1 v2 = p2 - p1 Then v3 = p1, will most certantly
		   be linerly independiant to v1 and v2.
		*/
		glm_vec4_sub(p2, p1, v1);
		glm_vec4_sub(p3, p1, v2);
		glm_vec4_copy(p1, v3);

		/* Setup U1 */
		{
			glm_vec4_copy(v1, u1);
		}

		/* Setup U2 */
		{
			vec4 proj;

			alpha = glm_vec4_dot(v2, u1) / glm_vec4_dot(u1, u1);
			glm_vec4_scale(u1, alpha, proj);

			glm_vec4_sub(v2, proj, u2);
		}

		/* Setup U3 */
		{
			vec4 proj1, proj2;

			alpha = glm_vec4_dot(v3, u1) / glm_vec4_dot(u1, u1);
			glm_vec4_scale(u1, alpha, proj1);

			alpha = glm_vec4_dot(v3, u2) / glm_vec4_dot(u2, u2);
			glm_vec4_scale(u2, alpha, proj2);

			glm_vec4_sub(v3, proj1, u3);
			glm_vec4_sub(u3, proj2, u3);
		}

		glm_vec4_copy(u3, normal);
		glm_vec4_normalize(normal);
		return;
	}
}

float *generate_normals_surface(float *d, unsigned char m)
{
	float *n;
	n = malloc((*d + 1) * sizeof(float));
	*n = *d;
	for (int i = 0; i < *d; i += 3 * m)
	{

		vec4 norm_vec;

		__calculate_normal(
			(d + 1) + i, (d + 1) + i + m, (d + 1) + i + 2 * m, norm_vec, m);
		glm_vec3_copy(norm_vec, (n + 1) + i);
		glm_vec3_copy(norm_vec, (n + 1) + i + m);
		glm_vec3_copy(norm_vec, (n + 1) + i + 2 * m);
	}
	return n;
}