/* ****************************************************************
 *                            geo.c
 * ****************************************************************
 *  MODULE PURPOSE:
 *      This module contains routines that preform common
 *      geometric manipulations for image synthesis.
 *
 *  MODULE CONTENTS:
 *      geo_dot         - vector dot product
 *      geo_cross       - vector cross product
 *      geo_norm        - vector normalization
 *      geo_line        - generate a line
 *      geo_rfl         - compute the reflected vector
 *      geo_rfr         - compute the refracted vector
 *      geo_H           - compute H vector
 *      geo_Ht          - compute the H' vector
 */
#include <stdio.h>
#include <math.h>
#include "geo.h"

/* ****************************************************************
 *  double geo_dot (v0, v1)
 *  DIR_VECT    *v0, *v1    (in) normalized direction vectors
 *
 *  Returns the dot product of two direction vectors, ref
 *   Eq.(\*(rK).  Note that the dot product is the cosine
 *   between the direction vectors because the direction
 *   vectors are normalized.
 */
double geo_dot (v0, v1)
DIR_VECT        *v0, *v1;
{   return (v0->i * v1->i) + (v0->j * v1->j) + (v0->k * v1->k);
}

/* ****************************************************************
 *  DIR_VECT geo_cross (v0, v1)
 *   DIR_VECT   *v0, *v1    (in)  normalized direction vectors
 *
 *  Returns the cross product of two direction vectors, refer
 *   to Eq.(\*(rL).
 */
DIR_VECT geo_cross (v0, v1)
DIR_VECT        *v0, *v1;
{   DIR_VECT        v2;
    v2.i = (v0->j * v1->k) - (v0->k * v1->j);
    v2.j = (v0->k * v1->i) - (v0->i * v1->k);
    v2.k = (v0->i * v1->j) - (v0->j * v1->i);
    return v2;
}

/* ****************************************************************
 *  double geo_norm (*v0)
 *   DIR_VECT   *v0         (mod) vector to be normalized
 *
 *  Returns the length of the vector before normalization.
 *   Returns zero if the vector could not be normalized.
 *   Note that the input vector is modified.
 */
double geo_norm (v0)
DIR_VECT        *v0;
{   double          len;
    if ((len = geo_dot(v0, v0)) <= 0.0) return FALSE;
    len = sqrt((double)len);
    v0->i /= len;
    v0->j /= len;
    v0->k /= len;
    return TRUE;
}

/* ****************************************************************
 *  double geo_line (*p0, *p1, *v)
 *   POINT3     *p0, *p1    (in)  from, to points
 *   LINE       *v          (out) generated line
 *
 *  Returns the distance from p0 to p1.  Returns 0.0
 *   on error (p0=p1).  The line is expressed in parametric
 *   form with a start point (p0) and a normalized direction
 *   vector. The vector direction of the line is normalized.
 */
double geo_line (p0, p1, v)
POINT3          *p0, *p1;
LINE            *v;
{   v->start = *p0;
    v->dir.i = p1->x - p0->x;
    v->dir.j = p1->y - p0->y;
    v->dir.k = p1->z - p0->z;
    return geo_norm(&v->dir);
}

/* ****************************************************************
 *  DIR_VECT *geo_rfl (V, N)
 *   DIR_VECT   *V          (in)  incident vector
 *   DIR_VECT   *N          (in)  surface normal
 *
 *  Returns the reflected direction vector.  The reflected
 *   direction is computed using the method given by Whitted
 *   (1980), refer to Eq.(\*(rM).
 */
DIR_VECT *geo_rfl (V, N)
DIR_VECT        *V;
DIR_VECT        *N;
{   double              N_dot_V;
    static DIR_VECT     rfl;
    N_dot_V = geo_dot (N,V);
    rfl.i = (2.0 * N_dot_V * N->i) - V->i;
    rfl.j = (2.0 * N_dot_V * N->j) - V->j;
    rfl.k = (2.0 * N_dot_V * N->k) - V->k;
    return &rfl;
}

/* ****************************************************************
 *  DIR_VECT *geo_rfr (V, N, ni, nt)
 *   DIR_VECT   *V          (in)  incident vector
 *   DIR_VECT   *N          (in)  surface normal
 *   double     ni          (in)  index of refraction for the
 *                              material on the front of the
 *                              interface (same side as N)
 *   double     nt          (in)  index of refraction for the
 *                              material on the back of the
 *                              interface (opposite size as N)
 *
 *  Returns the refracted vector, if there complete internal
 *   refracted (no refracted vector) then a NULL vector is
 *   NULL is returned.  The vector is computed using the
 *   method given by Hall (1983) with enhancements as
 *   described in Appendix III.
 */
DIR_VECT *geo_rfr (V, N, ni, nt)
DIR_VECT        *V;
DIR_VECT        *N;
double          ni;
double          nt;
{   static DIR_VECT T;      /* the refracted vector */
    DIR_VECT    sin_T;      /* sin vect of the refracted vect */
    DIR_VECT    cos_V;      /* cos vect of the incident vect */
    double      len_sin_T;  /* length of sin T squared */
    double      n_mult;     /* ni over nt */
    double      N_dot_V;
    double      N_dot_T;
    if ((N_dot_V = geo_dot(N,V)) > 0.0) n_mult = ni / nt;
    else n_mult = nt / ni;
    cos_V.i = N_dot_V * N->i;
    cos_V.j = N_dot_V * N->j;
    cos_V.k = N_dot_V * N->k;
    sin_T.i = (n_mult) * (cos_V.i - V->i);
    sin_T.j = (n_mult) * (cos_V.j - V->j);
    sin_T.k = (n_mult) * (cos_V.k - V->k);
    if ((len_sin_T = geo_dot(&sin_T, &sin_T)) >= 1.0)
	return NULL;    /* internal reflection */
    N_dot_T = sqrt(1.0 - len_sin_T);
    if (N_dot_V < 0.0) N_dot_T = -N_dot_T;
    T.i = sin_T.i - (N->i * N_dot_T);
    T.j = sin_T.j - (N->j * N_dot_T);
    T.k = sin_T.k - (N->k * N_dot_T);
    return &T;
}

/* ****************************************************************
 *  DIR_VECT *geo_H (L, V)
 *   DIR_VECT   *L          (in) incident vector
 *   DIR_VECT   *V          (in) reflected vector
 *
 *  Returns H, NULL on error (if L+H = 0).  Refer
 *   to Eq.(\*(rO).
 */
DIR_VECT *geo_H (L, V)
DIR_VECT        *L;
DIR_VECT        *V;
{   static DIR_VECT         H;
    H.i = L->i + V->i;
    H.j = L->j + V->j;
    H.k = L->k + V->k;
    if (!geo_norm(&H)) return NULL;
    return &H;
}

/* ****************************************************************
 *  DIR_VECT *geo_Ht(L, T, ni, nt)
 *   DIR_VECT   *L          (in) incident vector
 *   DIR_VECT   *T          (in) transmitted hector
 *   double     ni          (in) incident index
 *   double     nt          (in) transmitted index
 *
 *  Returns H' oriented to the same side of the surface
 *   as L computed using the method suggested by
 *   Hall (1983).   Returns NULL on error (if the angle
 *   between V and L is less than the critical angle).
 */
DIR_VECT *geo_Ht (L, T, ni, nt)
DIR_VECT        *L;
DIR_VECT        *T;
double          ni;
double          nt;
{
    float                   L_dot_T;
    float                   divisor;
    static DIR_VECT         Ht;

    L_dot_T = -(geo_dot(L,T));
    /* check for special cases */
    if (ni == nt) {
	if (L_dot_T == 1.0) return L;
	else return NULL;
    }
    if (ni < nt) {
	if (L_dot_T < ni/nt) return NULL;
	divisor = (nt / ni) - 1.0;
	Ht.i = -(((L->i + T->i) / divisor) + T->i);
	Ht.j = -(((L->j + T->j) / divisor) + T->j);
	Ht.k = -(((L->k + T->k) / divisor) + T->k);
    }
    else
    {   if (L_dot_T < nt/ni) return NULL;
	divisor = (ni / nt) - 1.0;
	Ht.i = ((L->i + T->i) / divisor) + L->i;
	Ht.j = ((L->j + T->j) / divisor) + L->j;
	Ht.k = ((L->k + T->k) / divisor) + L->k;
    }
    (void)geo_norm(&Ht);
    return &Ht;
}
/* ************************************************************* */