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

#include "SurfaceGeometry.h"

#define TWOPI           6.283185307179586476925287
#define PI              3.141592653589793238462643

typedef struct {GLfloat x,y,z;} recVec;
typedef struct {GLfloat r,g,b;} recColor;
float *sin1, *cos1, *sin2, *cos2, *sin3, *cos3;
long maxI, maxJ, uvRatio; double dp;
/* Code based on work by Paul Bourke */
void initTabs(unsigned int subdivisions, unsigned int xyRatio){
   maxI = subdivisions * xyRatio; maxJ = subdivisions;
   uvRatio = xyRatio;
   dp = TWOPI/maxI;
   float * grab(){return (float *)malloc((maxI+1)*sizeof(float));}
   sin1 = grab(), cos1 = grab(),  sin2 = grab(), cos2 = grab(),
     sin3 = grab(), cos3 = grab();
   {int j = maxI+1; while(j--) {sin1[j] = sin(dp*j); cos1[j] = cos(dp*j);
      sin2[j] = sin(2*dp*j); cos2[j] = cos(2*j*dp);
	  sin3[j] = sin(3*dp*j); cos3[j] = cos(3*dp*j);}};}

typedef struct{double x, y, z;} vc;
void pr(vc x, char* m){printf("%s = %12.9f %12.9f %12.9f\n", m, x.x, x.y, x.z);}
vc vsub(vc a, vc b){return(vc){a.x - b.x, a.y - b.y, a.z - b.z};}
double ip(vc p, vc q){return p.x*q.x + p.y*q.y + p.z*q.z;}
vc sp(double x, vc v){return (vc){x*v.x, x*v.y, x*v.z};}
vc nm(vc q){return sp(1/sqrt(ip(q, q)), q);}

recVec sh(vc x){return (recVec){x.x, x.y, x.z};}
typedef struct{recVec p, n;} sur;
static sur Eval(int u, int v, double g, double gd) {
	vc vadd(vc a, vc b){return(vc){a.x + b.x, a.y + b.y, a.z + b.z};}
	vc X(vc a, vc b){return (vc){a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x};}
    vc r = (vc){cos2[u] - cos1[u]/2, sin2[u] + sin1[u]/2, sin3[u]/2};
	vc Rd = (vc){-2*sin2[u] + sin1[u]/2, 2*cos2[u] + cos1[u]/2, 3*cos3[u]};
	double spd2 = ip(Rd, Rd), spd = sqrt(spd2);
	vc rdd = (vc){-4*cos2[u] + cos1[u]/2, -4*sin2[u] - sin1[u]/2, -9*sin3[u]};
    vc rd = sp(1/spd, Rd);
    vc rn = nm(vsub(rdd, sp(ip(rdd, rd), rd)));
    vc stk = vadd(sp(sin1[v], rn), sp(cos1[v], X(rd, rn)));
	return (sur){sh(sp(2.5, vadd(r, sp(g, stk)))),
	     sh(nm(vadd(sp(spd - g*ip(rdd, rn)/spd2*sin1[v], stk), sp(-gd, rd))))};
}

void BuildGeometry (){
	int i;
	recVec rowpos[maxJ], rownorm[maxJ]; recColor memcol;
	static float tim=0;++tim;
	for (i=0; i<= maxI; i++) {
		recColor c = (recColor){sin1[i], sin1[(i+maxI/3)%maxI], sin1[(i+maxI/3*2)%maxI]};
	    float const b = 6, h=.5;
	    double x = (fmod((.04*tim+i*dp), TWOPI)-PI)*b;
		double const bmpSiz = .25;
	    double g = bmpSiz*(.7+h*exp(-x*x)), gd = bmpSiz*(-2*x*exp(-x*x)*b*h);
	//	printf("%d g=%e gd=%e\n", i, g, gd);
  		glBegin(GL_TRIANGLE_STRIP);
		{int j; for (j = 0; j <= maxJ; j++) {
			sur dtm = Eval(i, j*uvRatio, g, gd);
			if(i) {
				glColor3fv (&memcol.r);
				glNormal3fv (&rownorm[j].x);
				glVertex3fv (&rowpos[j].x);
			    glColor3fv (&c.r);
				glNormal3fv (&dtm.n.x);
				glVertex3fv (&dtm.p.x);}
			rownorm[j] = dtm.n; rowpos[j] = dtm.p;}}
		memcol = c;
		glEnd ();}
	glutPostRedisplay();
}