#include "vogle.h"

#ifdef	TC

extern	double	cos();
extern	double	sin();
extern	double	asin();
extern	double	sqrt();
extern	double	fabs();

#else 

#include <math.h>

#endif

/*
 * NOTE: the words hither and yon are used in this file instead of near and far
 * as they are keywords on some PC compilers (groan). Change them back at your 
 * on peril.
 */

#define	SQ(a)	((a)*(a))
#define COT(a)	((float)(cos((double)(a)) / sin((double)(a))))

/*
 * polarview
 *
 * Specify the viewer's position in polar coordinates by giving
 * the distance from the viewpoint to the world origin, 
 * the azimuthal angle in the x-y plane, measured from the y-axis,
 * the incidence angle in the y-z plane, measured from the z-axis,
 * and the twist angle about the line of sight.
 *
 */
void polarview(float dist, float azim, float inc, float twist)
{
  if (!vdevice.initialised) 
    verror("polarview: vogle not initialised");
  
  translate(0.0, 0.0, -dist);
  rotate(-twist, 'z');
  rotate(-inc, 'x');
  rotate(-azim, 'z');
}

/*
 * up
 *
 *	set the up vector
 */
void up(float x, float y, float z)
{
  vdevice.upset = 1;
  
  vdevice.upvector[V_X] = x;
  vdevice.upvector[V_Y] = y;
  vdevice.upvector[V_Z] = z;
}

/*
 * normallookat
 *
 *	do the standard lookat transformation.
 */
static void normallookat(float vx, float vy, float vz, 
			 float px, float py, float pz)
{
  float	l2, l3, sintheta, sinphi, costheta, cosphi;
  Matrix	tmp;
  
  l2 = sqrt((double)(SQ((px - vx)) + SQ((pz - vz))));
  l3 = sqrt((double)(SQ((px - vx)) + SQ((py - vy)) + SQ((pz - vz))));
  
  if (l3 != 0.0) {
    sinphi = (vy - py) / l3;
    cosphi = l2 / l3;
    
    /*
     * Rotate about X by phi
     */
    identmatrix(tmp);
    tmp[1][1] = tmp[2][2] = cosphi;
    tmp[1][2] = sinphi;
    tmp[2][1] = -sinphi;
    multmatrix(tmp);
  }
  
  if (l2 != 0.0) {
    sintheta = (px - vx) / l2;
    costheta = (vz - pz) / l2;
    
    /*
     * Rotate about Y by theta
     */
    identmatrix(tmp);
    tmp[0][0] = tmp[2][2] = costheta;
    tmp[0][2] = -sintheta;
    tmp[2][0] = sintheta;
    multmatrix(tmp);
  }
}

/*
 * lookatwithup
 *
 *	do the standard lookat transformation using an up vector as well.
 */
static void lookatwithup(float vx, float vy, float vz, 
			 float px, float py, float pz)
{
  Vector	t, u, s;
  Matrix	tmp;
  double	dy, dz, lt, lu;
  
  t[V_X] = vx - px;
  t[V_Y] = vy - py;
  t[V_Z] = vz - pz;
  
  u[V_X] = vdevice.upvector[V_X];
  u[V_Y] = vdevice.upvector[V_Y];
  u[V_Z] = vdevice.upvector[V_Z];
  
  lt = sqrt(t[V_X] * t[V_X] + t[V_Y] * t[V_Y] + t[V_Z] * t[V_Z]);
  
  if (lt == 0.0)
    verror("vogle: can't have eyepoint and reference point the same in lookat.\n");
  
  lu = sqrt(u[V_X] * u[V_X] + u[V_Y] * u[V_Y] + u[V_Z] * u[V_Z]);
  
  if (lu == 0.0)
    verror("vogle: invalid up vector in lookat.\n");
  
  /*
   * normalise t and u
   */
  t[V_X] /= lt;
  t[V_Y] /= lt;
  t[V_Z] /= lt;
  
  u[V_X] /= lu;
  u[V_Y] /= lu;
  u[V_Z] /= lu;
  
  dz = t[V_X] * u[V_X] + t[V_Y] * u[V_Y] + t[V_Z] * u[V_Z];
  
  if (fabs(dz) >= 1.0)
    verror("vogle: up vector and direction of view are the same.\n");
  
  dy = sqrt(1.0 - dz * dz);
  
  /*
   * calculate view up
   */
  u[V_X] = (u[V_X] - dz * t[V_X]) / dy;
  u[V_Y] = (u[V_Y] - dz * t[V_Y]) / dy;
  u[V_Z] = (u[V_Z] - dz * t[V_Z]) / dy;
  
  /*
   * calculate side vector (cross product of u and t)
   */
  s[V_X] = u[V_Y] * t[V_Z] - u[V_Z] * t[V_Y];
  s[V_Y] = u[V_Z] * t[V_X] - u[V_X] * t[V_Z];
  s[V_Z] = u[V_X] * t[V_Y] - u[V_Y] * t[V_X];
  
  identmatrix(tmp);
  
  tmp[0][0] = s[V_X];
  tmp[1][0] = s[V_Y];
  tmp[2][0] = s[V_Z];
  
  tmp[0][1] = u[V_X];
  tmp[1][1] = u[V_Y];
  tmp[2][1] = u[V_Z];
  
  tmp[0][2] = t[V_X];
  tmp[1][2] = t[V_Y];
  tmp[2][2] = t[V_Z];
  
  multmatrix(tmp);
}

/*
 * lookat
 *
 * Specify the viewer's position by giving a viewpoint and a 
 * reference point in world coordinates. A twist about the line
 * of sight may also be given. 
 */
void lookat(float vx, float vy, float vz, 
	    float px, float py, float pz, float twist)
{
  if (!vdevice.initialised) 
    verror("lookat: vogle not initialised");
  
  rotate(-twist, 'z');
  
  if (vdevice.upset)
    lookatwithup(vx, vy, vz, px, py, pz);
  else
    normallookat(vx, vy, vz, px, py, pz);
  
  translate(-vx, -vy, -vz);
}

/*
 * perspective
 *
 * Specify a perspective viewing pyramid in world coordinates by
 * giving a field of view, aspect ratio, and the locations of the 
 * near(hither) and far(yon) clipping planes in the z direction.
 */
void perspective(float fov, float aspect, float hither, float yon)
{
  Matrix		mat;
  
  if (!vdevice.initialised)
    verror("perspective: vogle not initialised");
  
  if (aspect == 0.0)
    verror("perspective: can't have zero aspect ratio!");
  
  if ((yon - hither) == 0.0)
    verror("perspective: near clipping plane same as far one.");
  
  if (fov == 0.0 || fov == 180.0)
    verror("perspective: bad field of view passed.");
  
  identmatrix(mat);
  
  mat[0][0] = COT((D2R * fov / 2.0)) / aspect;
  mat[1][1] = COT((D2R * fov / 2.0));
  
  mat[2][2] = -(yon + hither) / (yon - hither);
  mat[2][3] = -1;
  mat[3][2] = -2.0 * yon * hither / (yon - hither);
  mat[3][3] = 0;
  
  loadmatrix(mat);
}

/*
 * window
 *
 * Specify a perspective viewing pyramid in world coordinates by
 * giving a rectangle at the near clipping plane and the location
 * of the far clipping plane.
 *
 */
void window(float left, float right, float bottom, float top, 
	    float hither, float yon)
{
  Matrix		mat;
  
  if (!vdevice.initialised)
    verror("window: vogle not initialised");
  
  if ((right - left) == 0.0)
    verror("window: left clipping plane same as right one.");
  
  if ((top - bottom) == 0.0)
    verror("window: bottom clipping plane same as top one.");
  
  if ((yon - hither) == 0.0)
    verror("window: near clipping plane same as far one.");
  
  identmatrix(mat);
  
  mat[0][0] = 2.0 * hither / (right - left);
  mat[1][1] = 2.0 * hither / (top - bottom);
  mat[2][0] = (right + left) / (right - left);
  mat[2][1] = (top + bottom) / (top - bottom);
  mat[2][2] = -(yon + hither) / (yon - hither);
  mat[2][3] = -1.0;
  mat[3][2] = -2.0 * yon * hither / (yon - hither);
  mat[3][3] = 0.0;
  
  loadmatrix(mat);
}

/*
 * ortho
 *
 * Define a three dimensional viewing box by giving the left,
 * right, bottom and top clipping plane locations and the distances
 * along the line of sight to the near and far clipping planes.
 *
 */
void ortho(float left, float right, float bottom, float top, 
	   float hither, float yon)
{
  Matrix		mat;
  
  if (!vdevice.initialised)
    verror("ortho: vogle not initialised");
  
  if ((right - left) == 0.0)
    verror("ortho: left clipping plane same as right one.");
  
  if ((top - bottom) == 0.0)
    verror("ortho: bottom clipping plane same as top one.");
  
  if ((yon - hither) == 0.0)
    verror("ortho: near clipping plane same as far one.");
  
  identmatrix(mat);
  
  mat[0][0] = 2.0 / (right - left);
  mat[1][1] = 2.0 / (top - bottom);
  mat[2][2] = -2.0 / (yon - hither);
  mat[3][0] = -(right + left) / (right - left);
  mat[3][1] = -(top + bottom) / (top - bottom);
  mat[3][2] = -(yon + hither) / (yon - hither);
  
  loadmatrix(mat);
}

/*
 * ortho2
 *
 * Specify a two dimensional viewing rectangle. 
 *
 */
void ortho2(float left, float right, float bottom, float top)
{
  Matrix	mat;
  
  if (!vdevice.initialised) 
    verror("ortho2: vogle not initialised");
  
  identmatrix(mat);
  
  if ((right - left) == 0.0)
    verror("ortho2: left clipping plane same as right one.");
  
  if ((top - bottom) == 0.0)
    verror("ortho2: bottom clipping plane same as top one.");
  
  mat[0][0] = 2.0 / (right - left);
  mat[1][1] = 2.0 / (top - bottom);
  mat[2][2] = -1.0;
  mat[3][0] = -(right + left) / (right - left);
  mat[3][1] = -(top + bottom) / (top - bottom);
  
  loadmatrix(mat);
}
