/*
 * $Id$
 *
 * $Date$
 * $Revision$
 *
 * (C) 1999 by Hyperion
 * All rights reserved
 *
 * This file is part of the MiniGL library project
 * See the file Licence.txt for more details
 *
 */

#include <math.h>
#include <stdio.h>
#include "sysinc.h"
#include "mgl/minigl.h"

#define VEC_NORM(v)                                                       \
{                                                                         \
	GLfloat recM = fast_reciprocal_sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]); \
	v[0] *= recM;                                                           \
	v[1] *= recM;                                                           \
	v[2] *= recM;                                                           \
}

#define VEC_CROSS(v, a, b)                             \
	v[0] = a[1] * b[2] - a[2] * b[1];                    \
	v[1] = a[2] * b[0] - a[0] * b[2];                    \
	v[2] = a[0] * b[1] - a[1] * b[0];

#define VEC_SUB(v, a, b)                               \
	v[0] = a[0] - b[0];                                  \
	v[1] = a[1] - b[1];                                  \
	v[2] = a[2] - b[2];

#define VEC_ADD(v, a, b)                               \
	v[0] = a[0] + b[0];                                  \
	v[1] = a[1] + b[1];                                  \
	v[2] = a[2] + b[2];

#define VEC_PRINT(v) printf("<%f, %f, %f>\n", v[0], v[1], v[2]);

INLINE GLboolean mInvertMatrix(const GLdouble *matrixIn, GLdouble *matrixOut);
INLINE void mMatMultGeneral(const GLdouble *pA, const GLdouble *pB, GLdouble *pC);

static char rcsid[] UNUSED = "$Id$ ";

const GLubyte * GLUErrorString(GLenum err)
{
	switch(err) {
		case GL_NO_ERROR:          return (const GLubyte * )"no error";
		case GL_INVALID_ENUM:      return (const GLubyte * )"invalid enum";
		case GL_INVALID_VALUE:     return (const GLubyte * )"invalid value";
		case GL_INVALID_OPERATION: return (const GLubyte * )"invalid operation";
		case GL_STACK_OVERFLOW:    return (const GLubyte * )"stack overflow";
		case GL_STACK_UNDERFLOW:   return (const GLubyte * )"stack underflow";
		case GL_OUT_OF_MEMORY:     return (const GLubyte * )"out of memory";
		default:                   return (const GLubyte * )"unknown error";
	}
}

void GLULookAt(GLfloat ex, GLfloat ey, GLfloat ez, GLfloat cx, GLfloat cy, GLfloat cz, GLfloat ux, GLfloat uy, GLfloat uz)
{
	GLfloat u[3], v[3], w[3];
	GLfloat m[16];

	w[0] = ex - cx;     w[1] = ey - cy;     w[2] = ez - cz;
	v[0] = ux;          v[1] = uy;          v[2] = uz;

	VEC_NORM(w);
	VEC_CROSS(u, v, w);
	VEC_CROSS(v, w, u);
	VEC_NORM(u);
	VEC_NORM(v);

	m[ 0] = u[0];   m[ 1] = v[0];   m[ 2] = w[0];   m[ 3] = 0.0;
	m[ 4] = u[1];   m[ 5] = v[1];   m[ 6] = w[1];   m[ 7] = 0.0;
	m[ 8] = u[2];   m[ 9] = v[2];   m[10] = w[2];   m[11] = 0.0;
	m[12] = 0.0;    m[13] = 0.0;    m[14] = 0.0;    m[15] = 1.0;

	glMultMatrixf(m);
	glTranslatef(-ex, -ey, -ez);
}

void GLUPerspective(GLfloat fovy, GLfloat aspect, GLfloat znear, GLfloat zfar)
{
	GLfloat xmin, xmax, ymin, ymax;

	ymax = znear * tan(fovy * 0.008726646);
	ymin = -ymax;
	xmin = ymin * aspect;
	xmax = ymax * aspect;

	glFrustum(xmin, xmax, ymin, ymax, znear, zfar);
}

void GLUOrtho2D(GLcontext context, GLdouble left, GLdouble right, GLdouble bottom, GLdouble top)
{
	glOrtho(left, right, bottom, top, -1, 1);
}

GLint GLUBuild2DMipmaps(
	GLcontext ctx, GLenum target, GLint components, GLint width, GLint height, GLenum format, GLenum type, const void *data
)
{
	dprintf("*** Implement me\n");
	glTexImage2D(target, 0, components, width, height, 0, format, type, data);

	return 0;
}

GLint GLUProject(
	GLcontext ctx, GLdouble objX, GLdouble objY, GLdouble objZ, const GLdouble *model, const GLdouble *proj, const GLint *view,
	GLdouble* winX, GLdouble* winY, GLdouble* winZ
)
{
	GLdouble fullMatrix[16], objT[3];
	mMatMultGeneral(proj, model, fullMatrix);

	#define a(x) (fullMatrix[OF_##x])

	objT[0] = a(11)*objX + a(12)*objY + a(13)*objZ + a(14);
	objT[1] = a(21)*objX + a(22)*objY + a(23)*objZ + a(24);
	objT[2] = a(31)*objX + a(32)*objY + a(33)*objZ + a(34);

	*winX = view[0] + view[2] * (objT[0] + 1) / 2;
	*winY = view[1] + view[3] * (objT[1] + 1) / 2;
	*winZ = (objT[2] + 1) / 2;

	#undef a

	return GL_TRUE;
}

GLint GLUUnProject(
	GLcontext ctx, GLdouble winX, GLdouble winY, GLdouble winZ, const GLdouble *model, const GLdouble *proj, const GLint *view,
	GLdouble*	objX, GLdouble*	objY, GLdouble*	objZ
)
{
	GLdouble fullMatrix[16], inverseMatrix[16], winT[3];
	mMatMultGeneral(proj, model, fullMatrix);
	if (mInvertMatrix(fullMatrix, inverseMatrix) == GL_FALSE) {
		return GL_FALSE;
	}

	if (view[2] == 0 || view[3] == 0) {
		return GL_FALSE;
	}
	winT[0] = 2 * (winX - view[0]) / view[2] - 1;
	winT[1] = 2 * (winY - view[1]) / view[3] - 1;
	winT[2] = 2 * winZ - 1;

	#define a(x) (inverseMatrix[OF_##x])
	*objX = a(11)*winT[0] + a(12)*winT[1] + a(13)*winT[2] + a(14);
	*objY = a(21)*winT[0] + a(22)*winT[1] + a(23)*winT[2] + a(24);
	*objZ = a(31)*winT[0] + a(32)*winT[1] + a(33)*winT[2] + a(34);

	#undef a

	return GL_TRUE;
}

void GLUPickMatrix(GLcontext ctx, GLdouble x, GLdouble y, GLdouble delX, GLdouble delY, GLint *viewport)
{
	GLfloat scaleX, scaleY;
	GLfloat translateX, translateY;
	GLfloat m[16];

	// Need to scale and move the viewport to the specified region
	scaleX = viewport[2] / delX;
	scaleY = viewport[3] / delY;
	translateX = scaleX + (viewport[0] - x) * 2.0 / delX;
	translateY = scaleY + (viewport[1] - y) * 2.0 / delY;

	m[OF_11] = scaleX;   m[OF_12] = 0.0;      m[OF_13] = 0.0;    m[OF_14] = translateX;
	m[OF_21] = 0.0;      m[OF_22] = scaleY;   m[OF_23] = 0.0;    m[OF_24] = translateY;
	m[OF_31] = 0.0;      m[OF_32] = 0.0;      m[OF_33] = 1.0;    m[OF_34] = 0.0;
	m[OF_41] = 0.0;      m[OF_42] = 0.0;      m[OF_43] = 0.0;    m[OF_44] = 1.0;

	glMultMatrixf(m);
}

// ---- Matrix functions ---- Required for gluProject()/gluUnProject

/*
 * Compute determinant of 3x3 matrix resulting from pA dropping
 * column i and line j
 * (Adapted from the Matrix and Quaternion FAQ, at
 * http://skal.planet-d.net/demo/matrixfaq.htm)
 */
INLINE GLdouble mDeterminant3(const GLdouble *pA, int i, int j)
{
	GLdouble a[9];
	int di, dj, si, sj;

	for (di = 0; di < 3; di++) {
		for (dj = 0; dj < 3; dj++) {
			si = di + ((di >= i) ? 1 : 0);
			sj = dj + ((dj >= j) ? 1 : 0);
			a[di * 3 + dj] = pA[sj * 4 + si]; // really ?
		}
	}

	return a[0] * ( a[4] * a[8] - a[7] * a[5])
	     - a[1] * ( a[3] * a[8] - a[6] * a[5])
	     + a[2] * ( a[3] * a[7] - a[6] * a[4]);
}


INLINE void mLoadIdentity(GLdouble *pA)
{
	#define a(x) pA[OF_##x] = 0.f;
	#define b(x) pA[OF_##x] = 1.f;

	b(11); a(21); a(31); a(41);
	a(12); b(22); a(32); a(42);
	a(13); a(23); b(33); a(43);
	a(14); a(24); a(34); b(44);

	#undef a
	#undef b
}

INLINE GLdouble mDeterminant4(const GLdouble *pA)
{
	GLdouble res = 0.0, i = 1;
	int n;

	for (n = 0; n < 4; n++) {
		res += mDeterminant3(pA, 0, n) * i;
		i *= -1.0;
	}

	return res;
}

INLINE GLboolean mInvertMatrix(const GLdouble *matrixIn, GLdouble *matrixOut)
{
	int i, j, sign;

	GLdouble mdet = mDeterminant4(matrixIn);

	if (fast_fabs(mdet) < 0.5e-9) {
		return GL_FALSE;
	}
	else {
		for (i = 0; i < 4; i++) {
			for (j = 0; j < 4; j++) {
				sign = 1 - ((i + j) % 2) * 2;
				matrixOut[i*4+j] = (mDeterminant3(matrixIn, i, j) * (GLdouble)sign) / mdet;
			}
		}
	}

	return GL_TRUE;
}

/*
** General case: Matrix multiply
** Multiply matrix A by float field, storing the result in matrix C
*/
INLINE void mMatMultGeneral(const GLdouble *pA, const GLdouble *pB, GLdouble *pC)
{
	#define a(x) (pA[OF_##x])
	#define b(x) (pB[OF_##x])
	#define c(x) (pC[OF_##x])

	GLdouble a11 = a(11),  b11 = b(11);
	GLdouble a12 = a(12),  b12 = b(12);
	GLdouble a13 = a(13),  b13 = b(13);
	GLdouble a14 = a(14),  b14 = b(14);
	GLdouble a21 = a(21),  b21 = b(21);
	GLdouble a22 = a(22),  b22 = b(22);
	GLdouble a23 = a(23),  b23 = b(23);
	GLdouble a24 = a(24),  b24 = b(24);
	GLdouble a31 = a(31),  b31 = b(31);
	GLdouble a32 = a(32),  b32 = b(32);
	GLdouble a33 = a(33),  b33 = b(33);
	GLdouble a34 = a(34),  b34 = b(34);
	GLdouble a41 = a(41),  b41 = b(41);
	GLdouble a42 = a(42),  b42 = b(42);
	GLdouble a43 = a(43),  b43 = b(43);
	GLdouble a44 = a(44),  b44 = b(44);


	c(11) = a11*b11 + a12*b21 + a13*b31 + a14*b41;
	c(12) = a11*b12 + a12*b22 + a13*b32 + a14*b42;
	c(13) = a11*b13 + a12*b23 + a13*b33 + a14*b43;
	c(14) = a11*b14 + a12*b24 + a13*b34 + a14*b44;
	c(21) = a21*b11 + a22*b21 + a23*b31 + a24*b41;
	c(22) = a21*b12 + a22*b22 + a23*b32 + a24*b42;
	c(23) = a21*b13 + a22*b23 + a23*b33 + a24*b43;
	c(24) = a21*b14 + a22*b24 + a23*b34 + a24*b44;
	c(31) = a31*b11 + a32*b21 + a33*b31 + a34*b41;
	c(32) = a31*b12 + a32*b22 + a33*b32 + a34*b42;
	c(33) = a31*b13 + a32*b23 + a33*b33 + a34*b43;
	c(34) = a31*b14 + a32*b24 + a33*b34 + a34*b44;
	c(41) = a41*b11 + a42*b21 + a43*b31 + a44*b41;
	c(42) = a41*b12 + a42*b22 + a43*b32 + a44*b42;
	c(43) = a41*b13 + a42*b23 + a43*b33 + a44*b43;
	c(44) = a41*b14 + a42*b24 + a43*b34 + a44*b44;

	#undef a
	#undef b
	#undef c
}