/*	
	Miscellaneous routines
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#ifndef WIN32
#include <sys/time.h>
#include <sys/resource.h>
#else
#include <windows.h>
#include <mmsystem.h>
#endif
#include "paulslib.h"

/*-------------------------------------------------------------------------
	Compute a point along a spline curve
	The knots u[] should have been calculated before using this
	There are n+1 control points
	t is the degree, typically 3 or 4
	v ranges from 0 to n-t+2, start of the curve to the end
*/
void SplinePoint(int *u,int n,int t,double v,XYZ *control,XYZ *output)
{
  int k;
  double b;

  output->x = 0;
  output->y = 0;
  output->z = 0;

  for (k=0;k<=n;k++) {
    b = SplineBlend(k,t,u,v);
    output->x += control[k].x * b;
    output->y += control[k].y * b;
    output->z += control[k].z * b;
  }
}

/*-------------------------------------------------------------------------
   Calculate the blending value for the spline recursively
	If the numerator and denominator are 0 the result is 0
*/
double SplineBlend(int k,int t,int *u,double v)
{
  double value;

  if (t == 1) {
    if ((u[k] <= v) && (v < u[k+1]))
      value = 1;
    else
      value = 0;
  } else {
    if ((u[k+t-1] == u[k]) && (u[k+t] == u[k+1]))
      value = 0;
    else if (u[k+t-1] == u[k])
      value = (u[k+t] - v) / (u[k+t] - u[k+1]) * SplineBlend(k+1,t-1,u,v);
    else if (u[k+t] == u[k+1])
      value = (v - u[k]) / (u[k+t-1] - u[k]) * SplineBlend(k,t-1,u,v);
    else
      value = (v - u[k]) / (u[k+t-1] - u[k]) * SplineBlend(k,t-1,u,v) +
              (u[k+t] - v) / (u[k+t] - u[k+1]) * SplineBlend(k+1,t-1,u,v);
  }

  return(value);
}

/*-------------------------------------------------------------------------
	Create all the points along a spline curve
	Control points "inp", "n" of them.
   Knots "knots", degree "t".
   Ouput curve "outp", "res" of them.
*/
void SplineCurve(XYZ *inp,int nn,int *knots,int t,XYZ *outp,int res)
{
   int i;
   double interval,increment;

   interval = 0;
   increment = (nn - t + 2) / (double)(res - 1);
   for (i=0;i<res-1;i++) {
      SplinePoint(knots,nn,t,interval,inp,&(outp[i]));
      interval += increment;
   }
	outp[res-1] = inp[nn];
}

/*-------------------------------------------------------------------------
   Figure out the knots, u[]
	There are n + t + 1 knots
*/
void SplineKnots(int *u,int n,int t)
{
  int j;

  for (j=0;j<=n+t;j++) {
    if (j < t)
      u[j] = 0;
    else if (j <= n)
      u[j] = j - t + 1;
    else if (j > n)
      u[j] = n - t + 2;
  }
}

/*-------------------------------------------------------------------------
   Three control point Bezier interpolation
	mu ranges from 0 to 1, start to end of the curve
*/
XYZ Bezier3(XYZ p1,XYZ p2,XYZ p3,double mu)
{
   double mum1,mum12,mu2;
   XYZ p;

   mu2 = mu * mu;
   mum1 = 1 - mu;
   mum12 = mum1 * mum1;
   p.x = p1.x * mum12 + 2 * p2.x * mum1 * mu + p3.x * mu2;
   p.y = p1.y * mum12 + 2 * p2.y * mum1 * mu + p3.y * mu2;
   p.z = p1.z * mum12 + 2 * p2.z * mum1 * mu + p3.z * mu2;

   return(p);
}

/*-------------------------------------------------------------------------
   Four control point Bezier interpolation
	mu ranges from 0 to 1, start to end of curve
*/
XYZ Bezier4(XYZ p1,XYZ p2,XYZ p3,XYZ p4,double mu)
{
   double mum1,mum13,mu3;
   XYZ p;

   mum1 = 1 - mu;
   mum13 = mum1 * mum1 * mum1;
   mu3 = mu * mu * mu;

   p.x = mum13*p1.x + 3*mu*mum1*mum1*p2.x + 3*mu*mu*mum1*p3.x + mu3*p4.x;
   p.y = mum13*p1.y + 3*mu*mum1*mum1*p2.y + 3*mu*mu*mum1*p3.y + mu3*p4.y;
   p.z = mum13*p1.z + 3*mu*mum1*mum1*p2.z + 3*mu*mu*mum1*p3.z + mu3*p4.z;

   return(p);
}

/*-------------------------------------------------------------------------
	General Bezier curve
	Number of control points is n+1
	0 <= mu < 1
*/
XYZ Bezier(XYZ *p,int n,double mu)
{
   int k,kn,nn,nkn;
   double blend,muk,munk;
   XYZ b = {0.0,0.0,0.0};

	muk = 1;
	munk = pow(1-mu,(double)n);

   for (k=0;k<=n;k++) {
		nn = n;
		kn = k;
		nkn = n - k;
		blend = muk * munk;
		muk *= mu;
		munk /= (1-mu);
		while (nn >= 1) {
			blend *= nn;
			nn--;
			if (kn > 1) {
				blend /= (double)kn;
				kn--;	
			}
			if (nkn > 1) {
				blend /= (double)nkn;
				nkn--;	
			}
		}
      b.x += p[k].x * blend;
      b.y += p[k].y * blend;
      b.z += p[k].z * blend;
   }

   return(b);
}

/*
   Calculate blending function for Bezier curves/surfaces
*/
double BezierBlend(int k,double mu,int n)
{
   int nn,kn,nkn;
   double blend=1;

   nn = n;
   kn = k;
   nkn = n - k;

   while (nn >= 1) {
      blend *= nn;
      nn--;
      if (kn > 1) {
         blend /= (double)kn;
         kn--;
      }
      if (nkn > 1) {
         blend /= (double)nkn;
         nkn--;
      }
   }
   if (k > 0)
      blend *= pow(mu,(double)k);
   if (n-k > 0)
      blend *= pow(1-mu,(double)(n-k));

   return(blend);
}

/*
	Piecewise cubic bezier curve as defined by Adobe in Postscript
	The two end points are p0 and p3
	Their associated control points are p1 and p2
*/
XYZ CubicBezier(XYZ p0,XYZ p1,XYZ p2,XYZ p3,double mu)
{
	XYZ a,b,c,p;

	c.x = 3 * (p1.x - p0.x);
   c.y = 3 * (p1.y - p0.y);
   c.z = 3 * (p1.z - p0.z);
	b.x = 3 * (p2.x - p1.x) - c.x;
	b.y = 3 * (p2.y - p1.y) - c.y;
	b.z = 3 * (p2.z - p1.z) - c.z;
	a.x = p3.x - p0.x - c.x - b.x;
   a.y = p3.y - p0.y - c.y - b.y;
   a.z = p3.z - p0.z - c.z - b.z;

	p.x = a.x * mu * mu * mu + b.x * mu * mu + c.x * mu + p0.x;
   p.y = a.y * mu * mu * mu + b.y * mu * mu + c.y * mu + p0.y;
   p.z = a.z * mu * mu * mu + b.z * mu * mu + c.z * mu + p0.z;

	return(p);
}

/*-------------------------------------------------------------------------
   Return the angle between two vectors on a plane
   The angle is from vector 1 to vector 2, positive anticlockwise
   The result is between -pi -> pi
*/
double Angle2D(double x1, double y1, double x2, double y2)
{
   double dtheta,theta1,theta2;

   theta1 = atan2(y1,x1);
   theta2 = atan2(y2,x2);
   dtheta = theta2 - theta1;
   while (dtheta > PI)
      dtheta -= TWOPI;
   while (dtheta < -PI)
      dtheta += TWOPI;

   return(dtheta);
}

/*-------------------------------------------------------------------------
	Dot product of two vectors in 3 space p1 dot p2
*/
double DotProduct(XYZ p1,XYZ p2)
{
   return(p1.x*p2.x + p1.y*p2.y + p1.z*p2.z);
}

/*-------------------------------------------------------------------------
   Return the angle in radians between two vectors (0..pi)
*/
double VectorAngle(XYZ p1,XYZ p2)
{
	double m1,m2;
   double costheta;

	m1 = MODULUS(p1);
   m2 = MODULUS(p2);
	if (m1*m2 <= EPS)
		return(0.0);
	else
		costheta = (p1.x*p2.x + p1.y*p2.y + p1.z*p2.z) / (m1*m2);

	if (costheta <= -1)
		return(PI);
	else if (costheta >= 1)
		return(0.0);
	else
		return(acos(costheta));
}

/*-------------------------------------------------------------------------
   Cross product between two vectors p = p1 x p2
*/
XYZ CrossProduct(XYZ p1,XYZ p2)
{
   XYZ p;

   p.x = p1.y * p2.z - p1.z * p2.y;
   p.y = p1.z * p2.x - p1.x * p2.z;
   p.z = p1.x * p2.y - p1.y * p2.x;

   return(p);
}

/*-------------------------------------------------------------------------
	Return the distance between two points
*/
double VectorLength(XYZ p1,XYZ p2)
{
	XYZ d;
	
	d.x = p1.x - p2.x;
	d.y = p1.y - p2.y;
	d.z = p1.z - p2.z;

	return(sqrt(d.x*d.x + d.y*d.y + d.z*d.z));
}

/*-------------------------------------------------------------------------
   Set thelength of a vector
*/
void SetVectorLength(XYZ *p,double len)
{
	Normalise(p);
	p->x *= len;
   p->y *= len;
   p->z *= len;
}

/*-------------------------------------------------------------------------
	Calculate the length of a vector
*/
double Modulus(XYZ p)
{
    return(sqrt(p.x * p.x + p.y * p.y + p.z * p.z));
}

/*-------------------------------------------------------------------------
	Normalise a vector
*/
void Normalise(XYZ *p)
{
   double length;

   length = p->x * p->x + p->y * p->y + p->z * p->z;
   if (length > 0) {
		length = sqrt(length);
      p->x /= length;
      p->y /= length;
      p->z /= length;
   } else {
		p->x = 0;
		p->y = 0;
		p->z = 0;
	}	
}

/*-------------------------------------------------------------------------
	Multiply a vector by a constant
*/
XYZ VectorMul(XYZ p,double c)
{
	XYZ p1;

	p1.x = c * p.x;
	p1.y = c * p.y;
	p1.z = c * p.z;

	return(p1);
}

/*-------------------------------------------------------------------------
   Create a vector
*/
void MakeVector(XYZ *p,double x,double y,double z)
{
	p->x = x;
	p->y = y;
	p->z = z;
}

/*-------------------------------------------------------------------------
   Subtract two vectors p = p2 - p1
*/
XYZ VectorSub(XYZ p1,XYZ p2)
{
   XYZ p;

   p.x = p2.x - p1.x;
   p.y = p2.y - p1.y;
   p.z = p2.z - p1.z;

   return(p);
}

/*-------------------------------------------------------------------------
   Add two vectors p = p2 + p1
*/
XYZ VectorAdd(XYZ p1,XYZ p2)
{
   XYZ p;

   p.x = p2.x + p1.x;
   p.y = p2.y + p1.y;
   p.z = p2.z + p1.z;

	return(p);
}

/*-------------------------------------------------------------------------
  Return TRUE if two vectors are "equal" else FALSE
*/
int VectorEqual(XYZ v1,XYZ v2)
{
  if (ABS(v1.x - v2.x) > EPSILON)
    return(FALSE);
  if (ABS(v1.y - v2.y) > EPSILON)
    return(FALSE);
  if (ABS(v1.z - v2.z) > EPSILON)
    return(FALSE);
  return(TRUE);
}

/*-------------------------------------------------------------------------
  Invert a vector
*/
XYZ VectorInvert(XYZ v)
{
	XYZ vi;

	vi.x = -v.x;
	vi.y = -v.y;
	vi.z = -v.z;

	return(vi);
}

/*-------------------------------------------------------------------------
   Force a vector v2 to be perpendicular to v1
*/
void VectorPerp(XYZ v1,XYZ *v2)
{
	XYZ vr;
	double len;

	len = Modulus(*v2);
   vr = CrossProduct(v1,*v2);
   *v2 = CrossProduct(vr,v1);
   Normalise(v2);
	v2->x *= len;
   v2->y *= len;
   v2->z *= len;
}

/*-------------------------------------------------------------------------
   Rotate a 2D vector clockwise about the origin by theta (radians)
*/
XY Vector2DRotate(XY p,double theta)
{
	XY q;

	q.x =  p.x * cos(theta) + p.y * sin(theta);
	q.y = -p.x * sin(theta) + p.y * cos(theta);

	return(q);
}

/*-------------------------------------------------------------------------
	Rotate vectors around each axis
	Clockwise looking into the origin from along the positive axis
*/
XYZ RotateX(XYZ p,double theta)
{
	XYZ q;

	q.x = p.x;
	q.y = p.y * cos(theta) + p.z * sin(theta);
	q.z = -p.y * sin(theta) + p.z * cos(theta);
	return(q);
}
XYZ RotateY(XYZ p,double theta)
{
   XYZ q;

   q.x = p.x * cos(theta) - p.z * sin(theta);
   q.y = p.y;
   q.z = p.x * sin(theta) + p.z * cos(theta);
   return(q);
}
XYZ RotateZ(XYZ p,double theta)
{
   XYZ q;

   q.x = p.x * cos(theta) + p.y * sin(theta);
   q.y = -p.x * sin(theta) + p.y * cos(theta);
	q.z = p.z;
   return(q);
}

/*-------------------------------------------------------------------------
	Return the midpoint between two vectors
*/
XYZ MidPoint(XYZ p1,XYZ p2)
{
	XYZ p;

	p.x = (p1.x + p2.x) / 2;
	p.y = (p1.y + p2.y) / 2;
	p.z = (p1.z + p2.z) / 2;

	return(p);
}

/*-------------------------------------------------------------------------
	Return the centroid of n points in 3 space
*/
XYZ Centroid(XYZ *p,int n)
{
	int i;
	XYZ sum;

	sum.x = 0;
	sum.y = 0;
	sum.z = 0;
	for (i=0;i<n;i++) {
		sum.x += p[i].x;
		sum.y += p[i].y;
		sum.z += p[i].z;
	}
	sum.x /= n;
	sum.y /= n;
	sum.z /= n;

	return(sum);
}

/*-------------------------------------------------------------------------
   This function is due to Rob Fowler.  Given dy and dx between 2 points
   A and B, we calculate a number in [0.0, 8.0) which is a monotonic
   function of the direction from A to B. 

   (0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0) correspond to
   (  0,  45,  90, 135, 180, 225, 270, 315, 360) degrees, measured
   counter-clockwise from the positive x axis.
*/
double FowlerAngle(double dy,double dx)
{
    double adx, ady;    /* Absolute Values of Dx and Dy */
    int    code;        /* Angular Region Classification Code */

    adx = (dx < 0) ? -dx : dx;  /* Compute the absolute values. */
    ady = (dy < 0) ? -dy : dy;

    code = (adx < ady) ? 1 : 0;
    if (dx < 0)  code += 2;
    if (dy < 0)  code += 4;

    switch (code) {
    case 0: return (dx==0) ? 0 : ady/adx;  /* [  0, 45] */
    case 1: return (2.0 - (adx/ady));      /* ( 45, 90] */
    case 3: return (2.0 + (adx/ady));      /* ( 90,135) */
    case 2: return (4.0 - (ady/adx));      /* [135,180] */
    case 6: return (4.0 + (ady/adx));      /* (180,225] */
    case 7: return (6.0 - (adx/ady));      /* (225,270) */
    case 5: return (6.0 + (adx/ady));      /* [270,315) */
    case 4: return (8.0 - (ady/adx));      /* [315,360) */
    }
	 return(0.0);
}

/*-------------------------------------------------------------------------
	Calculate the unit normal at p given two other points 
	p1,p2 on the surface. The normal points in the direction 
	of p1 crossproduct p2
*/
XYZ CalcNormal(XYZ p,XYZ p1,XYZ p2)
{
	XYZ n,pa,pb;

   pa.x = p1.x - p.x;
   pa.y = p1.y - p.y;
   pa.z = p1.z - p.z;
	pb.x = p2.x - p.x;
	pb.y = p2.y - p.y;
	pb.z = p2.z - p.z;
   n.x = pa.y * pb.z - pa.z * pb.y;
   n.y = pa.z * pb.x - pa.x * pb.z;
   n.z = pa.x * pb.y - pa.y * pb.x;
	Normalise(&n);

	return(n);
}

/*
   Calculate the line segment PaPb that is the shortest route between
   two lines P1P2 and P3P4. Calculate also the values of mua and mub where
      Pa = P1 + mua (P2 - P1)
      Pb = P3 + mub (P4 - P3)
   Return FALSE if no solution exists.
*/
int LineLineIntersect(
   XYZ p1,XYZ p2,XYZ p3,XYZ p4,XYZ *pa,XYZ *pb,
   double *mua, double *mub)
{
   XYZ p13,p43,p21;
   double d1343,d4321,d1321,d4343,d2121;
   double numer,denom;

   p13.x = p1.x - p3.x;
   p13.y = p1.y - p3.y;
   p13.z = p1.z - p3.z;
   p43.x = p4.x - p3.x;
   p43.y = p4.y - p3.y;
   p43.z = p4.z - p3.z;
   if (ABS(p43.x) < EPS && ABS(p43.y) < EPS && ABS(p43.z) < EPS)
      return(FALSE);
   p21.x = p2.x - p1.x;
   p21.y = p2.y - p1.y;
   p21.z = p2.z - p1.z;
   if (ABS(p21.x) < EPS && ABS(p21.y) < EPS && ABS(p21.z) < EPS)
      return(FALSE);

   d1343 = p13.x * p43.x + p13.y * p43.y + p13.z * p43.z;
   d4321 = p43.x * p21.x + p43.y * p21.y + p43.z * p21.z;
   d1321 = p13.x * p21.x + p13.y * p21.y + p13.z * p21.z;
   d4343 = p43.x * p43.x + p43.y * p43.y + p43.z * p43.z;
   d2121 = p21.x * p21.x + p21.y * p21.y + p21.z * p21.z;

   denom = d2121 * d4343 - d4321 * d4321;
   if (ABS(denom) < EPS)
      return(FALSE);
   numer = d1343 * d4321 - d1321 * d4343;

   *mua = numer / denom;
   *mub = d1343 + d4321 * (*mua) / d4343;

   pa->x = p1.x + *mua * p21.x;
   pa->y = p1.y + *mua * p21.y;
   pa->z = p1.z + *mua * p21.z;
   pb->x = p3.x + *mub * p43.x;
   pb->y = p3.y + *mub * p43.y;
   pb->z = p3.z + *mub * p43.z;

   return(TRUE);
}

/*-------------------------------------------------------------------------
	Find the closest point in a cloud of points
	Return the index to the closest point, the closest point, the distance
*/
int FindClosest(XYZ p,XYZ *poly,int n,double *dmin,XYZ *pmin)
{
	int k,kmin;
	double d;

   *dmin = 1e32;
   kmin = 0;
   for (k=0;k<n;k++) {
      d = VectorLength(poly[k],p);
      if (d < *dmin) {
         *dmin = d;
         kmin = k;
      }
   }

	*pmin = poly[kmin];
	return(kmin);
}

/*-------------------------------------------------------------------------
   Clip a 3 vertex facet in place
   The 3 point facet is defined by vertices p[0],p[1],p[2], "p[3]"
      There must be a fourth point as a 4 point facet may result
   The normal to the plane is n
   A point on the plane is p0
   The side of the plane containing the normal is clipped away
   Return the number of vertices in the clipped polygon
*/
int ClipFacet(XYZ *p,XYZ n,XYZ p0)
{
   double A,B,C,D;
   double l;
   double side[3];
   XYZ q;

   /*
      Determine the equation of the plane as
      Ax + By + Cz + D = 0
   */
   l = sqrt(n.x*n.x + n.y*n.y + n.z*n.z);
   A = n.x / l;
   B = n.y / l;
   C = n.z / l;
   D = -(n.x*p0.x + n.y*p0.y + n.z*p0.z);

   /*
      Evaluate the equation of the plane for each vertex
      If side < 0 then it is on the side to be retained
      else it is to be clippped
   */
   side[0] = A*p[0].x + B*p[0].y + C*p[0].z + D;
   side[1] = A*p[1].x + B*p[1].y + C*p[1].z + D;
   side[2] = A*p[2].x + B*p[2].y + C*p[2].z + D;

   /* Are all the vertices are on the clipped side */
   if (side[0] >= 0 && side[1] >= 0 && side[2] >= 0)
      return(0);

   /* Are all the vertices on the not-clipped side */
   if (side[0] <= 0 && side[1] <= 0 && side[2] <= 0)
      return(3);

   /* Is p0 the only point on the clipped side */
   if (side[0] > 0 && side[1] < 0 && side[2] < 0) {
      q.x = p[0].x - side[0] * (p[2].x - p[0].x) / (side[2] - side[0]);
      q.y = p[0].y - side[0] * (p[2].y - p[0].y) / (side[2] - side[0]);
      q.z = p[0].z - side[0] * (p[2].z - p[0].z) / (side[2] - side[0]);
      p[3] = q;
      q.x = p[0].x - side[0] * (p[1].x - p[0].x) / (side[1] - side[0]);
      q.y = p[0].y - side[0] * (p[1].y - p[0].y) / (side[1] - side[0]);
      q.z = p[0].z - side[0] * (p[1].z - p[0].z) / (side[1] - side[0]);
      p[0] = q;
      return(4);
   }

   /* Is p1 the only point on the clipped side */
   if (side[1] > 0 && side[0] < 0 && side[2] < 0) {
      p[3] = p[2];
      q.x = p[1].x - side[1] * (p[2].x - p[1].x) / (side[2] - side[1]);
      q.y = p[1].y - side[1] * (p[2].y - p[1].y) / (side[2] - side[1]);
      q.z = p[1].z - side[1] * (p[2].z - p[1].z) / (side[2] - side[1]);
      p[2] = q;
      q.x = p[1].x - side[1] * (p[0].x - p[1].x) / (side[0] - side[1]);
      q.y = p[1].y - side[1] * (p[0].y - p[1].y) / (side[0] - side[1]);
      q.z = p[1].z - side[1] * (p[0].z - p[1].z) / (side[0] - side[1]);
      p[1] = q;
      return(4);
   }

   /* Is p2 the only point on the clipped side */
   if (side[2] > 0 && side[0] < 0 && side[1] < 0) {
      q.x = p[2].x - side[2] * (p[0].x - p[2].x) / (side[0] - side[2]);
      q.y = p[2].y - side[2] * (p[0].y - p[2].y) / (side[0] - side[2]);
      q.z = p[2].z - side[2] * (p[0].z - p[2].z) / (side[0] - side[2]);
      p[3] = q;
      q.x = p[2].x - side[2] * (p[1].x - p[2].x) / (side[1] - side[2]);
      q.y = p[2].y - side[2] * (p[1].y - p[2].y) / (side[1] - side[2]);
      q.z = p[2].z - side[2] * (p[1].z - p[2].z) / (side[1] - side[2]);
      p[2] = q;
      return(4);
   }

   /* Is p0 the only point on the not-clipped side */
   if (side[0] < 0 && side[1] > 0 && side[2] > 0) {
      q.x = p[0].x - side[0] * (p[1].x - p[0].x) / (side[1] - side[0]);
      q.y = p[0].y - side[0] * (p[1].y - p[0].y) / (side[1] - side[0]);
      q.z = p[0].z - side[0] * (p[1].z - p[0].z) / (side[1] - side[0]);
      p[1] = q;
      q.x = p[0].x - side[0] * (p[2].x - p[0].x) / (side[2] - side[0]);
      q.y = p[0].y - side[0] * (p[2].y - p[0].y) / (side[2] - side[0]);
      q.z = p[0].z - side[0] * (p[2].z - p[0].z) / (side[2] - side[0]);
      p[2] = q;
      return(3);
   }

   /* Is p1 the only point on the not-clipped side */
   if (side[1] < 0 && side[0] > 0 && side[2] > 0) {
      q.x = p[1].x - side[1] * (p[0].x - p[1].x) / (side[0] - side[1]);
      q.y = p[1].y - side[1] * (p[0].y - p[1].y) / (side[0] - side[1]);
      q.z = p[1].z - side[1] * (p[0].z - p[1].z) / (side[0] - side[1]);
      p[0] = q;
      q.x = p[1].x - side[1] * (p[2].x - p[1].x) / (side[2] - side[1]);
      q.y = p[1].y - side[1] * (p[2].y - p[1].y) / (side[2] - side[1]);
      q.z = p[1].z - side[1] * (p[2].z - p[1].z) / (side[2] - side[1]);
      p[2] = q;
      return(3);
   }

   /* Is p2 the only point on the not-clipped side */
   if (side[2] < 0 && side[0] > 0 && side[1] > 0) {
      q.x = p[2].x - side[2] * (p[1].x - p[2].x) / (side[1] - side[2]);
      q.y = p[2].y - side[2] * (p[1].y - p[2].y) / (side[1] - side[2]);
      q.z = p[2].z - side[2] * (p[1].z - p[2].z) / (side[1] - side[2]);
      p[1] = q;
      q.x = p[2].x - side[2] * (p[0].x - p[2].x) / (side[0] - side[2]);
      q.y = p[2].y - side[2] * (p[0].y - p[2].y) / (side[0] - side[2]);
      q.z = p[2].z - side[2] * (p[0].z - p[2].z) / (side[0] - side[2]);
      p[0] = q;
      return(3);
   }

   /* Shouldn't get here */
   return(-1);
}

/*-------------------------------------------------------------------------
	Clip a line segment p1,p2 to a bounding box (axis aligned)
	Return FALSE if the line segment isn't included in the box,
	The box is defined by its minimum and maximum corners, bmin and bmax
*/
int ClipLine2Box(XYZ *p1,XYZ *p2,XYZ bmin,XYZ bmax)
{
	double mu;

	/* Clip to bmin.x and bmax.x */
	if (p1->x > p2->x)
		SwapXYZ(p1,p2);
	if (p1->x >= bmax.x || p2->x <= bmin.x)
		return(FALSE);
	if (p1->x < bmin.x && p2->x > bmin.x) {
		mu = (bmin.x - p1->x) / (p2->x - p1->x);
		p1->x = bmin.x;
		p1->y = p1->y + mu * (p2->y - p1->y);
		p1->z = p1->z + mu * (p2->z - p1->z);
	}
	if (p1->x < bmax.x && p2->x > bmax.x) {
		mu = (bmax.x - p1->x) / (p2->x - p1->x);
		p2->x = bmax.x;
		p2->y = p1->y + mu * (p2->y - p1->y);
		p2->z = p1->z + mu * (p2->z - p1->z);
	}
	
   /* Clip to bmin.y and bmax.y */
   if (p1->y > p2->y)
      SwapXYZ(p1,p2);
	if (p1->y >= bmax.y || p2->y <= bmin.y)
		return(FALSE);
   if (p1->y < bmin.y && p2->y > bmin.y) {
      mu = (bmin.y - p1->y) / (p2->y - p1->y);
      p1->x = p1->x + mu * (p2->x - p1->x);
      p1->y = bmin.y;
      p1->z = p1->z + mu * (p2->z - p1->z);
   }
   if (p1->y < bmax.y && p2->y > bmax.y) {
      mu = (bmax.y - p1->y) / (p2->y - p1->y);
      p2->x = p1->x + mu * (p2->x - p1->x);
      p2->y = bmax.y;
      p2->z = p1->z + mu * (p2->z - p1->z);
   } 

   /* Clip to bmin.z and bmax.z */
   if (p1->z > p2->z)
      SwapXYZ(p1,p2);
	if (p1->z >= bmax.z || p2->z <= bmin.z)
		return(FALSE);
   if (p1->z < bmin.z && p2->z > bmin.z) {
      mu = (bmin.z - p1->z) / (p2->z - p1->z);
      p1->x = p1->x + mu * (p2->x - p1->x);
      p1->y = p1->y + mu * (p2->y - p1->y);
      p1->z = bmin.z;
   }
   if (p1->z < bmax.z && p2->z > bmax.z) {
      mu = (bmax.z - p1->z) / (p2->z - p1->z);
      p2->x = p1->x + mu * (p2->x - p1->x);
      p2->y = p1->y + mu * (p2->y - p1->y);
      p2->z = bmax.z;
   }

	if ((p1->x < bmin.x && p2->x < bmin.x) ||
       (p1->y < bmin.y && p2->y < bmin.y) ||
		 (p1->z < bmin.z && p2->z < bmin.z) ||
		 (p1->x > bmax.x && p2->x > bmax.x) ||
		 (p1->y > bmax.y && p2->y > bmax.y) ||
		 (p1->z > bmax.z && p2->z > bmax.z))
		return(FALSE);

	return(TRUE);
}

/*-------------------------------------------------------------------------
	Return true if the point p is inside the box with 
	opposite vertices pmin,pmax
*/
int PointInBox(XYZ p,XYZ pmin,XYZ pmax)
{
	if (p.x <= pmin.x || p.x >= pmax.x)
   	return(FALSE);
   if (p.y <= pmin.y || p.y >= pmax.y)
 		return(FALSE);
   if (p.z <= pmin.z || p.z >= pmax.z)
		return(FALSE); 
	return(TRUE);
}

/*-------------------------------------------------------------------------
   Create a contour slice through a 3 vertex facet "p"
   Given the normal of the cutting plane "n" and a point on the plane "p0"
   Return
       0 if the contour plane doesn't cut the facet
       2 if it does cut the facet, the contour line segment is p1->p2
      -1 for an unexpected occurence
   If a vertex touches the contour plane nothing need to be drawn!?
   Note: the following has been written as a "stand alone" piece of
   code that will work but is far from efficient....
*/
int ContourFacet(XYZ *p,XYZ n,XYZ p0,XYZ *p1,XYZ *p2)
{
   double A,B,C,D;
   double l;
   double side[3];

   /*
      Determine the equation of the plane as
      Ax + By + Cz + D = 0
   */
   l = sqrt(n.x*n.x + n.y*n.y + n.z*n.z);
   A = n.x / l;
   B = n.y / l;
   C = n.z / l;
   D = -(n.x*p0.x + n.y*p0.y + n.z*p0.z);

   /*
      Evaluate the equation of the plane for each vertex
      If side < 0 then it is on the side to be retained
      else it is to be clippped
   */
   side[0] = A*p[0].x + B*p[0].y + C*p[0].z + D;
   side[1] = A*p[1].x + B*p[1].y + C*p[1].z + D;
   side[2] = A*p[2].x + B*p[2].y + C*p[2].z + D;

   /* Are all the vertices on one side */
   if (side[0] >= 0 && side[1] >= 0 && side[2] >= 0)
      return(0);
   if (side[0] <= 0 && side[1] <= 0 && side[2] <= 0)
      return(0);

   /* Is p0 the only point on a side by itself */
   if ((SIGN(side[0]) != SIGN(side[1])) && (SIGN(side[0]) != SIGN(side[2]))) {
      p1->x = p[0].x - side[0] * (p[2].x - p[0].x) / (side[2] - side[0]);
      p1->y = p[0].y - side[0] * (p[2].y - p[0].y) / (side[2] - side[0]);
      p1->z = p[0].z - side[0] * (p[2].z - p[0].z) / (side[2] - side[0]);
      p2->x = p[0].x - side[0] * (p[1].x - p[0].x) / (side[1] - side[0]);
      p2->y = p[0].y - side[0] * (p[1].y - p[0].y) / (side[1] - side[0]);
      p2->z = p[0].z - side[0] * (p[1].z - p[0].z) / (side[1] - side[0]);
      return(2);
   }

   /* Is p1 the only point on a side by itself */
   if ((SIGN(side[1]) != SIGN(side[0])) && (SIGN(side[1]) != SIGN(side[2]))) {
      p1->x = p[1].x - side[1] * (p[2].x - p[1].x) / (side[2] - side[1]);
      p1->y = p[1].y - side[1] * (p[2].y - p[1].y) / (side[2] - side[1]);
      p1->z = p[1].z - side[1] * (p[2].z - p[1].z) / (side[2] - side[1]);
      p2->x = p[1].x - side[1] * (p[0].x - p[1].x) / (side[0] - side[1]);
      p2->y = p[1].y - side[1] * (p[0].y - p[1].y) / (side[0] - side[1]);
      p2->z = p[1].z - side[1] * (p[0].z - p[1].z) / (side[0] - side[1]);
      return(2);
   }

   /* Is p2 the only point on a side by itself */
   if ((SIGN(side[2]) != SIGN(side[0])) && (SIGN(side[2]) != SIGN(side[1]))) {
      p1->x = p[2].x - side[2] * (p[0].x - p[2].x) / (side[0] - side[2]);
      p1->y = p[2].y - side[2] * (p[0].y - p[2].y) / (side[0] - side[2]);
      p1->z = p[2].z - side[2] * (p[0].z - p[2].z) / (side[0] - side[2]);
      p2->x = p[2].x - side[2] * (p[1].x - p[2].x) / (side[1] - side[2]);
      p2->y = p[2].y - side[2] * (p[1].y - p[2].y) / (side[1] - side[2]);
      p2->z = p[2].z - side[2] * (p[1].z - p[2].z) / (side[1] - side[2]);
      return(2);
   }

	/* Shouldn't get here */
	return(-1);
}

/*-------------------------------------------------------------------------
   Determine whether or not the line segment p1,p2
   Intersects the 3 vertex facet bounded by pa,pb,pc
   Return true/false and the intersection point p
   The equation of the line is p = p1 + mu (p2 - p1)
   The equation of the plane is a x + b y + c z + d = 0
                                n.x x + n.y y + n.z z + d = 0
*/
int LineFacet(XYZ p1,XYZ p2,XYZ pa,XYZ pb,XYZ pc,XYZ *p)
{
   double d;
   double a1,a2,a3;
   double total,denom,mu;
   XYZ n,pa1,pa2,pa3;

   /* Make sure the arguments are OK */
   if (VectorLength(p1,p2) < EPS)
      return(FALSE);
	if (VectorLength(pa,pb) < EPS || 
		 VectorLength(pb,pc) < EPS || 
		 VectorLength(pc,pa) < EPS)
		return(FALSE);

   /* Calculate the parameters for the plane */
   n.x = (pb.y - pa.y)*(pc.z - pa.z) - (pb.z - pa.z)*(pc.y - pa.y);
   n.y = (pb.z - pa.z)*(pc.x - pa.x) - (pb.x - pa.x)*(pc.z - pa.z);
   n.z = (pb.x - pa.x)*(pc.y - pa.y) - (pb.y - pa.y)*(pc.x - pa.x);
	if (ABS(n.x) < EPS && ABS(n.y) < EPS && ABS(n.z) < EPS)
		return(FALSE);
	Normalise(&n);
   d = - n.x * pa.x - n.y * pa.y - n.z * pa.z;

   /* Calculate the position on the line that intersects the plane */
   denom = n.x * (p2.x - p1.x) + n.y * (p2.y - p1.y) + n.z * (p2.z - p1.z);
   if (ABS(denom) < EPS)         /* Line and plane don't intersect */
      return(FALSE);
   mu = - (d + n.x * p1.x + n.y * p1.y + n.z * p1.z) / denom;
   p->x = p1.x + mu * (p2.x - p1.x);
   p->y = p1.y + mu * (p2.y - p1.y);
   p->z = p1.z + mu * (p2.z - p1.z);
   if (mu < 0 || mu > 1)   /* Intersection not along line segment */
      return(FALSE);

   /* Determine whether or not the intersection point is bounded by pa,pb,pc */
   pa1.x = pa.x - p->x;
   pa1.y = pa.y - p->y;
   pa1.z = pa.z - p->z;
   Normalise(&pa1);
   pa2.x = pb.x - p->x;
   pa2.y = pb.y - p->y;
   pa2.z = pb.z - p->z;
   Normalise(&pa2);
   pa3.x = pc.x - p->x;
   pa3.y = pc.y - p->y;
   pa3.z = pc.z - p->z;
   Normalise(&pa3);
   a1 = pa1.x*pa2.x + pa1.y*pa2.y + pa1.z*pa2.z;
   a2 = pa2.x*pa3.x + pa2.y*pa3.y + pa2.z*pa3.z;
   a3 = pa3.x*pa1.x + pa3.y*pa1.y + pa3.z*pa1.z;
   total = (acos(a1) + acos(a2) + acos(a3)) * RTOD;
   if (ABS(total - 360) > EPS)
      return(FALSE);

   return(TRUE);
}

/*-------------------------------------------------------------------------
	Determine the intersection point of two line segments
	Return FALSE if the lines don't intersect
*/
int LineIntersect(
double x1, double y1, 
double x2, double y2, 
double x3, double y3,
double x4, double y4,
double *x, double *y)
{
	double mua,mub;
	double denom,numera,numerb;

	denom  = (y4-y3) * (x2-x1) - (x4-x3) * (y2-y1);
	numera = (x4-x3) * (y1-y3) - (y4-y3) * (x1-x3);
	numerb = (x2-x1) * (y1-y3) - (y2-y1) * (x1-x3);
	
	/* Are the line coincident? */
	if (ABS(numera) < EPS && ABS(numerb) < EPS && ABS(denom) < EPS) {
		*x = (x1 + x2) / 2;
		*y = (y1 + y2) / 2;
		return(TRUE);
	}

	/* Are the line parallel */
	if (ABS(denom) < EPS) {
		*x = 0;
		*y = 0;
		return(FALSE);
	}

	/* Is the intersection along the the segments */
	mua = numera / denom;
	mub = numerb / denom;
	if (mua < 0 || mua > 1 || mub < 0 || mub > 1) {
		*x = 0;
		*y = 0;
		return(FALSE);
	}
	*x = x1 + mua * (x2 - x1);
	*y = y1 + mua * (y2 - y1);
	return(TRUE);
}

/*-------------------------------------------------------------------------
	Return the distance of a point to a line and the closest point
	The point in question is p3, the line is between p1 and p2.
*/
double PointLine2D(XY p3,XY p1,XY p2,XY *close,double *mu)
{
   double dx,dy,dd;

   dx = p2.x - p1.x;
   dy = p2.y - p1.y;
   dd = dx * dx + dy * dy;

   /* Are the two points p1 and p2 coincident? */
   if (dd < EPS) {
      *mu = 0;
      close->x = p1.x;
      close->y = p1.y;
   } else {
      *mu = ((p3.x - p1.x) * dx + (p3.y - p1.y) * dy) / dd;
      close->x = p1.x + (*mu) * dx;
      close->y = p1.y + (*mu) * dy;
   }
   dx = p3.x - close->x;
   dy = p3.y - close->y;
   return(sqrt(dx * dx + dy * dy));
}
double PointLine3D(XYZ p3,XYZ p1,XYZ p2,XYZ *close,double *mu)
{
   double dx,dy,dz,dd;

   dx = p2.x - p1.x;
   dy = p2.y - p1.y;
	dz = p2.z - p1.z;
   dd = dx * dx + dy * dy + dz * dz;

   /* Are the two points p1 and p2 coincident? */
   if (dd < EPS) {
      *mu = 0;
      close->x = p1.x;
      close->y = p1.y;
		close->z = p1.z;
   } else {
      *mu = ((p3.x - p1.x)*dx + (p3.y - p1.y)*dy + (p3.z - p1.z)*dz) / dd;
      close->x = p1.x + (*mu) * dx;
      close->y = p1.y + (*mu) * dy;
		close->z = p1.z + (*mu) * dz;
   }
   dx = p3.x - close->x;
   dy = p3.y - close->y;
	dz = p3.z - close->z;
   return(sqrt(dx * dx + dy * dy + dz * dz));
}

/*
	Return the clockwise status of a curve, clockwise or anticlockwise
	n vertices making up curve p
          CLOCKWISE == 1
          ANTICLOCKWISE == -1
*/
int ClockWise(XY *p,int n)
{
	int i,j;
	double area=0;

	if (n < 3)
		return(0);

	for (i=0;i<n;i++) {
		j = (i + 1) % n;
		area += (p[i].x*p[j].y - p[j].x*p[i].y);
	}
	if (area < 0)
		return(CLOCKWISE);
	else if (area > 0)
		return(ANTICLOCKWISE);
	else
		return(0);
} 

/*
	Return whether a polygon in 2D is concave or convex
   return 0 for incomputables eg: colinear points
          CLOCKWISE == 1
          ANTICLOCKWISE == -1
   It is assumed that the polygon is simple 
	(does not intersect itself or have holes)
*/
int ConvexPolygon(XY *p,int n)
{
   int i,j,k;
   int flag = 0;
   double z;

   if (n < 3)
      return(0);

   for (i=0;i<n;i++) {
      j = (i + 1) % n;
      k = (i + 2) % n;
      z  = (p[j].x - p[i].x) * (p[k].y - p[j].y);
      z -= (p[j].y - p[i].y) * (p[k].x - p[j].x);
      if (z < 0)
         flag |= 1;
      else if (z > 0)
         flag |= 2;
      if (flag == 3)
	  		return(CONCAVE);
   }
   if (flag != 0)
      return(CONVEX);
   else
      return(0);
}

/*
	Cartesian to polar (degrees) coordinates
*/
void XYZ2Polar(
	double x,double y,double z,
	double *range,double *above,double *about)
{
	double dr;
	
	*range = sqrt(x*x + y*y + z*z);
	if (*range < EPSILON) {
		*above = 0.0;
		*about = 0.0;
	} else {
		*above = asin(z / *range);
		dr = *range * cos(*above);
		if (ABS(dr) < EPSILON) {
			*about = 0.0;
		} else {
			*about = asin(y / dr);
		}
	}
	
	*about *= RTOD;
	*above *= RTOD;
	*about += 90;
}

/*
	Spherical coordinates theta,phi from vector
*/
void InverseSpherical(XYZ p,double *theta,double *phi)
{
   Normalise(&p);
   *theta = PI + atan2(p.y,p.x);
   *phi = atan2(p.z,sqrt(p.x*p.x+p.y*p.y));
}

/*
	Polar (degrees) to cartesian coordinates
*/
void polar2xyz(
	double range,double above,double about,
	double *x,double *y,double *z)
{
	about -= 90;
	*x = range * cos(above*DTOR) * cos(about*DTOR);
	*y = range * cos(above*DTOR) * sin(about*DTOR);
	*z = range * sin(above*DTOR);
}

/*
   Triangulation subroutine
   Takes as input NV vertices in array pxyz
   Returned is a list of ntri triangular faces in the array v
   These triangles are arranged in a consistent clockwise order.
	The triangle array 'v' should be malloced to 3 * nv
	The vertex array pxyz must be big enough to hold 3 more points
	The vertex array must be sorted in increasing x values

	qsort(p,nv,sizeof(XYZ),XYZCompare);
      :
	int XYZCompare(void *v1,void *v2)
	{
   	XYZ *p1,*p2;
   	p1 = v1;
   	p2 = v2;
   	if (p1->x < p2->x)
   	   return(-1);
   	else if (p1->x > p2->x)
   	   return(1);
   	else
   	   return(0);
	}
*/
int Triangulate(int nv,XYZ *pxyz,ITRIANGLE *v,int *ntri)
{
   int *complete = NULL;
   IEDGE *edges = NULL;
   int nedge = 0;
	int trimax,emax = 200;
	int status = 0;

   int inside;
   int i,j,k;
   double xp,yp,x1,y1,x2,y2,x3,y3,xc,yc,r;
   double xmin,xmax,ymin,ymax,xmid,ymid;
   double dx,dy,dmax;

   /* Allocate memory for the completeness list, flag for each triangle */
	trimax = 4 * nv;
   if ((complete = (int *)malloc(trimax*sizeof(int))) == NULL) {
		status = 1;
      goto skip;
   }

   /* Allocate memory for the edge list */
   if ((edges = (IEDGE *)malloc(emax*(int)sizeof(EDGE))) == NULL) {
      status = 2;
      goto skip;
   }

   /*
      Find the maximum and minimum vertex bounds.
      This is to allow calculation of the bounding triangle
   */
   xmin = pxyz[0].x;
   ymin = pxyz[0].y;
   xmax = xmin;
   ymax = ymin;
   for (i=1;i<nv;i++) {
      if (pxyz[i].x < xmin) xmin = pxyz[i].x;
      if (pxyz[i].x > xmax) xmax = pxyz[i].x;
      if (pxyz[i].y < ymin) ymin = pxyz[i].y;
      if (pxyz[i].y > ymax) ymax = pxyz[i].y;
   }
   dx = xmax - xmin;
   dy = ymax - ymin;
   dmax = (dx > dy) ? dx : dy;
   xmid = (xmax + xmin) / 2.0;
   ymid = (ymax + ymin) / 2.0;

   /*
      Set up the supertriangle
      This is a triangle which encompasses all the sample points.
      The supertriangle coordinates are added to the end of the
      vertex list. The supertriangle is the first triangle in
      the triangle list.
   */
   pxyz[nv+0].x = xmid - 20 * dmax;
   pxyz[nv+0].y = ymid - dmax;
   pxyz[nv+0].z = 0.0;
   pxyz[nv+1].x = xmid;
   pxyz[nv+1].y = ymid + 20 * dmax;
   pxyz[nv+1].z = 0.0;
   pxyz[nv+2].x = xmid + 20 * dmax;
   pxyz[nv+2].y = ymid - dmax;
   pxyz[nv+2].z = 0.0;
   v[0].p1 = nv;
   v[0].p2 = nv+1;
   v[0].p3 = nv+2;
   complete[0] = FALSE;
   *ntri = 1;

   /*
      Include each point one at a time into the existing mesh
   */
   for (i=0;i<nv;i++) {

      xp = pxyz[i].x;
      yp = pxyz[i].y;
      nedge = 0;

      /*
         Set up the edge buffer.
         If the point (xp,yp) lies inside the circumcircle then the
         three edges of that triangle are added to the edge buffer
         and that triangle is removed.
      */
      for (j=0;j<(*ntri);j++) {
         if (complete[j]) 
				continue;
         x1 = pxyz[v[j].p1].x;
         y1 = pxyz[v[j].p1].y;
         x2 = pxyz[v[j].p2].x;
         y2 = pxyz[v[j].p2].y;
         x3 = pxyz[v[j].p3].x;
         y3 = pxyz[v[j].p3].y;
         inside = CircumCircle(xp,yp,x1,y1,x2,y2,x3,y3,&xc,&yc,&r);
         // was if (xc + r < xp)
			if (xc < xp && ((xp-xc)*(xp-xc)) > r)
            complete[j] = TRUE;
         if (inside) {
            /* Check that we haven't exceeded the edge list size */
            if (nedge+3 >= emax) {
					emax += 100;
					if ((edges = (IEDGE *)realloc(edges,emax*(int)sizeof(EDGE))) == NULL) {
						status = 3;
						goto skip;
					}
            }
            edges[nedge+0].p1 = v[j].p1;
            edges[nedge+0].p2 = v[j].p2;
            edges[nedge+1].p1 = v[j].p2;
            edges[nedge+1].p2 = v[j].p3;
            edges[nedge+2].p1 = v[j].p3;
            edges[nedge+2].p2 = v[j].p1;
            nedge += 3;
            v[j] = v[(*ntri)-1];
            complete[j] = complete[(*ntri)-1];
            (*ntri)--;
				j--;
         }
		} 

      /*
         Tag multiple edges
         Note: if all triangles are specified anticlockwise then all
               interior edges are opposite pointing in direction.
      */
      for (j=0;j<nedge-1;j++) {
         for (k=j+1;k<nedge;k++) {
            if ((edges[j].p1 == edges[k].p2) && (edges[j].p2 == edges[k].p1)) {
               edges[j].p1 = -1;
               edges[j].p2 = -1;
               edges[k].p1 = -1;
               edges[k].p2 = -1;
            }
            /* Shouldn't need the following, see note above */
            if ((edges[j].p1 == edges[k].p1) && (edges[j].p2 == edges[k].p2)) {
               edges[j].p1 = -1;
               edges[j].p2 = -1;
               edges[k].p1 = -1;
               edges[k].p2 = -1;
            }
         }
      }

      /*
         Form new triangles for the current point
         Skipping over any tagged edges.
         All edges are arranged in clockwise order.
      */
      for (j=0;j<nedge;j++) {
         if (edges[j].p1 < 0 || edges[j].p2 < 0) 
				continue;
         if ((*ntri) >= trimax) {
            status = 4;
				goto skip;
        	}
         v[*ntri].p1 = edges[j].p1;
         v[*ntri].p2 = edges[j].p2;
         v[*ntri].p3 = i;
         complete[*ntri] = FALSE;
			(*ntri)++;
      }
   }

   /*
      Remove triangles with supertriangle vertices
      These are triangles which have a vertex number greater than nv
   */
   for (i=0;i<(*ntri);i++) {
      if (v[i].p1 >= nv || v[i].p2 >= nv || v[i].p3 >= nv) {
         v[i] = v[(*ntri)-1];
         (*ntri)--;
			i--;
      }
   }

skip:
   free(edges);
   free(complete);
	return(status);
}

/*
   Return TRUE if a point (xp,yp) is inside the circumcircle made up
   of the points (x1,y1), (x2,y2), (x3,y3)
   The circumcircle centre is returned in (xc,yc) and the radius r
   NOTE: A point on the edge is inside the circumcircle
*/
int CircumCircle(double xp,double yp,
   double x1,double y1,double x2,double y2,double x3,double y3,
   double *xc,double *yc,double *rsqr)
{
   double m1,m2,mx1,mx2,my1,my2;
   double dx,dy,drsqr;
   double fabsy1y2 = fabs(y1-y2);
   double fabsy2y3 = fabs(y2-y3);

   /* Check for coincident points */
   if (fabsy1y2 < EPSILON && fabsy2y3 < EPSILON)
       return(FALSE);

   if (fabsy1y2 < EPSILON) {
      m2 = - (x3-x2) / (y3-y2);
      mx2 = (x2 + x3) / 2.0;
      my2 = (y2 + y3) / 2.0;
      *xc = (x2 + x1) / 2.0;
      *yc = m2 * (*xc - mx2) + my2;
   } else if (fabsy2y3 < EPSILON) {
      m1 = - (x2-x1) / (y2-y1);
      mx1 = (x1 + x2) / 2.0;
      my1 = (y1 + y2) / 2.0;
      *xc = (x3 + x2) / 2.0;
      *yc = m1 * (*xc - mx1) + my1;
   } else {
      m1 = - (x2-x1) / (y2-y1);
      m2 = - (x3-x2) / (y3-y2);
      mx1 = (x1 + x2) / 2.0;
      mx2 = (x2 + x3) / 2.0;
      my1 = (y1 + y2) / 2.0;
      my2 = (y2 + y3) / 2.0;
      *xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
      if (fabsy1y2 > fabsy2y3) {
         *yc = m1 * (*xc - mx1) + my1;
      }
      else {
         *yc = m2 * (*xc - mx2) + my2;
      }
   }

   dx = x2 - *xc;
   dy = y2 - *yc;
   *rsqr = dx*dx + dy*dy;

   dx = xp - *xc;
   dy = yp - *yc;
   drsqr = dx*dx + dy*dy;

   return((drsqr <= *rsqr) ? TRUE : FALSE);
}

/*
	Rotate a point p by angle theta around an arbitrary normal r
	Return the rotated point.
	Positive angles are anticlockwise looking down the axis 
	towards the origin.
   Assume right hand coordinate system.  
*/
XYZ ArbitraryRotate(XYZ p,double theta,XYZ r)
{
	XYZ q = {0.0,0.0,0.0};
	double costheta,sintheta;

	Normalise(&r);
	costheta = cos(theta);
	sintheta = sin(theta);

	q.x += (costheta + (1 - costheta) * r.x * r.x) * p.x;
	q.x += ((1 - costheta) * r.x * r.y - r.z * sintheta) * p.y;
	q.x += ((1 - costheta) * r.x * r.z + r.y * sintheta) * p.z;

	q.y += ((1 - costheta) * r.x * r.y + r.z * sintheta) * p.x;
	q.y += (costheta + (1 - costheta) * r.y * r.y) * p.y;
	q.y += ((1 - costheta) * r.y * r.z - r.x * sintheta) * p.z;

	q.z += ((1 - costheta) * r.x * r.z - r.y * sintheta) * p.x;
	q.z += ((1 - costheta) * r.y * r.z + r.x * sintheta) * p.y;
	q.z += (costheta + (1 - costheta) * r.z * r.z) * p.z;

	return(q);
}

/*
   Rotate a point p by angle theta around an arbitrary line segment p1-p2
   Return the rotated point.
   Positive angles are anticlockwise looking down the axis
   towards the origin.
   Assume right hand coordinate system.  
*/
XYZ ArbitraryRotate2(XYZ p,double theta,XYZ p1,XYZ p2)
{
   XYZ q = {0.0,0.0,0.0};
   double costheta,sintheta;
	XYZ r;
	
	r.x = p2.x - p1.x;
   r.y = p2.y - p1.y;
   r.z = p2.z - p1.z;
	p.x -= p1.x;
	p.y -= p1.y;
	p.z -= p1.z;
   Normalise(&r);

   costheta = cos(theta);
   sintheta = sin(theta);

   q.x += (costheta + (1 - costheta) * r.x * r.x) * p.x;
   q.x += ((1 - costheta) * r.x * r.y - r.z * sintheta) * p.y;
   q.x += ((1 - costheta) * r.x * r.z + r.y * sintheta) * p.z;

   q.y += ((1 - costheta) * r.x * r.y + r.z * sintheta) * p.x;
   q.y += (costheta + (1 - costheta) * r.y * r.y) * p.y;
   q.y += ((1 - costheta) * r.y * r.z - r.x * sintheta) * p.z;

   q.z += ((1 - costheta) * r.x * r.z - r.y * sintheta) * p.x;
   q.z += ((1 - costheta) * r.y * r.z + r.x * sintheta) * p.y;
   q.z += (costheta + (1 - costheta) * r.z * r.z) * p.z;

	q.x += p1.x;
	q.y += p1.y;
	q.z += p1.z;
   return(q);
}
XYZ ArbitraryRotate3(XYZ p,double theta,XYZ p1,XYZ p2)
{
   XYZ u,q1,q2;
   double d;

   /* Step 1 */
   q1.x = p.x - p1.x;
   q1.y = p.y - p1.y;
   q1.z = p.z - p1.z;

   u.x = p2.x - p1.x;
   u.y = p2.y - p1.y;
   u.z = p2.z - p1.z;
   Normalise(&u);
   d = sqrt(u.y*u.y + u.z*u.z);

   /* Step 2 */
   if (ABS(d) > 0.000001) {
      q2.x = q1.x;
      q2.y = q1.y * u.z / d - q1.z * u.y / d;
      q2.z = q1.y * u.y / d + q1.z * u.z / d;
   } else {
      q2 = q1;
   }

   /* Step 3 */
   q1.x = q2.x * d - q2.z * u.x;
   q1.y = q2.y;
   q1.z = q2.x * u.x + q2.z * d;

   /* Step 4 */
   q2.x = q1.x * cos(theta) - q1.y * sin(theta);
   q2.y = q1.x * sin(theta) + q1.y * cos(theta);
   q2.z = q1.z;

   /* Inverse of step 3 */
   q1.x =   q2.x * d + q2.z * u.x;
   q1.y =   q2.y;
   q1.z = - q2.x * u.x + q2.z * d;

   /* Inverse of step 2 */
   if (ABS(d) > 0.000001) {
      q2.x =   q1.x;
      q2.y =   q1.y * u.z / d + q1.z * u.y / d;
      q2.z = - q1.y * u.y / d + q1.z * u.z / d;
   } else {
      q2 = q1;
   }

   /* Inverse of step 1 */
   q1.x = q2.x + p1.x;
   q1.y = q2.y + p1.y;
   q1.z = q2.z + p1.z;
   return(q1);
}

double LinearInterpolate(double y1,double y2,double mu)
{
   return(y1 * (1 - mu) + y2 * mu);
}

double CosineInterpolate(double y1,double y2,double mu)
{
   double mu2;

   mu2 = (1 - cos(mu*PI)) / 2;
   return(y1 * (1 - mu2) + y2 * mu2);
}

double CubicInterpolate(double y0,double y1,double y2,double y3,double mu)
{
   double a0,a1,a2,a3;

   a0 = y3 - y2 - y0 + y1;
   a1 = y0 - y1 - a0;
   a2 = y2 - y0;
   a3 = y1;

   return(a0*mu*mu*mu + a1*mu*mu + a2*mu + a3);
}

/* 
	Tension: 1 is high, 0 normal, -1 is low
	Bias: 0 is even, positive is towards first segment, 
   negative towards the other 
*/
double HermiteInterpolate(double y0,double y1,double y2,double y3,
   double mu,double tension,double bias)
{
   double m0,m1,mu2,mu3;
   double a0,a1,a2,a3;

	mu2 = mu * mu;
	mu3 = mu2 * mu;
   m0  = (y1 - y0) * (1 + bias) * (1 - tension) / 2;
   m0 += (y2 - y1) * (1 - bias) * (1 - tension) / 2;
   m1  = (y2 - y1) * (1 + bias) * (1 - tension) / 2;
   m1 += (y3 - y2) * (1 - bias) * (1 - tension) / 2;
   a0 =  2*mu3 - 3*mu2 + 1;
   a1 =    mu3 - 2*mu2 + mu;
   a2 =    mu3 -   mu2;
   a3 = -2*mu3 + 3*mu2;

   return(a0*y1 + a1*m0 + a2*m1 + a3*y2);
}

/*
	Derivation from CONREC
   d               ! matrix of data to contour
   ilb,iub,jlb,jub ! index bounds of data matrix
   x               ! data matrix column coordinates
   y               ! data matrix row coordinates
   nc              ! number of contour levels
   z               ! contour levels in increasing order
*/
void Contour(double **d,int ilb,int iub,int jlb,int jub,
	double *x,double *y,int nc,double *z,
   void (*ConrecLine)(double,double,double,double,double))
{
#define xsect(p1,p2) (h[p2]*xh[p1]-h[p1]*xh[p2])/(h[p2]-h[p1])
#define ysect(p1,p2) (h[p2]*yh[p1]-h[p1]*yh[p2])/(h[p2]-h[p1])

   int m1,m2,m3,case_value;
   double dmin,dmax,x1=0,x2=0,y1=0,y2=0;
   int i,j,k,m;
   double h[5];
   int sh[5];
   double xh[5],yh[5];
   int im[4] = {0,1,1,0},jm[4]={0,0,1,1};
   int castab[3][3][3] = {
     { {0,0,8},{0,2,5},{7,6,9} },
     { {0,3,4},{1,3,1},{4,3,0} },
     { {9,6,7},{5,2,0},{8,0,0} }
   };
	double temp1,temp2;

   for (j=(jub-1);j>=jlb;j--) {
      for (i=ilb;i<=iub-1;i++) {
         temp1 = MIN(d[i][j],d[i][j+1]);
         temp2 = MIN(d[i+1][j],d[i+1][j+1]);
         dmin  = MIN(temp1,temp2);
         temp1 = MAX(d[i][j],d[i][j+1]);
         temp2 = MAX(d[i+1][j],d[i+1][j+1]);
         dmax  = MAX(temp1,temp2);
         if (dmax < z[0] || dmin > z[nc-1]) 
				continue;
         for (k=0;k<nc;k++) {
            if (z[k] < dmin || z[k] > dmax) 
					continue;
            for (m=4;m>=0;m--) {
               if (m > 0) {
                  h[m]  = d[i+im[m-1]][j+jm[m-1]]-z[k];
                  xh[m] = x[i+im[m-1]];
                  yh[m] = y[j+jm[m-1]];
               } else {
                  h[0]  = 0.25 * (h[1]+h[2]+h[3]+h[4]);
                  xh[0] = 0.50 * (x[i]+x[i+1]);
                  yh[0] = 0.50 * (y[j]+y[j+1]);
               }
               if (h[m] > 0.0) 
                  sh[m] = 1;
               else if (h[m] < 0.0) 
                  sh[m] = -1;
               else
                  sh[m] = 0;
            }

			   /*
               Note: at this stage the relative heights of the corners and the
               centre are in the h array, and the corresponding coordinates are
               in the xh and yh arrays. The centre of the box is indexed by 0
               and the 4 corners by 1 to 4 as shown below.
               Each triangle is then indexed by the parameter m, and the 3
               vertices of each triangle are indexed by parameters m1,m2,and m3.
               It is assumed that the centre of the box is always vertex 2
               though this isimportant only when all 3 vertices lie exactly on
               the same contour level, in which case only the side of the box
               is drawn.
                  vertex 4 +-------------------+ vertex 3
                           | \               / |
                           |   \    m-3    /   |
                           |     \       /     |
                           |       \   /       |
                           |  m=2    X   m=2   |       the centre is vertex 0
                           |       /   \       |
                           |     /       \     |
                           |   /    m=1    \   |
                           | /               \ |
                  vertex 1 +-------------------+ vertex 2
            */
            
		   	/* Scan each triangle in the box */
            for (m=1;m<=4;m++) {
               m1 = m;
               m2 = 0;
               if (m != 4)
                  m3 = m + 1;
               else
                  m3 = 1;
               if ((case_value = castab[sh[m1]+1][sh[m2]+1][sh[m3]+1]) == 0)
						continue;
               switch (case_value) {
               case 1: /* Line between vertices 1 and 2 */
                   x1 = xh[m1];
                   y1 = yh[m1];
                   x2 = xh[m2];
                   y2 = yh[m2];
                   break;
               case 2: /* Line between vertices 2 and 3 */
                   x1 = xh[m2];
                   y1 = yh[m2];
                   x2 = xh[m3];
                   y2 = yh[m3];
                   break;
               case 3: /* Line between vertices 3 and 1 */
                   x1 = xh[m3];
                   y1 = yh[m3];
                   x2 = xh[m1];
                   y2 = yh[m1];
                   break;
               case 4: /* Line between vertex 1 and side 2-3 */
                   x1 = xh[m1];
                   y1 = yh[m1];
                   x2 = xsect(m2,m3);
                   y2 = ysect(m2,m3);
                   break;
               case 5: /* Line between vertex 2 and side 3-1 */
                   x1 = xh[m2];
                   y1 = yh[m2];
                   x2 = xsect(m3,m1);
                   y2 = ysect(m3,m1);
                   break;
               case 6: /* Line between vertex 3 and side 1-2 */
                   x1 = xh[m1];
                   y1 = yh[m1];
                   x2 = xsect(m1,m2);
                   y2 = ysect(m1,m2);
                   break;
               case 7: /* Line between sides 1-2 and 2-3 */
                   x1 = xsect(m1,m2);
                   y1 = ysect(m1,m2);
                   x2 = xsect(m2,m3);
                   y2 = ysect(m2,m3);
                   break;
               case 8: /* Line between sides 2-3 and 3-1 */
                   x1 = xsect(m2,m3);
                   y1 = ysect(m2,m3);
                   x2 = xsect(m3,m1);
                   y2 = ysect(m3,m1);
                   break;
               case 9: /* Line between sides 3-1 and 1-2 */
                   x1 = xsect(m3,m1);
                   y1 = ysect(m3,m1);
                   x2 = xsect(m1,m2);
                   y2 = ysect(m1,m2);
                   break;
               default:
                   break;
               }

	  			   /* Finally draw the line */
			      ConrecLine(x1,y1,x2,y2,z[k]); 
            } /* m */
			} /* k - contour */
		} /* i */
	} /* j */
}

/*
   Calculate the intersection of a ray and a sphere
	The line segment is defined from p1 to p2
	The sphere is of radius r and centered at sc
	There are potentially two points of intersection given by
	p = p1 + mu1 (p2 - p1)
	p = p1 + mu2 (p2 - p1)
	Return FALSE if the ray doesn't intersect the sphere.
*/
int RaySphere(XYZ p1,XYZ p2,XYZ sc,double r,double *mu1,double *mu2)
{
   double a,b,c;
   double bb4ac;
   XYZ dp;

   dp.x = p2.x - p1.x;
   dp.y = p2.y - p1.y;
   dp.z = p2.z - p1.z;
   a = dp.x * dp.x + dp.y * dp.y + dp.z * dp.z;
   b = 2 * (dp.x * (p1.x - sc.x) + dp.y * (p1.y - sc.y) + dp.z * (p1.z - sc.z));
   c = sc.x * sc.x + sc.y * sc.y + sc.z * sc.z;
   c += p1.x * p1.x + p1.y * p1.y + p1.z * p1.z;
   c -= 2 * (sc.x * p1.x + sc.y * p1.y + sc.z * p1.z);
   c -= r * r;
   bb4ac = b * b - 4 * a * c;
   if (ABS(a) < EPS || bb4ac < 0) {
		*mu1 = 0;
		*mu2 = 0;
      return(FALSE);
	}

   *mu1 = (-b + sqrt(bb4ac)) / (2 * a);
   *mu2 = (-b - sqrt(bb4ac)) / (2 * a);

   return(TRUE);
}

/*
   Determine the intersection of a line with a plane
   Return false if the line doesn't intersect
*/
int LinePlane(XYZ p1,XYZ p2,PLANE plane,double *mu,XYZ *p)
{
   double numer,denom;
   XYZ dp;

   dp.x = p1.x - p2.x;
   dp.y = p1.y - p2.y;
   dp.z = p1.z - p2.z;
   denom = plane.a * dp.x + plane.b * dp.y + plane.c * dp.z;
   if (ABS(denom) < EPS)
      return(FALSE);
   numer = plane.a * p1.x + plane.b * p1.y + plane.c * p1.z + plane.d;

   *mu = numer / denom;
   p->x = p1.x + (*mu) * (p2.x - p1.x);
   p->y = p1.y + (*mu) * (p2.y - p1.y);
   p->z = p1.z + (*mu) * (p2.z - p1.z);

   return(TRUE);
}
/*
	Same as the above except for a line segment
	Return true if the line segent intersects the plane, return intersection point
*/
int LineSegmentPlane(XYZ p1,XYZ p2,PLANE plane,XYZ *p)
{
   double numer,denom;
   XYZ dp;
  	double mu;

   dp.x = p1.x - p2.x;
   dp.y = p1.y - p2.y;
   dp.z = p1.z - p2.z;
   denom = plane.a * dp.x + plane.b * dp.y + plane.c * dp.z;
   if (ABS(denom) < EPS)
      return(FALSE);
   numer = plane.a * p1.x + plane.b * p1.y + plane.c * p1.z + plane.d;

   mu = numer / denom;
	if (mu < 0 || mu > 1)
		return(FALSE);
   p->x = p1.x + mu * (p2.x - p1.x);
   p->y = p1.y + mu * (p2.y - p1.y); 
   p->z = p1.z + mu * (p2.z - p1.z); 

   return(TRUE);
}

/* FILE IO ---------------------------------------------------------*/

/*
   Encode an RLE buffer
   Return the length in the output buffer
*/
int RLECompress(unsigned char *output,unsigned char *input,int length)
{
   int count = 0,index,i;
   unsigned char pixel;
   int out = 0;

   while (count < length) {
      index = count;
      pixel = input[index++];
      while (index < length && index - count < 127 && input[index] == pixel)
         index++;
      if (index - count == 1) {
         /*
            Failed to "replicate" the current pixel. See how many to copy.
            Avoid a replicate run of only 2-pixels after a literal run. There
            is no gain in this, and there is a risK of loss if the run after
            the two identical pixels is another literal run. So search for
            3 identical pixels.
         */
         while (index < length && 
					 index - count < 127 &&
                (input[index] != input[index-1] || 
					  (index > 1 && 
					  input[index] != input[index-2])))
            index++;
         /*
            Check why this run stopped. If it found two identical pixels, reset
            the index so we can add a run. Do this twice: the previous run
            tried to detect a replicate run of at least 3 pixels. So we may be
            able to back up two pixels if such a replicate run was found.
         */
         while (index < length && input[index] == input[index-1])
            index--;
         output[out++] = (unsigned char)(count - index);
         for (i=count;i<index;i++)
            output[out++] = input[i];
      } else {
         output[out++] = (unsigned char)(index - count);
         output[out++] = pixel;
      } /* if */
      count=index;
   } /* while */
   return(out);
}

/*
   Uncompress a RLE buffer
*/
void RLEUncompress(unsigned char *output,unsigned char *input,int length)
{
   signed char count;

   while (length > 0) {
      count = (signed char)*input++;
      if (count > 0) {                      /* replicate run */
         memset(output,*input++,count);
      } else if (count < 0) {               /* literal run */
         count = (signed char)-count;
         memcpy(output,input,count);
         input += count;
      } /* if */
      output += count;
      length -= count;
   } /* while */
}

/*
	Read a line
*/
int ReadLine(FILE *fptr,char *s,int lmax)
{
	int i=0,c;

	s[0] = '\0';
	while ((c = fgetc(fptr)) != '\n' && c != '\r') {
		if (c == EOF)
			return(FALSE);
		s[i] = c;
		i++;
		s[i] = '\0';
		if (i >= lmax)
			break;
	}
	return(TRUE);
}

/*
  Read until a particular character is encountered
*/
int ReadUntil(FILE *fptr,int find,char *s)
{
  int c,i = 0;

  i = 0;
  s[i] = 0;
  while ((c = fgetc(fptr)) != find && c != EOF) {
    s[i] = c;
    i++;
    s[i] = '\0';
  }

  if (c == EOF)
    return(FALSE);
  return(TRUE);
}

/*
	Trim any spaces off the end of a string
*/
void RightTrim(char *s)
{
   int i;

   i = strlen(s) - 1;
   while (i >= 0 && s[i] == ' ')
      s[i--] = '\0';
}

/*
  Read until a particular string is encountered
*/
void ReadToString(FILE *fptr,char *s)
{
  char ss[100];

  while (fscanf(fptr,"%s",ss) == 1) {
    if (strcmp(ss,s) == 0)
      return;
  }
}

/*
   Read a string of n characters
*/
int ReadString(FILE *fptr,char *s,int n)
{
   int i,c;

   s[0] = '\0';
   for (i=0;i<n;i++) {
      c = fgetc(fptr);
      if (c == EOF)
         return(FALSE);
      s[i] = c;
      s[i+1] = '\0';
   }
   return(TRUE);
}

/*
   Read a possibly byte swapped short integer
*/
int ReadShort(FILE *fptr,short int *n,int swap)
{
   unsigned char *cptr,tmp;

   if (fread(n,2,1,fptr) != 1)
      return(FALSE);
	if (swap) {
   	cptr = (unsigned char *)n;
   	tmp = cptr[0];
   	cptr[0] = cptr[1];
   	cptr[1] = tmp;
	}
   return(TRUE);
}

/*
   Write a possibly byte swapped short integer
*/
int WriteShort(FILE *fptr,short n,int swap)
{
   unsigned char *cptr,tmp;

	if (!swap) {
   	if (fwrite(&n,2,1,fptr) != 1)
      	return(FALSE);
   } else {
      cptr = (unsigned char *)(&n);
      tmp = cptr[0];
      cptr[0] = cptr[1];
      cptr[1] = tmp;
      if (fwrite(&n,2,1,fptr) != 1)
         return(FALSE);
   }
   return(TRUE);
}

/*
   Read a possibly byte swapped unsigned short integer
*/
int ReadUShort(FILE *fptr,short unsigned *n,int swap)
{
   unsigned char *cptr,tmp;

   if (fread(n,2,1,fptr) != 1)
      return(FALSE);
	if (swap) {
   	cptr = (unsigned char *)n;
   	tmp = cptr[0];
   	cptr[0] = cptr[1];
   	cptr[1] =tmp;
	}
   return(TRUE);
}

/*
   Write a possibly byte swapped unsigned short integer
*/
int WriteUShort(FILE *fptr,short unsigned n,int swap)
{
   unsigned char *cptr,tmp;

	if (!swap) {
   	if (fwrite(&n,2,1,fptr) != 1)
      	return(FALSE);
   } else {
      cptr = (unsigned char *)(&n);
      tmp = cptr[0];
      cptr[0] = cptr[1];
      cptr[1] =tmp;
      if (fwrite(&n,2,1,fptr) != 1)
         return(FALSE);
   }
   return(TRUE);
}

/*
   Read a possibly byte swapped integer
*/
int ReadInt(FILE *fptr,int *n,int swap)
{
   unsigned char *cptr,tmp;

   if (fread(n,4,1,fptr) != 1)
      return(FALSE);
   if (swap) {
      cptr = (unsigned char *)n;
      tmp = cptr[0];
      cptr[0] = cptr[3];
      cptr[3] = tmp;
      tmp = cptr[1];
      cptr[1] = cptr[2];
      cptr[2] = tmp;
   }
   return(TRUE);
}

/*
   Read a possibly byte swapped unsigned integer
*/
int ReadUInt(FILE *fptr,unsigned int *n,int swap)
{
   unsigned char *cptr,tmp;

   if (fread(n,4,1,fptr) != 1)
      return(FALSE);
   if (swap) {
      cptr = (unsigned char *)n;
      tmp = cptr[0];
      cptr[0] = cptr[3];
      cptr[3] = tmp;
      tmp = cptr[1];
      cptr[1] = cptr[2];
      cptr[2] = tmp;
   }
   return(TRUE);
}

/*
   Write a possibly byte swapped integer
*/
int WriteInt(FILE *fptr,int n,int swap)
{
   unsigned char *cptr,tmp;

	if (!swap) {
   	if (fwrite(&n,4,1,fptr) != 1)
      	return(FALSE);
   } else {
      cptr = (unsigned char *)(&n);
      tmp = cptr[0];
      cptr[0] = cptr[3];
      cptr[3] = tmp;
      tmp = cptr[1];
      cptr[1] = cptr[2];
      cptr[2] = tmp;
      if (fwrite(&n,4,1,fptr) != 1)
         return(FALSE);
   }
   return(TRUE);
}

/*
   Write a possibly byte swapped unsigned integer
*/
int WriteUInt(FILE *fptr,unsigned int n,int swap)
{
   unsigned char *cptr,tmp;

   if (!swap) {
      if (fwrite(&n,4,1,fptr) != 1)
         return(FALSE);
   } else {
      cptr = (unsigned char *)(&n);
      tmp = cptr[0];
      cptr[0] = cptr[3];
      cptr[3] = tmp;
      tmp = cptr[1];
      cptr[1] = cptr[2];
      cptr[2] = tmp;
      if (fwrite(&n,4,1,fptr) != 1)
         return(FALSE);
   }
   return(TRUE);
}

/*
   Read a possibly byte swapped long integer
*/
int ReadLong(FILE *fptr,long *n,int swap)
{
   unsigned char *cptr,tmp;

   if (fread(n,4,1,fptr) != 1)
      return(FALSE);
	if (swap) {
   	cptr = (unsigned char *)n;
   	tmp = cptr[0];
   	cptr[0] = cptr[3];
   	cptr[3] = tmp;
   	tmp = cptr[1];
   	cptr[1] = cptr[2];
   	cptr[2] = tmp;
	}
   return(TRUE);
}

/*
   Write a possibly byte swapped long integer
	arning: machine dependent, this assumes longs are 4 bytes, might be 8!
*/
int WriteLong(FILE *fptr,long n,int swap)
{
   unsigned char *cptr,tmp;

	if (!swap) {
   	if (fwrite(&n,4,1,fptr) != 1)
      	return(FALSE);
   } else {
      cptr = (unsigned char *)(&n);
      tmp = cptr[0];
      cptr[0] = cptr[3];
      cptr[3] = tmp;
      tmp = cptr[1];
      cptr[1] = cptr[2];
      cptr[2] = tmp;
      if (fwrite(&n,4,1,fptr) != 1)
         return(FALSE);
   }
   return(TRUE);
}

/*
   Write a possibly byte swapped unsigned long integer
   Warning: machine dependent, this assumes longs are 4 bytes, might be 8!
*/
int WriteULong(FILE *fptr,unsigned long n,int swap)
{
   unsigned char *cptr,tmp;

   if (!swap) {
      if (fwrite(&n,4,1,fptr) != 1)
         return(FALSE);
   } else {
      cptr = (unsigned char *)(&n);
      tmp = cptr[0];
      cptr[0] = cptr[3];
      cptr[3] = tmp;
      tmp = cptr[1];
      cptr[1] = cptr[2];
      cptr[2] = tmp;
      if (fwrite(&n,4,1,fptr) != 1)
         return(FALSE);
   }
   return(TRUE);
}

/*
   Write a possibly byte swapped double precision number
   Assume IEEE
*/
int WriteDouble(FILE *fptr,double n,int swap)
{
   unsigned char *cptr,tmp;

	if (!swap) {
   	if (fwrite(&n,8,1,fptr) != 1)
      	return(FALSE);
	} else {
   	cptr = (unsigned char *)(&n);
   	tmp = cptr[0];
   	cptr[0] = cptr[7];
   	cptr[7] =tmp;
   	tmp = cptr[1];
   	cptr[1] = cptr[6];
   	cptr[6] = tmp;
   	tmp = cptr[2];
   	cptr[2] = cptr[5];
   	cptr[5] =tmp;
   	tmp = cptr[3];
   	cptr[3] = cptr[4];
   	cptr[4] = tmp;
      if (fwrite(&n,8,1,fptr) != 1)
         return(FALSE);
	}
   return(TRUE);
}

/*
   Read a possibly byte swapped double precision number
   Assume IEEE
*/
int ReadDouble(FILE *fptr,double *n,int swap)
{
   unsigned char *cptr,tmp;

   if (fread(n,8,1,fptr) != 1)
      return(FALSE);
   if (swap) {
      cptr = (unsigned char *)n;
      tmp = cptr[0];
      cptr[0] = cptr[7];
      cptr[7] =tmp;
      tmp = cptr[1];
      cptr[1] = cptr[6];
      cptr[6] = tmp;
      tmp = cptr[2];
      cptr[2] = cptr[5];
      cptr[5] =tmp;
      tmp = cptr[3];
      cptr[3] = cptr[4];
      cptr[4] = tmp;
   }
   return(TRUE);
}

/*
   Read a possibly byte swapped floating point number
   Assume IEEE format
*/
int ReadFloat(FILE *fptr,float *n,int swap)
{
   unsigned char *cptr,tmp;

   if (fread(n,4,1,fptr) != 1)
      return(FALSE);
	if (swap) {
  	 	cptr = (unsigned char *)n;
   	tmp = cptr[0];
   	cptr[0] = cptr[3];
   	cptr[3] =tmp;
   	tmp = cptr[1];
   	cptr[1] = cptr[2];
   	cptr[2] = tmp;
	}
   return(TRUE);
}

/*
   Write a possibly byte swapped float
   Warning: machine dependent, this assumes floats are 4 bytes, might be 8!
*/
int WriteFloat(FILE *fptr,float n,int swap)
{
   unsigned char *cptr,tmp;

   if (!swap) {
      if (fwrite(&n,4,1,fptr) != 1)
         return(FALSE);
   } else {
      cptr = (unsigned char *)(&n);
      tmp = cptr[0];
      cptr[0] = cptr[3];
      cptr[3] = tmp;
      tmp = cptr[1];
      cptr[1] = cptr[2];
      cptr[2] = tmp;
      if (fwrite(&n,4,1,fptr) != 1)
         return(FALSE);
   }
   return(TRUE);
}

/*
   Read past n bytes
*/
int ReadSkip(FILE *fptr,int n)
{
   int i,c;

   for (i=0;i<n;i++) {
      c = fgetc(fptr);
      if (c == EOF)
         return(FALSE);
   }
   return(TRUE);
}

/*
   Read and ignore from a file until some character is reached
*/
void SkipUntil(FILE *fptr,int c)
{
  int cc;

  while ((cc = fgetc(fptr)) != c && cc != EOF)
    ;
}

/*
	Convert a string to lower case
*/
void StrToLower(char *s)
{
  unsigned int i;

  for (i=0;i<strlen(s);i++) {
    if (s[i] >= 'A' && s[i] <= 'Z')
      s[i] += 'a' - 'A';
  }
}

/*
   Convert a string to upper case
*/
void StrToUpper(char *s)
{
  unsigned int i;

  for (i=0;i<strlen(s);i++) {
    if (s[i] >= 'a' && s[i] <= 'z')
      s[i] += 'A' - 'a';
  }
}

/*
	Echo a file, if it exists, to another file
*/
void EchoFile(char *fname,FILE *fp)
{
  int c;
  FILE *fptr;

  if ((fptr = fopen(fname,"rb")) != NULL) {
     while ((c = fgetc(fptr)) != EOF)
        fputc(c,fp);
     fclose(fptr);
   }
}

/*
	Copy parts of a string into another
*/
int IndexCopy(char *sin,char *sout,unsigned int start,unsigned int stop)
{
   unsigned int i,j=0;

   sout[0] = '\0';
   if (strlen(sin) < stop)
      return(FALSE);
   for (i=start;i<=stop;i++) {
      sout[j] = sin[i];
      j++;
      sout[j] = '\0';
   }
   return(TRUE);
}

/* MISC FUNCTIONS ---------------------------------------------------*/

/*
   Clip a double to bounds returning true if necessary
*/
int ClipDouble(double *v,double upper,double lower)
{
	if (*v < lower) {
		*v = lower;
		return(TRUE);		
	}
	if (*v > upper) {
		*v = upper;
		return(TRUE);
	}
	return(FALSE);
}

/*
   Clip a int to bounds returning true if clipping was performed
*/
int ClipInteger(int *v,int lower,int upper)
{
   if (*v < lower) {
      *v = lower;   
      return(TRUE); 
   }
   if (*v > upper) {
      *v = upper;
      return(TRUE);
   }
   return(FALSE);
}

/*
   Swap two doubles
*/
void SwapDouble(double *d1,double *d2)
{
   double tmp;

   tmp = *d1;
   *d1 = *d2;
   *d2 = tmp;
}

/*
	Swap two XYZ points
*/
void SwapXYZ(XYZ *p1,XYZ *p2)
{
	XYZ tmp;

	tmp = *p1;
	*p1 = *p2;
	*p2 = tmp;
}

/*
	Calculate the bounding box of a vector
*/
void MinMaxXYZ(XYZ p,XYZ *min,XYZ *max)
{
   min->x = MIN(min->x,p.x);
   min->y = MIN(min->y,p.y);
   min->z = MIN(min->z,p.z);
   max->x = MAX(max->x,p.x);
   max->y = MAX(max->y,p.y);
   max->z = MAX(max->z,p.z);
}

/*
	Clean a string
*/
void StringClean(char *s)
{
	unsigned int i;
	int c;

	/* Strip non printing charagters from the end */
	while (strlen(s) > 0 && ((c = s[strlen(s)-1]) <= ' ' || c > '~'))
		s[strlen(s)-1] = '\0';

	/* Strip non printing characters from start */
	while (strlen(s) > 0 && (s[0] <= ' ' || s[0] > '~')) {
		for (i=1;i<=strlen(s);i++) 
			s[i-1] = s[i];
	}
}

/* COLOUR STUFF -----------------------------------------------------*/

/*
	Return the colour between c1 and c2 linearly given mu (0..1)
*/
COLOUR RampColour2(double mu,COLOUR c1,COLOUR c2)
{
	COLOUR c;

	if (mu < 0)
		return(c1);
	if (mu > 1)
		return(c2);

	c.r = c1.r + mu * (c2.r - c1.r);
   c.g = c1.g + mu * (c2.g - c1.g);
   c.b = c1.b + mu * (c2.b - c1.b);

	return(c);
}

/*
   Return a colour from one of a number of colour ramps
   type == 1  blue -> cyan -> green -> magenta -> red 
           2  blue -> red 
           3  grey scale 
           4  red -> yellow -> green -> cyan -> blue -> magenta -> red 
           5  green -> yellow 
			  6  green -> magenta
			  7  blue -> green -> red -> green -> blue
			  8  black -> white -> black
			  9  red -> blue -> cyan -> magenta
          10  blue -> cyan -> green -> yellow -> red -> white
			 11  dark brown -> lighter brown (Mars colours, 2 colour transition)
			 12  3 colour transition mars colours
			 13  landscape colours, green -> brown -> yellow
          14  yellow -> red
			 15  blue -> cyan -> green -> yellow -> brown -> white
          16  blue -> green -> red       (Chromadepth for black background)
          17  yellow -> magenta -> cyan  (Chromadepth for white background)
          18  blue -> cyan
          19  blue -> white
          20  landscape colours, modified, green -> brown -> yellow
          21  yellowish to blueish
          22  yellow to blue
			 23  red to white to blue
          24  blue -> yellow -> red
			 25  blue -> cyan -> yellow -> red
   v should lie between vmin and vmax otherwise it is clipped
   The colour components range from 0 to 1
*/
COLOUR GetColour(double v,double vmin,double vmax,int type)
{
   double dv,vmid;
	COLOUR c = {1.0,1.0,1.0};
	COLOUR c1,c2,c3;
	double ratio;

	if (vmax < vmin) {
		dv = vmin;
		vmin = vmax;
		vmax = dv;
	}
	if (vmax - vmin < 0.000001) {
		vmin -= 1;
		vmax += 1;
	}
   if (v < vmin)
      v = vmin;
   if (v > vmax)
      v = vmax;
   dv = vmax - vmin;

	switch (type) {
	case 1:
   	if (v < (vmin + 0.25 * dv)) {
      	c.r = 0;
      	c.g = 4 * (v - vmin) / dv;
			c.b = 1;
   	} else if (v < (vmin + 0.5 * dv)) {
      	c.r = 0;
			c.g = 1;
      	c.b = 1 + 4 * (vmin + 0.25 * dv - v) / dv;
   	} else if (v < (vmin + 0.75 * dv)) {
      	c.r = 4 * (v - vmin - 0.5 * dv) / dv;
			c.g = 1;
      	c.b = 0;
   	} else {
			c.r = 1;
      	c.g = 1 + 4 * (vmin + 0.75 * dv - v) / dv;
      	c.b = 0;
   	}
		break;
	case 2:
		c.r = (v - vmin) / dv;
		c.g = 0;
		c.b = (vmax - v) / dv;
		break;
	case 3:
   	c.r = (v - vmin) / dv;
   	c.b = c.r;
   	c.g = c.r;
		break;
	case 4:
      if (v < (vmin + dv / 6.0)) {
         c.r = 1; 
         c.g = 6 * (v - vmin) / dv;
         c.b = 0;
      } else if (v < (vmin + 2.0 * dv / 6.0)) {
         c.r = 1 + 6 * (vmin + dv / 6.0 - v) / dv;
         c.g = 1;
         c.b = 0;
      } else if (v < (vmin + 3.0 * dv / 6.0)) {
         c.r = 0;
         c.g = 1;
         c.b = 6 * (v - vmin - 2.0 * dv / 6.0) / dv;
      } else if (v < (vmin + 4.0 * dv / 6.0)) {
         c.r = 0;
         c.g = 1 + 6 * (vmin + 3.0 * dv / 6.0 - v) / dv;
         c.b = 1;
      } else if (v < (vmin + 5.0 * dv / 6.0)) {
         c.r = 6 * (v - vmin - 4.0 * dv / 6.0) / dv;
         c.g = 0;
         c.b = 1;
      } else {
         c.r = 1;
         c.g = 0;
         c.b = 1 + 6 * (vmin + 5.0 * dv / 6.0 - v) / dv;
      }
		break;
   case 5:
      c.r = (v - vmin) / (vmax - vmin);
      c.g = 1;
      c.b = 0;
		break;
   case 6:
      c.r = (v - vmin) / (vmax - vmin);
      c.g = (vmax - v) / (vmax - vmin);
      c.b = c.r;
		break;
   case 7:
      if (v < (vmin + 0.25 * dv)) {
         c.r = 0;
         c.g = 4 * (v - vmin) / dv;
         c.b = 1 - c.g;
      } else if (v < (vmin + 0.5 * dv)) {
			c.r = 4 * (v - vmin - 0.25 * dv) / dv;
         c.g = 1 - c.r;
         c.b = 0;
      } else if (v < (vmin + 0.75 * dv)) {
         c.g = 4 * (v - vmin - 0.5 * dv) / dv;
			c.r = 1 - c.g;
         c.b = 0;
      } else {
         c.r = 0;
         c.b = 4 * (v - vmin - 0.75 * dv) / dv;
			c.g = 1 - c.b;
      }
      break;
   case 8:
      if (v < (vmin + 0.5 * dv)) {
         c.r = 2 * (v - vmin) / dv;
         c.g = c.r;
         c.b = c.r;
      } else {
         c.r = 1 - 2 * (v - vmin - 0.5 * dv) / dv;
			c.g = c.r;
         c.b = c.r;
      }
      break;
   case 9:
      if (v < (vmin + dv / 3)) {
         c.b = 3 * (v - vmin) / dv;
			c.g = 0;
         c.r = 1 - c.b;
      } else if (v < (vmin + 2 * dv / 3)) {
         c.r = 0;
         c.g = 3 * (v - vmin - dv / 3) / dv;
         c.b = 1;
      } else {
         c.r = 3 * (v - vmin - 2 * dv / 3) / dv;
         c.g = 1 - c.r;
			c.b = 1;
      }
      break;
   case 10:
      if (v < (vmin + 0.2 * dv)) {
         c.r = 0;
         c.g = 5 * (v - vmin) / dv;
         c.b = 1;
      } else if (v < (vmin + 0.4 * dv)) {
         c.r = 0;
         c.g = 1;
         c.b = 1 + 5 * (vmin + 0.2 * dv - v) / dv;
      } else if (v < (vmin + 0.6 * dv)) {
         c.r = 5 * (v - vmin - 0.4 * dv) / dv;
         c.g = 1;
         c.b = 0;
      } else if (v < (vmin + 0.8 * dv)) {
         c.r = 1;
         c.g = 1 - 5 * (v - vmin - 0.6 * dv) / dv;
         c.b = 0;
      } else {
         c.r = 1;
         c.g = 5 * (v - vmin - 0.8 * dv) / dv;
         c.b = 5 * (v - vmin - 0.8 * dv) / dv;
      }
      break;
   case 11:
		c1.r = 200 / 255.0; c1.g =  60 / 255.0; c1.b =   0 / 255.0;
		c2.r = 250 / 255.0; c2.g = 160 / 255.0; c2.b = 110 / 255.0;
      c.r = (c2.r - c1.r) * (v - vmin) / dv + c1.r;
      c.g = (c2.g - c1.g) * (v - vmin) / dv + c1.g;
      c.b = (c2.b - c1.b) * (v - vmin) / dv + c1.b;
      break;
	case 12:
		c1.r =  55 / 255.0; c1.g =  55 / 255.0; c1.b =  45 / 255.0;
      /* c2.r = 200 / 255.0; c2.g =  60 / 255.0; c2.b =   0 / 255.0; */
		c2.r = 235 / 255.0; c2.g =  90 / 255.0; c2.b =  30 / 255.0;
		c3.r = 250 / 255.0; c3.g = 160 / 255.0; c3.b = 110 / 255.0;
		ratio = 0.4;
		vmid = vmin + ratio * dv;
		if (v < vmid) {
      	c.r = (c2.r - c1.r) * (v - vmin) / (ratio*dv) + c1.r;
      	c.g = (c2.g - c1.g) * (v - vmin) / (ratio*dv) + c1.g;
      	c.b = (c2.b - c1.b) * (v - vmin) / (ratio*dv) + c1.b;
		} else {
         c.r = (c3.r - c2.r) * (v - vmid) / ((1-ratio)*dv) + c2.r;
         c.g = (c3.g - c2.g) * (v - vmid) / ((1-ratio)*dv) + c2.g;
         c.b = (c3.b - c2.b) * (v - vmid) / ((1-ratio)*dv) + c2.b;
		}
		break;
	case 13:
      c1.r =   0 / 255.0; c1.g = 255 / 255.0; c1.b =   0 / 255.0;
      c2.r = 255 / 255.0; c2.g = 150 / 255.0; c2.b =   0 / 255.0;
      c3.r = 255 / 255.0; c3.g = 250 / 255.0; c3.b = 240 / 255.0;
      ratio = 0.3;
      vmid = vmin + ratio * dv;
      if (v < vmid) {
         c.r = (c2.r - c1.r) * (v - vmin) / (ratio*dv) + c1.r;
         c.g = (c2.g - c1.g) * (v - vmin) / (ratio*dv) + c1.g;
         c.b = (c2.b - c1.b) * (v - vmin) / (ratio*dv) + c1.b;
      } else {
         c.r = (c3.r - c2.r) * (v - vmid) / ((1-ratio)*dv) + c2.r;
         c.g = (c3.g - c2.g) * (v - vmid) / ((1-ratio)*dv) + c2.g;
         c.b = (c3.b - c2.b) * (v - vmid) / ((1-ratio)*dv) + c2.b;
      }
		break;
   case 14:
      c.r = 1;
      c.g = 1 - (v - vmin) / dv;
      c.b = 0;
      break;
   case 15:
      if (v < (vmin + 0.25 * dv)) {
         c.r = 0;
         c.g = 4 * (v - vmin) / dv;
         c.b = 1;
      } else if (v < (vmin + 0.5 * dv)) {
         c.r = 0;
         c.g = 1;
         c.b = 1 - 4 * (v - vmin - 0.25 * dv) / dv;
      } else if (v < (vmin + 0.75 * dv)) {
			c.r = 4 * (v - vmin - 0.5 * dv) / dv;
         c.g = 1;
         c.b = 0;
      } else {
         c.r = 1;
			c.g = 1;
         c.b = 4 * (v - vmin - 0.75 * dv) / dv;
      }
      break;
   case 16:
      if (v < (vmin + 0.5 * dv)) {
         c.r = 0.0;
         c.g = 2 * (v - vmin) / dv;
         c.b = 1 - 2 * (v - vmin) / dv;
      } else {
         c.r = 2 * (v - vmin - 0.5 * dv) / dv;
         c.g = 1 - 2 * (v - vmin - 0.5 * dv) / dv;
         c.b = 0.0;
      }
      break;
   case 17:
      if (v < (vmin + 0.5 * dv)) {
         c.r = 1.0;
         c.g = 1 - 2 * (v - vmin) / dv;
         c.b = 2 * (v - vmin) / dv;
      } else {
         c.r = 1 - 2 * (v - vmin - 0.5 * dv) / dv;
         c.g = 2 * (v - vmin - 0.5 * dv) / dv;
         c.b = 1.0;
      }
      break;
   case 18:
      c.r = 0;
      c.g = (v - vmin) / (vmax - vmin);
      c.b = 1;
      break;
   case 19:
      c.r = (v - vmin) / (vmax - vmin);
      c.g = c.r;
      c.b = 1;
      break;
   case 20:
      c1.r =   0 / 255.0; c1.g = 160 / 255.0; c1.b =   0 / 255.0;
      c2.r = 180 / 255.0; c2.g = 220 / 255.0; c2.b =   0 / 255.0;
      c3.r = 250 / 255.0; c3.g = 220 / 255.0; c3.b = 170 / 255.0;
      ratio = 0.3;
      vmid = vmin + ratio * dv;
      if (v < vmid) {
         c.r = (c2.r - c1.r) * (v - vmin) / (ratio*dv) + c1.r;
         c.g = (c2.g - c1.g) * (v - vmin) / (ratio*dv) + c1.g;
         c.b = (c2.b - c1.b) * (v - vmin) / (ratio*dv) + c1.b;
      } else {
         c.r = (c3.r - c2.r) * (v - vmid) / ((1-ratio)*dv) + c2.r;
         c.g = (c3.g - c2.g) * (v - vmid) / ((1-ratio)*dv) + c2.g;
         c.b = (c3.b - c2.b) * (v - vmid) / ((1-ratio)*dv) + c2.b;
      }
      break;
   case 21:
      c1.r = 255 / 255.0; c1.g = 255 / 255.0; c1.b = 200 / 255.0;
      c2.r = 150 / 255.0; c2.g = 150 / 255.0; c2.b = 255 / 255.0;
      c.r = (c2.r - c1.r) * (v - vmin) / dv + c1.r;
      c.g = (c2.g - c1.g) * (v - vmin) / dv + c1.g;
      c.b = (c2.b - c1.b) * (v - vmin) / dv + c1.b;
      break;
   case 22:
      c.r = 1 - (v - vmin) / dv;
      c.g = 1 - (v - vmin) / dv;
      c.b = (v - vmin) / dv;
      break;
   case 23:
      if (v < (vmin + 0.5 * dv)) {
         c.r = 1;
         c.g = 2 * (v - vmin) / dv;
         c.b = c.g;
      } else {
         c.r = 1 - 2 * (v - vmin - 0.5 * dv) / dv;
         c.g = c.r;
         c.b = 1;
      }
      break;
   case 24:
      if (v < (vmin + 0.5 * dv)) {
         c.r = 2 * (v - vmin) / dv;
         c.g = c.r;
         c.b = 1 - c.r;
      } else {
         c.r = 1;
         c.g = 1 - 2 * (v - vmin - 0.5 * dv) / dv;
         c.b = 0;
      }
      break;
   case 25:
      if (v < (vmin + dv / 3)) {
         c.r = 0;
         c.g = 3 * (v - vmin) / dv;
         c.b = 1;
      } else if (v < (vmin + 2 * dv / 3)) {
         c.r = 3 * (v - vmin - dv / 3) / dv;
         c.g = 1 - c.r;
         c.b = 1;
      } else {
         c.r = 1;
         c.g = 0;
         c.b = 1 - 3 * (v - vmin - 2 * dv / 3) / dv;
      }
      break;
	}
	return(c);
}

/*
	Return an ascii character whose density relates to its value
   between two bounds
*/
int AsciiColour(double v,double vmin,double vmax)
{
	int steps = 15;
	int ramp[15] = 
		{' ','.',',',':',';','-','|','I','*','o','O','%','#','@','M'};

   if (v < vmin)
      v = vmin;
   if (v > vmax)
      v = vmax;

	return(ramp[(int)(steps * (v - vmin) / (vmax - vmin))]);
}

/*
	Return the greyscale value of a colour, between 0 (black) and 1 (white)
*/
double AbsColour(COLOUR c)
{
	return(sqrt(c.r*c.r + c.b*c.b + c.g*c.g) / SQRT3);
}

/*
   Return true if two colours are the same
*/
int EqualColour(COLOUR c1,COLOUR c2)
{
	if (ABS(c1.r-c2.r) > EPSILON)
		return(FALSE);
   if (ABS(c1.g-c2.g) > EPSILON)
      return(FALSE);
   if (ABS(c1.b-c2.b) > EPSILON)
      return(FALSE);
   return(TRUE);
}

/*
   Invert a colour
*/
void InvertColour(COLOUR *c)
{
   c->r = 1 - c->r;
	c->g = 1 - c->g;
	c->b = 1 - c->b;
}

/*
   Clip a colour to legal values
*/
void ClipColour(COLOUR *c)
{
   if (c->r < 0) c->r = 0;
	if (c->r > 1) c->r = 1;
   if (c->g < 0) c->g = 0;
   if (c->g > 1) c->g = 1;
   if (c->b < 0) c->b = 0;
   if (c->b > 1) c->b = 1;
}
void ClipRGBA(RGBA *c)
{
   if (c->r < 0) c->r = 0;
   if (c->r > 1) c->r = 1;
   if (c->g < 0) c->g = 0;
   if (c->g > 1) c->g = 1;
   if (c->b < 0) c->b = 0;
   if (c->b > 1) c->b = 1;
   if (c->a < 0) c->a = 0;
   if (c->a > 1) c->a = 1;
}

/*
   Calculate HSL from RGB
   Hue is in degrees
   Lightness is betweeen 0 and 1
   Saturation is between 0 and 1
*/
HSL RGB2HSL(COLOUR c1)
{
   double themin,themax,delta;
   HSL c2;

   themin = MIN(c1.r,MIN(c1.g,c1.b));
   themax = MAX(c1.r,MAX(c1.g,c1.b));
   delta = themax - themin;
   c2.l = (themin + themax) / 2;
   c2.s = 0;
   if (c2.l > 0 && c2.l < 1)
      c2.s = delta / (c2.l < 0.5 ? (2*c2.l) : (2-2*c2.l));
   c2.h = 0;
   if (delta > 0) {
      if (themax == c1.r && themax != c1.g)
         c2.h += (c1.g - c1.b) / delta;
      if (themax == c1.g && themax != c1.b)
         c2.h += (2 + (c1.b - c1.r) / delta);
      if (themax == c1.b && themax != c1.r)
         c2.h += (4 + (c1.r - c1.g) / delta);
      c2.h *= 60;
		if (c2.h < 0)
			c2.h += 360;
   }
   return(c2);
}

/*
   Calculate RGB from HSL, reverse of RGB2HSL()
   Hue is in degrees
   Lightness is between 0 and 1
   Saturation is between 0 and 1
*/
COLOUR HSL2RGB(HSL c1)
{
   COLOUR c2,sat,ctmp;

   while (c1.h < 0)
      c1.h += 360;
   while (c1.h > 360)
      c1.h -= 360;

   if (c1.h < 120) {
      sat.r = (120 - c1.h) / 60.0;
      sat.g = c1.h / 60.0;
      sat.b = 0;
   } else if (c1.h < 240) {
      sat.r = 0;
      sat.g = (240 - c1.h) / 60.0;
      sat.b = (c1.h - 120) / 60.0;
   } else {
      sat.r = (c1.h - 240) / 60.0;
      sat.g = 0;
      sat.b = (360 - c1.h) / 60.0;
   }
   sat.r = MIN(sat.r,1);
   sat.g = MIN(sat.g,1);
   sat.b = MIN(sat.b,1);

   ctmp.r = 2 * c1.s * sat.r + (1 - c1.s);
   ctmp.g = 2 * c1.s * sat.g + (1 - c1.s);
   ctmp.b = 2 * c1.s * sat.b + (1 - c1.s);

   if (c1.l < 0.5) {
      c2.r = c1.l * ctmp.r;
      c2.g = c1.l * ctmp.g;
      c2.b = c1.l * ctmp.b;
   } else {
      c2.r = (1 - c1.l) * ctmp.r + 2 * c1.l - 1;
      c2.g = (1 - c1.l) * ctmp.g + 2 * c1.l - 1;
      c2.b = (1 - c1.l) * ctmp.b + 2 * c1.l - 1;
   }

   return(c2);
}

/*
	Compute YCC from RGB
	Luminance is between 0 and 1
	Both chrominance channels are between -0.5 and 0.5
*/
YCC RGB2YCC(COLOUR c1)
{
	YCC c2;

	c2.y  =  0.2989 * c1.r + 0.5866 * c1.g + 0.1145 * c1.b;
	c2.cr = -0.1687 * c1.r - 0.3312 * c1.g + 0.5000 * c1.b;
	c2.cb =  0.5000 * c1.r - 0.4183 * c1.g - 0.0816 * c1.b;

	return(c2);
}

/*
   Compute RGB from YCC
	Reverse of RGB2YCC()
*/
COLOUR YCC2RGB(YCC c1)
{
	COLOUR c2;

	c2.r = c1.y                  + 1.4022 * c1.cb;
	c2.g = c1.y - 0.3456 * c1.cr - 0.7145 * c1.cb;
	c2.b = c1.y + 1.7710 * c1.cr;
	if (c2.r < 0) c2.r = 0;
   if (c2.g < 0) c2.g = 0;
   if (c2.b < 0) c2.b = 0;
   if (c2.r > 1) c2.r = 1;
   if (c2.g > 1) c2.g = 1;
   if (c2.b > 1) c2.b = 1;

	return(c2);
}

/*
   Calculate RGB from HSV, reverse of RGB2HSV()
   Hue is in degrees
   Lightness is between 0 and 1
   Saturation is between 0 and 1
*/
COLOUR HSV2RGB(HSV c1)
{
   COLOUR c2,sat;

   while (c1.h < 0)
      c1.h += 360;
   while (c1.h > 360)
      c1.h -= 360;

   if (c1.h < 120) {
      sat.r = (120 - c1.h) / 60.0;
      sat.g = c1.h / 60.0;
      sat.b = 0;
   } else if (c1.h < 240) {
      sat.r = 0;
      sat.g = (240 - c1.h) / 60.0;
      sat.b = (c1.h - 120) / 60.0;
   } else {
      sat.r = (c1.h - 240) / 60.0;
      sat.g = 0;
      sat.b = (360 - c1.h) / 60.0;
   }
   sat.r = MIN(sat.r,1);
   sat.g = MIN(sat.g,1);
   sat.b = MIN(sat.b,1);

   c2.r = (1 - c1.s + c1.s * sat.r) * c1.v;
   c2.g = (1 - c1.s + c1.s * sat.g) * c1.v;
   c2.b = (1 - c1.s + c1.s * sat.b) * c1.v;

   return(c2);
}

/*
   Calculate HSV from RGB
   Hue is in degrees
   Lightness is betweeen 0 and 1
   Saturation is between 0 and 1
*/
HSV RGB2HSV(COLOUR c1)
{
   double themin,themax,delta;
   HSV c2;

   themin = MIN(c1.r,MIN(c1.g,c1.b));
   themax = MAX(c1.r,MAX(c1.g,c1.b));
   delta = themax - themin;
   c2.v = themax;
   c2.s = 0;
   if (themax > 0)
      c2.s = delta / themax;
   c2.h = 0;
   if (delta > 0) {
      if (themax == c1.r && themax != c1.g)
         c2.h += (c1.g - c1.b) / delta;
      if (themax == c1.g && themax != c1.b)
         c2.h += (2 + (c1.b - c1.r) / delta);
      if (themax == c1.b && themax != c1.r)
         c2.h += (4 + (c1.r - c1.g) / delta);
      c2.h *= 60;
		if (c2.h < 0)
			c2.h += 360;
   }
   return(c2);
}

/*-------------------------------------------------------------------------
	Return the plane parameters a,b,c from vertices p[]
*/
void PlaneParam(XYZ *p,double *a,double *b,double *c,double *d)
{
   *a = p[0].y * (p[1].z - p[2].z) +
        p[1].y * (p[2].z - p[0].z) +
        p[2].y * (p[0].z - p[1].z);
   *b = p[0].z * (p[1].x - p[2].x) +
        p[1].z * (p[2].x - p[0].x) +
        p[2].z * (p[0].x - p[1].x);
   *c = p[0].x * (p[1].y - p[2].y) +
        p[1].x * (p[2].y - p[0].y) +
        p[2].x * (p[0].y - p[1].y);
   *d = - p[0].x * (p[1].y * p[2].z - p[2].y * p[1].z)
        - p[1].x * (p[2].y * p[0].z - p[0].y * p[2].z)
        - p[2].x * (p[0].y * p[1].z - p[1].y * p[0].z);
}
	
/*
	Given a quad defined by p[0], p[1], p[2], p[3]
	and two ratios mua and mub
	Return the interpolated point

     p[0]    -mux->
     +-----
     |     -----     p[1]
   |  |         -----+
  muy  |             |
   |    |             |
   v    +-----         |
        p[2]  -----     |
                   -----+
                        p[3]
*/
XY PointInPolygon(XY *p,double mux,double muy)
{
	double x1,y1,x2,y2,x3,y3,x4,y4;
	double denom,numer;
	XY q;

   x1 = p[0].x + mux * (p[1].x - p[0].x);
   y1 = p[0].y + mux * (p[1].y - p[0].y);
   x2 = p[2].x + mux * (p[3].x - p[2].x);
   y2 = p[2].y + mux * (p[3].y - p[2].y);
   x3 = p[0].x + muy * (p[2].x - p[0].x);
   y3 = p[0].y + muy * (p[2].y - p[0].y);
   x4 = p[1].x + muy * (p[3].x - p[1].x);
   y4 = p[1].y + muy * (p[3].y - p[1].y);
   numer = (x4 - x3) * (y1 - y3) - (y4 - y3) * (x1 - x3);
   denom = ((y4 - y3) * (x2 - x1) - (x4 - x3) * (y2 - y1));
	// mu = numer / denom;
	if (ABS(denom) < EPS) { // Should not happen in proper usage
		q.x = 0;
		q.y = 0;
	} else {
   	q.x = x1 + numer * (x2 - x1) / denom;
   	q.y = y1 + numer * (y2 - y1) / denom;
	}

	return(q);
}

/*-------------------------------------------------------------------------
	Return TRUE if the point "q" lies within the bounded polygon
	p[] with "np" vertices, assume the polygon is convex!
	In any case the distance "dist" is returned
*/
int PointIn3DPolygon(int np,XYZ *p,XYZ q,double *dist)
{
	int i,sign=0;
	double a,b,c,d;
	double dd,denom;
	XYZ n,p1[3];

	/* Determine the plane parameters (the normal is a,b,c) */
	PlaneParam(p,&n.x,&n.y,&n.z,&d);

	*dist = 0;
	denom = sqrt(n.x*n.x + n.y*n.y + n.z*n.z);
	if (ABS(denom) < EPSILON)		/* Degenerate plane */
		return(FALSE);
	*dist = (n.x * q.x + n.y * q.y + n.z * q.z + d) / denom;

	for (i=0;i<np;i++) {
		p1[0] = p[i]; 
		p1[1] = p[(i+1)%np]; 
		p1[2].x = p[i].x + n.x;
   	p1[2].y = p[i].y + n.y;
   	p1[2].z = p[i].z + n.z;
		PlaneParam(p1,&a,&b,&c,&d);
		dd = (a * q.x + b * q.y + c * q.z + d) / sqrt(a*a + b*b + c*c);
		if (i == 0)
			sign = SIGN(dd);
		else if (sign != SIGN(dd))
      	return(FALSE);
	}

	return(TRUE);
}

/*-------------------------------------------------------------------------
   Determines which side of a line the point (x,y) lies.
   The line goes from (x1,y1) to (x2,y2)
   Return codes are -1 for points to the left
                     0 for points on the line
                    +1 for points to the right
*/
int WhichSide(double x,double y,double x1,double y1,double x2,double y2)
{
    double dist;

    dist = (y - y1) * (x2 - x1) - (y2 - y1) * (x - x1);
    if (dist > 0)
       return(-1);
    else if (dist < 0)
       return(1);
    else
       return(0);
}

/*-------------------------------------------------------------------------
   Given a plane through the points p1,p2,p3 determine the intersection 
   of the perpendicular from the point (xp,yp) on the x-y plane.
*/
double PlanePoint(double xp,double yp,XYZ p1,XYZ p2,XYZ p3)
{
   double a,b,c,d;

   a = p1.y * (p2.z - p3.z) + p2.y * (p3.z - p1.z) + p3.y * (p1.z - p2.z);
   b = p1.z * (p2.x - p3.x) + p2.z * (p3.x - p1.x) + p3.z * (p1.x - p2.x);
   c = p1.x * (p2.y - p3.y) + p2.x * (p3.y - p1.y) + p3.x * (p1.y - p2.y);
   d = - p1.x * (p2.y * p3.z - p3.y * p2.z) 
        - p2.x * (p3.y * p1.z - p1.y * p3.z) 
        - p3.x * (p1.y * p2.z - p2.y * p1.z);
        
   if (ABS(c) > EPSILON)
      return(- (a * xp + b * yp + d) / c);
   else
      return(0.0);
}

/*-------------------------------------------------------------------------
   Returns TRUE if the point (xp,yp) lies inside the projection of the
   triangle (p1,p2,p3) onto the x-y plane
   A point is in the centre if it is on the same side of all the edges
   or if it lies on one of the edges.
*/
int InTriangle(double xp,double yp,XYZ p1,XYZ p2,XYZ p3)
{
   int side1,side2,side3;

   side1 = WhichSide(xp,yp,p1.x,p1.y,p2.x,p2.y);
   side2 = WhichSide(xp,yp,p2.x,p2.y,p3.x,p3.y);
   side3 = WhichSide(xp,yp,p3.x,p3.y,p1.x,p1.y);

   if (side1 == 0 && side2 == 0)
      return(TRUE);
   if (side2 == 0 && side3 == 0)
      return(TRUE);

   if (side1 == 0 && (side2 == side3))
      return(TRUE);
   if (side2 == 0 && (side1 == side3))
      return(TRUE);
   if (side3 == 0 && (side1 == side2))
      return(TRUE);
        
   if ((side1 == side2) && (side2 == side3))
      return(TRUE);
                
   return(FALSE);
}

/*-------------------------------------------------------------------------
	Calculate the length of a polygon
*/
double PolygonLength(XYZ *p,int np)
{
   int i;
   double length = 0;

   for (i=0;i<np-1;i++)
      length += VectorLength(p[i],p[i+1]);

   return(length);
}

/*-------------------------------------------------------------------------
   Calculate a curvature measure of a polygon
*/
double PolygonCurve(XYZ *p,int np)
{
   int i;
   double c,cmin=1e32,cmax=-1e32;
   XYZ p1,p2;

   if (np < 3)
      return(0.0);

   p1.x = p[1].x - p[0].x;
   p1.y = p[1].y - p[0].y;
   p1.z = p[1].z - p[0].z;

   for (i=2;i<np;i++) {
      p2.x = p[i].x - p[0].x;
      p2.y = p[i].y - p[0].y;
      p2.z = p[i].z - p[0].z;
      c = VectorAngle(p1,p2);
      cmax = MAX(c,cmax);
      cmin = MIN(c,cmin);
   }

   return((cmax - cmin)*RTOD);
}

/*-------------------------------------------------------------------------
   Given a grid cell and an isolevel, calculate the triangular
   facets requied to represent the isosurface through the cell.
   Return the number of triangular facets, the array "triangles"
   will be loaded up with the vertices at most 5 triangular facets.
   0 will be returned if the grid cell is either totally above
   of totally below the isolevel.
*/
int PolygoniseCube(GRIDCELL g,double iso,TRIANGLE *tri)
{
	int i,ntri = 0;
	int cubeindex;
	XYZ vertlist[12];
/*
	int edgeTable[256].  It corresponds to the 2^8 possible combinations of
	of the eight (n) vertices either existing inside or outside (2^n) of the
	surface.  A vertex is inside of a surface if the value at that vertex is
	less than that of the surface you are scanning for.  The table index is
	constructed bitwise with bit 0 corresponding to vertex 0, bit 1 to vert
	1.. bit 7 to vert 7.  The value in the table tells you which edges of
	the table are intersected by the surface.  Once again bit 0 corresponds
	to edge 0 and so on, up to edge 12. 
	Constructing the table simply consisted of having a program run thru
	the 256 cases and setting the edge bit if the vertices at either end of
	the edge had different values (one is inside while the other is out). 
	The purpose of the table is to speed up the scanning process.  Only the
	edges whose bit's are set contain vertices of the surface.
	Vertex 0 is on the bottom face, back edge, left side.  
	The progression of vertices is clockwise around the bottom face
	and then clockwise around the top face of the cube.  Edge 0 goes from
	vertex 0 to vertex 1, Edge 1 is from 2->3 and so on around clockwise to
	vertex 0 again. Then Edge 4 to 7 make up the top face, 4->5, 5->6, 6->7
	and 7->4.  Edge 8 thru 11 are the vertical edges from vert 0->4, 1->5,
	2->6, and 3->7.
	    4--------5     *---4----*
 	   /|       /|    /|       /|    
 	  / |      / |   7 |      5 |    
 	 /  |     /  |  /  8     /  9    
 	7--------6   | *----6---*   | 
 	|   |    |   | |   |    |   |
 	|   0----|---1 |   *---0|---*  
 	|  /     |  /  11 /     10 /   
 	| /      | /   | 3      | 1
 	|/       |/    |/       |/    
 	3--------2     *---2----*
*/
int edgeTable[256]={
0x0  , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0   };

/*
	int triTable[256][16] also corresponds to the 256 possible combinations
	of vertices.
	The [16] dimension of the table is again the list of edges of the cube
	which are intersected by the surface.  This time however, the edges are
	enumerated in the order of the vertices making up the triangle mesh of
	the surface.  Each edge contains one vertex that is on the surface. 
	Each triple of edges listed in the table contains the vertices of one
	triangle on the mesh.  The are 16 entries because it has been shown that
	there are at most 5 triangles in a cube and each "edge triple" list is
	terminated with the value -1. 
	For example triTable[3] contains 
	{1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
	This corresponds to the case of a cube whose vertex 0 and 1 are inside
	of the surface and the rest of the verts are outside (00000001 bitwise
	OR'ed with 00000010 makes 00000011 == 3).  Therefore, this cube is
	intersected by the surface roughly in the form of a plane which cuts
	edges 8,9,1 and 3.  This quadrilateral can be constructed from two
	triangles: one which is made of the intersection vertices found on edges
	1,8, and 3; the other is formed from the vertices on edges 9,8, and 1. 
	Remember, each intersected edge contains only one surface vertex.  The
	vertex triples are listed in counter clockwise order for proper facing.
*/
int triTable[256][16] =
{{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
{3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
{3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
{3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
{9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
{2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
{8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
{4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
{3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
{1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
{4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
{4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
{5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
{2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
{9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
{0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
{2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
{10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
{5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
{5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
{9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
{10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
{8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
{2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
{7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
{2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
{11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
{5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
{11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
{11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
{9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
{2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
{6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
{3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
{6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
{10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
{6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
{8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
{7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
{3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
{0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
{9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
{8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
{5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
{0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
{6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
{10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
{10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
{8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
{1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
{0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
{10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
{3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
{6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
{9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
{8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
{3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
{6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
{0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
{10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
{10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
{2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
{7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
{7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
{2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
{1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
{11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
{8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
{0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
{7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
{10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
{2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
{6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
{7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
{2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
{10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
{10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
{0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
{7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
{6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
{8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
{9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
{6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
{4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
{10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
{8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
{0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
{1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
{8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
{10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
{4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
{10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
{11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
{9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
{6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
{7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
{3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
{7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
{3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
{6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
{9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
{1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
{4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
{7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
{6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
{3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
{0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
{6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
{0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
{11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
{6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
{5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
{9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
{1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
{10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
{0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
{5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
{10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
{11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
{9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
{7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
{2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
{8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
{9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
{9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
{1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
{9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
{5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
{0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
{10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
{2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
{0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
{0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
{9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
{5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
{3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
{5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
{8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
{0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
{9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
{1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
{3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
{4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
{9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
{11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
{11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
{2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
{9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
{3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
{1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
{4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
{3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
{0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
{1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}};

   /*
      Determine the index into the edge table which
      tells us which vertices are inside of the surface
   */
   cubeindex = 0;
   if (g.val[0] < iso) cubeindex |= 1;
   if (g.val[1] < iso) cubeindex |= 2;
   if (g.val[2] < iso) cubeindex |= 4;
   if (g.val[3] < iso) cubeindex |= 8;
   if (g.val[4] < iso) cubeindex |= 16;
   if (g.val[5] < iso) cubeindex |= 32;
   if (g.val[6] < iso) cubeindex |= 64;
   if (g.val[7] < iso) cubeindex |= 128;

   /* Cube is entirely in/out of the surface */
   if (edgeTable[cubeindex] == 0)
      return(0);

   /* Find the vertices where the surface intersects the cube */
   if (edgeTable[cubeindex] & 1) {
      vertlist[0] = VertexInterp(iso,g.p[0],g.p[1],g.val[0],g.val[1]);
	}
   if (edgeTable[cubeindex] & 2) {
      vertlist[1] = VertexInterp(iso,g.p[1],g.p[2],g.val[1],g.val[2]);
	}
   if (edgeTable[cubeindex] & 4) {
      vertlist[2] = VertexInterp(iso,g.p[2],g.p[3],g.val[2],g.val[3]);
	}
   if (edgeTable[cubeindex] & 8) {
      vertlist[3] = VertexInterp(iso,g.p[3],g.p[0],g.val[3],g.val[0]);
	}
   if (edgeTable[cubeindex] & 16) {
      vertlist[4] = VertexInterp(iso,g.p[4],g.p[5],g.val[4],g.val[5]);
	}
   if (edgeTable[cubeindex] & 32) {
      vertlist[5] = VertexInterp(iso,g.p[5],g.p[6],g.val[5],g.val[6]);
	}
   if (edgeTable[cubeindex] & 64) {
      vertlist[6] = VertexInterp(iso,g.p[6],g.p[7],g.val[6],g.val[7]);
	}
   if (edgeTable[cubeindex] & 128) {
      vertlist[7] = VertexInterp(iso,g.p[7],g.p[4],g.val[7],g.val[4]);
	}
   if (edgeTable[cubeindex] & 256) {
      vertlist[8] = VertexInterp(iso,g.p[0],g.p[4],g.val[0],g.val[4]);
	}
   if (edgeTable[cubeindex] & 512) {
      vertlist[9] = VertexInterp(iso,g.p[1],g.p[5],g.val[1],g.val[5]);
	}
   if (edgeTable[cubeindex] & 1024) {
      vertlist[10] = VertexInterp(iso,g.p[2],g.p[6],g.val[2],g.val[6]);
	}
   if (edgeTable[cubeindex] & 2048) {
      vertlist[11] = VertexInterp(iso,g.p[3],g.p[7],g.val[3],g.val[7]);
	}

   /* Create the triangles */
   for (i=0;triTable[cubeindex][i]!=-1;i+=3) {
      tri[ntri].p[0] = vertlist[triTable[cubeindex][i  ]];
      tri[ntri].p[1] = vertlist[triTable[cubeindex][i+1]];
      tri[ntri].p[2] = vertlist[triTable[cubeindex][i+2]];
      ntri++;
   }

	return(ntri);
}

/*
   Polygonise a tetrahedron given its vertices within a cube
	This is an alternative algorithm to polygonisegrid.
	It results in a smoother surface but more triangular facets.

                      + 0
                     /|\
                    / | \
                   /  |  \
                  /   |   \
                 /    |    \
                /     |     \
               +-------------+ 1
              3 \     |     /
                 \    |    /
                  \   |   /
                   \  |  /
                    \ | /
                     \|/
                      + 2

	It's main purpose is still to polygonise a gridded dataset and
	would normally be called 6 times, one for each tetrahedron making
	up the grid cell.
	Given the grid labelling as in PolygniseGrid one would call
		PolygoniseTri(grid,iso,triangles,0,2,3,7);
		PolygoniseTri(grid,iso,triangles,0,2,6,7);
		PolygoniseTri(grid,iso,triangles,0,4,6,7);
		PolygoniseTri(grid,iso,triangles,0,6,1,2);
		PolygoniseTri(grid,iso,triangles,0,6,1,4);
		PolygoniseTri(grid,iso,triangles,5,6,1,4);
*/
int PolygoniseTri(GRIDCELL g,double iso,
	TRIANGLE *tri,int v0,int v1,int v2,int v3)
{
   int ntri = 0;
   int triindex;
  
   /*
      Determine which of the 16 cases we have given which vertices
      are above or below the isosurface
   */
   triindex = 0;
   if (g.val[v0] < iso) triindex |= 1;
   if (g.val[v1] < iso) triindex |= 2;
   if (g.val[v2] < iso) triindex |= 4;
   if (g.val[v3] < iso) triindex |= 8;

   /* Form the vertices of the triangles for each case */
   switch (triindex) {
   case 0x00:
   case 0x0F:
      break;
   case 0x0E:
   case 0x01:
      tri[0].p[0] = VertexInterp(iso,g.p[v0],g.p[v1],g.val[v0],g.val[v1]);
      tri[0].p[1] = VertexInterp(iso,g.p[v0],g.p[v2],g.val[v0],g.val[v2]);
      tri[0].p[2] = VertexInterp(iso,g.p[v0],g.p[v3],g.val[v0],g.val[v3]);
      ntri++;
      break;
   case 0x0D:
   case 0x02:
      tri[0].p[0] = VertexInterp(iso,g.p[v1],g.p[v0],g.val[v1],g.val[v0]);
      tri[0].p[1] = VertexInterp(iso,g.p[v1],g.p[v3],g.val[v1],g.val[v3]);
      tri[0].p[2] = VertexInterp(iso,g.p[v1],g.p[v2],g.val[v1],g.val[v2]);
      ntri++;
      break;
   case 0x0C:
   case 0x03:
      tri[0].p[0] = VertexInterp(iso,g.p[v0],g.p[v3],g.val[v0],g.val[v3]);
      tri[0].p[1] = VertexInterp(iso,g.p[v0],g.p[v2],g.val[v0],g.val[v2]);
      tri[0].p[2] = VertexInterp(iso,g.p[v1],g.p[v3],g.val[v1],g.val[v3]);
      ntri++;
      tri[1].p[0] = tri[0].p[2];
      tri[1].p[1] = VertexInterp(iso,g.p[v1],g.p[v2],g.val[v1],g.val[v2]);
      tri[1].p[2] = tri[0].p[1];
      ntri++;
      break;
   case 0x0B:
   case 0x04:
      tri[0].p[0] = VertexInterp(iso,g.p[v2],g.p[v0],g.val[v2],g.val[v0]);
      tri[0].p[1] = VertexInterp(iso,g.p[v2],g.p[v1],g.val[v2],g.val[v1]);
      tri[0].p[2] = VertexInterp(iso,g.p[v2],g.p[v3],g.val[v2],g.val[v3]);
      ntri++;
      break;
   case 0x0A:
   case 0x05:
      tri[0].p[0] = VertexInterp(iso,g.p[v0],g.p[v1],g.val[v0],g.val[v1]);
      tri[0].p[1] = VertexInterp(iso,g.p[v2],g.p[v3],g.val[v2],g.val[v3]);
      tri[0].p[2] = VertexInterp(iso,g.p[v0],g.p[v3],g.val[v0],g.val[v3]);
      ntri++;
      tri[1].p[0] = tri[0].p[0];
      tri[1].p[1] = VertexInterp(iso,g.p[v1],g.p[v2],g.val[v1],g.val[v2]);
      tri[1].p[2] = tri[0].p[1];
      ntri++;
      break;
   case 0x09:
   case 0x06:
      tri[0].p[0] = VertexInterp(iso,g.p[v0],g.p[v1],g.val[v0],g.val[v1]);
      tri[0].p[1] = VertexInterp(iso,g.p[v1],g.p[v3],g.val[v1],g.val[v3]);
      tri[0].p[2] = VertexInterp(iso,g.p[v2],g.p[v3],g.val[v2],g.val[v3]);
      ntri++;
      tri[1].p[0] = tri[0].p[0];
      tri[1].p[1] = VertexInterp(iso,g.p[v0],g.p[v2],g.val[v0],g.val[v2]);
      tri[1].p[2] = tri[0].p[2];
      ntri++;
      break;
   case 0x07:
   case 0x08:
      tri[0].p[0] = VertexInterp(iso,g.p[v3],g.p[v0],g.val[v3],g.val[v0]);
      tri[0].p[1] = VertexInterp(iso,g.p[v3],g.p[v2],g.val[v3],g.val[v2]);
      tri[0].p[2] = VertexInterp(iso,g.p[v3],g.p[v1],g.val[v3],g.val[v1]);
      ntri++;
      break;
   }

   return(ntri);
}

/*-------------------------------------------------------------------------
   Return the point between two points in the same ratio as
   isolevel is between valp1 and valp2
*/
XYZ VertexInterp(double isolevel,XYZ p1,XYZ p2,double valp1,double valp2)
{
   double mu;
   XYZ p;

   if (ABS(isolevel-valp1) < 0.00001)
      return(p1);
   if (ABS(isolevel-valp2) < 0.00001)
      return(p2);
   if (ABS(valp1-valp2) < 0.00001)
      return(p1);
   mu = (isolevel - valp1) / (valp2 - valp1);
   p.x = p1.x + mu * (p2.x - p1.x);
   p.y = p1.y + mu * (p2.y - p1.y);
   p.z = p1.z + mu * (p2.z - p1.z);

   return(p);
}

/*
	Compare two 3 or 4 vertex polygons
	Perform the following tests 
	If there is nothing interesting about their relationship return 0
	If the polygons are identical, return 1
*/
int PolygonCompare(POLY34 p1,POLY34 p2) 
{
	
	/* See if the polygons are identical */
	if (p1.np == 3 && p2.np == 3) {
		if (VectorEqual(p1.p[0],p2.p[0]) || 
			 VectorEqual(p1.p[0],p2.p[1]) ||
			 VectorEqual(p1.p[0],p2.p[2])) {
      	if (VectorEqual(p1.p[1],p2.p[0]) ||
          	 VectorEqual(p1.p[1],p2.p[1]) ||
          	 VectorEqual(p1.p[1],p2.p[2])) {
         	if (VectorEqual(p1.p[2],p2.p[0]) ||
             	 VectorEqual(p1.p[2],p2.p[1]) ||
             	 VectorEqual(p1.p[2],p2.p[2])) {
					return(1);
				}
			}
      }
	}
   if (p1.np == 4 && p2.np == 4) {
      if (VectorEqual(p1.p[0],p2.p[0]) ||
          VectorEqual(p1.p[0],p2.p[1]) ||
          VectorEqual(p1.p[0],p2.p[2]) ||
			 VectorEqual(p1.p[0],p2.p[3])) {
         if (VectorEqual(p1.p[1],p2.p[0]) ||
             VectorEqual(p1.p[1],p2.p[1]) ||
             VectorEqual(p1.p[1],p2.p[2]) ||
				 VectorEqual(p1.p[1],p2.p[3])) {
            if (VectorEqual(p1.p[2],p2.p[0]) ||
                VectorEqual(p1.p[2],p2.p[1]) ||
                VectorEqual(p1.p[2],p2.p[2]) ||
					 VectorEqual(p1.p[2],p2.p[3])) {
            	if (VectorEqual(p1.p[3],p2.p[0]) ||
                	 VectorEqual(p1.p[3],p2.p[1]) ||
                	 VectorEqual(p1.p[3],p2.p[2]) ||
						 VectorEqual(p1.p[3],p2.p[3])) {
               	return(1);
					}
            }
         }
      }
   }

	return(0);
}

/*
   Create a triangular facet approximation to a sphere
	Unit radius
   Return the number of facets created.
   The number of facets will be (4^iterations) * 8
*/
int CreateNSphere(TRIFACE *f,int iterations)
{
   int i,it;
   double a;
   XYZ p[6] = {{0,0,1},  {0,0,-1},  {-1,-1,0},  {1,-1,0},  {1,1,0}, {-1,1,0}};
   XYZ pa,pb,pc;
   int nt = 0,ntold;

   /* Create the level 0 object */
   a = 1 / sqrt(2.0);
   for (i=0;i<6;i++) {
      p[i].x *= a;
      p[i].y *= a;
   }
   f[0].p[0] = p[0]; f[0].p[1] = p[3]; f[0].p[2] = p[4];
   f[1].p[0] = p[0]; f[1].p[1] = p[4]; f[1].p[2] = p[5];
   f[2].p[0] = p[0]; f[2].p[1] = p[5]; f[2].p[2] = p[2];
   f[3].p[0] = p[0]; f[3].p[1] = p[2]; f[3].p[2] = p[3];
   f[4].p[0] = p[1]; f[4].p[1] = p[4]; f[4].p[2] = p[3];
   f[5].p[0] = p[1]; f[5].p[1] = p[5]; f[5].p[2] = p[4];
   f[6].p[0] = p[1]; f[6].p[1] = p[2]; f[6].p[2] = p[5];
   f[7].p[0] = p[1]; f[7].p[1] = p[3]; f[7].p[2] = p[2];
   nt = 8;

   if (iterations < 1)
      return(nt);

   /* Bisect each edge and move to the surface of a unit sphere */
   for (it=0;it<iterations;it++) {
      ntold = nt;
      for (i=0;i<ntold;i++) {
			pa = MidPoint(f[i].p[0],f[i].p[1]);
			pb = MidPoint(f[i].p[1],f[i].p[2]);
			pc = MidPoint(f[i].p[2],f[i].p[0]);
         Normalise(&pa);
         Normalise(&pb);
         Normalise(&pc);
         f[nt].p[0] = f[i].p[0]; f[nt].p[1] = pa;        f[nt].p[2] = pc; nt++;
         f[nt].p[0] = pa;        f[nt].p[1] = f[i].p[1]; f[nt].p[2] = pb; nt++;
         f[nt].p[0] = pb;        f[nt].p[1] = f[i].p[2]; f[nt].p[2] = pc; nt++;
         f[i].p[0] = pa;
         f[i].p[1] = pb;
         f[i].p[2] = pc;
      }
   }

   return(nt);
}

/* GEOMETRIC OUTPUT -----------------------------------------------------*/

/*
	Write a Radiance facet 
*/
void WriteRadFacet(FILE *fptr,XYZ *p,int n,char *mat)
{
   int i;

   if (n <= 2)
      return;

    fprintf(fptr,"%s polygon p\n",mat);
    fprintf(fptr,"0\n0\n%d",3*n);
    for (i=0;i<n;i++)
       fprintf(fptr,"\t%g %g %g\n",p[i].x,p[i].y,p[i].z);  
}

/*
   Write a DXF header
*/
void WriteDXFHeader(FILE *fptr,char *s)
{
   fprintf(fptr,"999\nDXF from %s\n",s);
   fprintf(fptr,"  0\nSECTION\n");
   fprintf(fptr,"  2\nHEADER\n");
   fprintf(fptr,"  9\n$ACADVER\n");
   fprintf(fptr,"  1\nAC1006\n");
   fprintf(fptr,"  9\n$INSBASE\n");
   fprintf(fptr,"  10\n0.0\n");
   fprintf(fptr,"  20\n0.0\n");
   fprintf(fptr,"  30\n0.0\n");
   fprintf(fptr,"  9\n$EXTMIN\n");
   fprintf(fptr,"  10\n0.0\n");
   fprintf(fptr,"  20\n0.0\n");
   fprintf(fptr,"  9\n$EXTMAX\n");
   fprintf(fptr,"  10\n1000.0\n");
   fprintf(fptr,"  20\n1000.0\n");
   fprintf(fptr,"  0\nENDSEC\n");

	fprintf(fptr,"  0\nSECTION\n");
	fprintf(fptr,"  2\nTABLES\n");
   fprintf(fptr,"  0\nTABLE\n");
   fprintf(fptr,"  2\nLTYPE\n");
   fprintf(fptr,"  70\n1\n");
   fprintf(fptr,"  0\nLTYPE\n");
   fprintf(fptr,"  2\niONTINUOUS\n");
   fprintf(fptr,"  70\n64\n");
   fprintf(fptr,"  3\nSolid line\n");
   fprintf(fptr,"  72\n65\n");
   fprintf(fptr,"  73\n0\n");
   fprintf(fptr,"  40\n0.000000\n");
   fprintf(fptr,"  0\nENDTAB\n");
   fprintf(fptr,"  0\nTABLE\n");
   fprintf(fptr,"  2\nLAYER\n");
   fprintf(fptr,"  70\n6\n");
   fprintf(fptr,"  0\nLAYER\n");
   fprintf(fptr,"  2\n1\n");
   fprintf(fptr,"  70\n64\n");
   fprintf(fptr,"  62\n7\n");
   fprintf(fptr,"  6\nCONTINUOUS\n");
   fprintf(fptr,"  0\nLAYER\n");
   fprintf(fptr,"  2\n2\n");
   fprintf(fptr,"  70\n64\n");
   fprintf(fptr,"  62\n7\n");
   fprintf(fptr,"  6\nCONTINUOUS\n");
   fprintf(fptr,"  0\nENDTAB\n");
   fprintf(fptr,"  0\nTABLE\n");
   fprintf(fptr,"  2\nSTYLE\n");
   fprintf(fptr,"  70\n0\n");
   fprintf(fptr,"  0\nENDTAB\n");
	fprintf(fptr,"  0\nENDSEC\n"); 

   fprintf(fptr,"  0\nSECTION\n");
	fprintf(fptr,"  2\nBLOCKS\n");
	fprintf(fptr,"  0\nENDSEC\n");

	fprintf(fptr,"  0\nSECTION\n");
   fprintf(fptr,"  2\nENTITIES\n");
}

/*
   Write a DXF footer
*/
void WriteDXFFooter(FILE *fptr)
{
    fprintf(fptr,"  0\nENDSEC\n");
    fprintf(fptr,"  0\nEOF\n");
}

/*
   Write a DXF facet
*/
void WriteDXFFacet(FILE *fptr,XYZ *p,int n)
{
	int i;

	if (n < 2)
		return;

	switch (n) {
	case 2:
		fprintf(fptr,"  0\nLINE\n");
		fprintf(fptr,"  8\nDefault\n");
		for (i=0;i<n;i++) {
         fprintf(fptr,"%3d\n%g\n",10+i,p[i].x);  
         fprintf(fptr,"%3d\n%g\n",20+i,p[i].y);  
         fprintf(fptr,"%3d\n%g\n",30+i,p[i].z);  
		}
		break;
	case 3:
	case 4:
		fprintf(fptr,"  0\n3DFACE\n");
		fprintf(fptr,"  8\nDefault\n");
		for (i=0;i<n;i++) {
			fprintf(fptr,"%3d\n%g\n",10+i,p[i].x);	
         fprintf(fptr,"%3d\n%g\n",20+i,p[i].y);  
         fprintf(fptr,"%3d\n%g\n",30+i,p[i].z);  
		}
		if (n == 3) {
         fprintf(fptr,"%3d\n%g\n",13,p[2].x);  
         fprintf(fptr,"%3d\n%g\n",23,p[2].y);  
         fprintf(fptr,"%3d\n%g\n",33,p[2].z);  
		}
		break;
	default:
		break;
	}
}

/* SYSTEM SPECIFIC  --------------------------------------------------*/

/*
	Time scale at microsecond resolution but returned as seconds
*/
double GetRunTime(void)
{
#ifndef WIN32
	double sec = 0;
	struct timeval tp;

   gettimeofday(&tp,NULL);
   sec = tp.tv_sec + tp.tv_usec / 1000000.0;
#else
	DWORD systime;
	double sec = 0;
	
	systime = timeGetTime();
	sec  = systime/1000.0;
#endif

	return(sec);
}

#ifndef WIN32
/*
	Set keyboard handling to be canonical or not
*/
void SetKeyboardIO(int on)
{
	static struct termios stored_settings;
   struct termios new_settings;

	if (on) {
   	tcgetattr(0,&stored_settings);
   	new_settings = stored_settings;

   	/* Disable canonical mode, and set buffer size to 1 byte */
   	new_settings.c_lflag &= (~ICANON);
   	new_settings.c_cc[VTIME] = 0;
   	new_settings.c_cc[VMIN] = 1;

   	tcsetattr(0,TCSANOW,&new_settings);
	} else {
		tcsetattr(0,TCSANOW,&stored_settings);
	}
}
#endif

/*
   Calculate the logarith with an arbitrary base
*/
double Logarithm(double n,double base)
{
   return(log(n) / log(base));
}

/*
   Return n1 to the power of n2
*/
long PowerInt(long n1,long n2)
{
   long i;
   long answer = 1;

   for (i=1;i<=n2;i++)
      answer *= n1;
   return(answer);
}

/*
   Integer square root
*/
long intsqrt(long l)
{
  long i,j=0,k;

  for (i=0x4000;i!=0;i>>=1) {
    k = (j + i) * (j + i);
    if (k <= l)
      j += i;
    if (k == l)
      break;
  }

  return(j);
}

/*
   Iterative definition of factorial
*/
double Factorial(int n)
{
   int i;
   double factor = 1;

   for (i=2;i<=n;i++)
      factor *= i;

   return(factor);
}

/*
   Series expansion for the moddified Bessel function of the 0 order
   This is correct to 16 decimal places
*/
double BesselI0(double x)
{
   double sum,xx,x2;

   sum = 1;
   x2 = x * x;
   xx = x2;

   sum += xx / 4.0;
   xx  *= x2;
   sum += xx / 64.0;
   xx *= x2;
   sum += xx / 2304.0;
   xx *= x2;
   sum += xx / 147456.0;
   xx *= x2;
   sum += xx / 14745600.0;
   xx *= x2;
   sum += xx / 2123366400.0;
   xx *= x2;
   sum += xx / 416179814400.0;
   xx *= x2;
   sum += xx / 106542032486400.0;

   return(sum);
}

double Bessel(double x,int n)
{
   double v,theta,dtheta,sum = 0;

   dtheta = PI / 10000;
   for (theta=0;theta<PI;theta+=dtheta) {
      v = cos(x * sin(theta) - n*theta);
      sum += v * dtheta;
   }

   return(sum/PI);
}

/* 
	n'th zero of the m'th bessel function 
*/
double BesselZero(int n,int m)
{
   int sign1,sign2;
   double x=0,y,ncount=0;
/*
double zeros[9][11] = {
   {-1.0,2.406, 5.520, 8.655,11.792,14.932,18.071,21.213,24.352,27.495,30.634},
   {-1.0,0.000, 3.832, 7.016,10.174,13.324,16.471,19.616,22.761,25.904,29.047},
   {-1.0,0.000, 5.136, 8.417,11.621,14.796,17.961,21.117,24.271,27.420,30.570},
   {-1.0,0.000, 6.381, 9.762,13.016,16.224,19.410,22.583,25.749,28.909,32.065},
   {-1.0,0.000, 7.589,11.065,14.374,17.616,20.828,24.019,27.200,30.371,33.538},
   {-1.0,0.000, 8.772,12.339,15.701,18.981,22.218,25.431,28.627,31.812,34.989},
   {-1.0,0.000, 9.937,13.589,17.005,20.321,23.587,26.820,30.035,33.233,36.423},
   {-1.0,0.000,11.087,14.822,18.288,21.642,24.935,28.192,31.423,34.638,37.839},
   {-1.0,0.000,12.226,16.038,19.556,22.945,26.268,29.545,32.797,36.025,39.242},
}; */

   if (n <= 0)
      return(0.0);

   y = Bessel(x,m);
   if (ABS(y) < 0.001) {
      if (n == 1)
         return(0.0);
      ncount++;
      x += 0.001;
      sign1 = SIGN(Bessel(x,m));
   } else {
      sign1 = SIGN(x);
   }
   x += 0.001;
   for (;x<100;x+=0.001) {
      y = SIGN(Bessel(x,m));
      sign2 = SIGN(y);
      if (sign2 != sign1) {
         sign1 = sign2;
         ncount++;
         if (ncount >= n)
            return(x);
      }
   }
   return(-1.0);
}

/*
	Natural log of the gamma function
	Derived from "Numerical Receipes in C"
	xx > 0
*/
double LnGamma(double xx)
{
	int j;
   double x,y,tmp,ser;
   double cof[6] = {
		76.18009172947146,    -86.50532032941677,
      24.01409824083091,    -1.231739572450155,
      0.1208650973866179e-2,-0.5395239384953e-5
	};

   y = x = xx;
   tmp = x + 5.5 - (x + 0.5) * log(x + 5.5);
   ser = 1.000000000190015;
   for (j=0;j<=5;j++) 
		ser += (cof[j] / ++y);
   return(log(2.5066282746310005 * ser / x) - tmp);
}

/* 2 x 2 matrices ----------------------------------------------- */

/*
	Normal mathematical convention is a_row_column
   | a1  b1 |      | a_0_0  a_0_1 |
   |        |  or  |              |
   | a2  b2 |      | a_1_0  a_1_1 |
*/
double Determinant22(double a[2][2])
{
	return(a[0][0]*a[1][1] - a[0][1]*a[1][0]);
}

double Det2x2(double a1,double b1,double a2,double b2)
{
  return(a1*b2 - b1*a2);
}

int Inverse22(double m[2][2],double i[2][2])
{
   double det;

   det = Determinant22(m);
   if (ABS(det) <= EPSILON)
      return(FALSE);

   i[0][0] =   m[1][1] / det;
   i[0][1] = - m[0][1] / det;
   i[1][0] = - m[1][0] / det;
   i[1][1] =   m[0][0] / det;
   return(TRUE);
}

/* 3 x 3 matrices ----------------------------------------------- */

/*
	| a1  b1  c1 |      | a_0_0  a_0_1  a_0_2 |
	|            |      |                     |
	| a2  b2  c2 |  or  | a_1_0  a_1_1  a_1_2 |
	|            |      |                     |
	| a3  b3  c3 |      | a_2_0  a_2_1  a_2_2 |
*/
double Determinant33(double a[3][3])
{
	return(a[0][0]*(a[1][1]*a[2][2] - a[1][2]*a[2][1]) +
          a[0][1]*(a[1][0]*a[2][2] - a[1][2]*a[2][0]) +
          a[0][2]*(a[1][0]*a[2][1] - a[2][0]*a[1][1]));
}

double Det3x3(
	double a1,double a2,double a3,
	double b1,double b2,double b3,
	double c1,double c2,double c3)
{
  return( a1 * Det2x2(b2,b3,c2,c3)
        - b1 * Det2x2(a2,a3,c2,c3)
        + c1 * Det2x2(a2,a3,b2,b3));
}

int Inverse33(double a[3][3],double b[3][3])
{
   double det;

   det = Determinant33(a);
   if (ABS(det) <= EPSILON)
      return(FALSE);

   b[0][0] =   (a[1][1]*a[2][2] - a[2][1]*a[1][2]) / det;
   b[0][1] = - (a[0][1]*a[2][2] - a[2][1]*a[0][2]) / det;
   b[0][2] =   (a[0][1]*a[1][2] - a[1][1]*a[0][2]) / det;

   b[1][0] = - (a[1][0]*a[2][2] - a[2][0]*a[1][2]) / det;
   b[1][1] =   (a[0][0]*a[2][2] - a[2][0]*a[0][2]) / det;
   b[1][2] = - (a[0][0]*a[1][2] - a[1][0]*a[0][2]) / det;

   b[2][0] =   (a[1][0]*a[2][1] - a[2][0]*a[1][1]) / det;
   b[2][1] = - (a[0][0]*a[2][1] - a[2][0]*a[0][1]) / det;
   b[2][2] =   (a[0][0]*a[1][1] - a[1][0]*a[0][1]) / det;
	
   return(TRUE);
}

/* 4 x 4 matrices ----------------------------------------------- */ 

/*
   | a1  b1  c1  d1 |      | a_0_0  a_0_1  a_0_2  a_0_3 |
   |                |      |                            |
   | a2  b2  c2  d2 |  or  | a_1_0  a_1_1  a_1_2  a_1_3 |
   |                |      |                            |
   | a3  b3  c3  d3 |      | a_2_0  a_2_1  a_2_2  a_2_3 |
   |                |      |                            |
   | a4  b4  c4  d4 |      | a_3_0  a_3_1  a_3_2  a_3_3 |
*/
double Determinant44(double a[4][4])
{
   double a1,a2,a3,a4,b1,b2,b3,b4,c1,c2,c3,c4,d1,d2,d3,d4;

   a1 = a[0][0]; b1 = a[0][1]; c1 = a[0][2]; d1 = a[0][3];
   a2 = a[1][0]; b2 = a[1][1]; c2 = a[1][2]; d2 = a[1][3];
   a3 = a[2][0]; b3 = a[2][1]; c3 = a[2][2]; d3 = a[2][3];
   a4 = a[3][0]; b4 = a[3][1]; c4 = a[3][2]; d4 = a[3][3];

   return(a1 * Det3x3(b2,b3,b4,c2,c3,c4,d2,d3,d4)
        - b1 * Det3x3(a2,a3,a4,c2,c3,c4,d2,d3,d4)
        + c1 * Det3x3(a2,a3,a4,b2,b3,b4,d2,d3,d4)
        - d1 * Det3x3(a2,a3,a4,b2,b3,b4,c2,c3,c4));
}

/*
   Calculate the adjoint of a 4x4 matrix
   Let  a_i_j  denote the minor determinant of matrix A obtained by
   deleting the ith row and jth column from A.
   Let b_i_j = (-1)^(i+j) a_j_i
   The matrix B = (b_i_j) is the adjoint of A
*/
void Adjoint44(double a[4][4],double b[4][4])
{
  double a1,a2,a3,a4,b1,b2,b3,b4;
  double c1,c2,c3,c4,d1,d2,d3,d4;

  a1 = a[0][0]; b1 = a[0][1]; c1 = a[0][2]; d1 = a[0][3];
  a2 = a[1][0]; b2 = a[1][1]; c2 = a[1][2]; d2 = a[1][3];
  a3 = a[2][0]; b3 = a[2][1]; c3 = a[2][2]; d3 = a[2][3];
  a4 = a[3][0]; b4 = a[3][1]; c4 = a[3][2]; d4 = a[3][3];

  b[0][0] =   Det3x3(b2,b3,b4,c2,c3,c4,d2,d3,d4);
  b[1][0] = - Det3x3(a2,a3,a4,c2,c3,c4,d2,d3,d4);
  b[2][0] =   Det3x3(a2,a3,a4,b2,b3,b4,d2,d3,d4);
  b[3][0] = - Det3x3(a2,a3,a4,b2,b3,b4,c2,c3,c4);

  b[0][1] = - Det3x3(b1,b3,b4,c1,c3,c4,d1,d3,d4);
  b[1][1] =   Det3x3(a1,a3,a4,c1,c3,c4,d1,d3,d4);
  b[2][1] = - Det3x3(a1,a3,a4,b1,b3,b4,d1,d3,d4);
  b[3][1] =   Det3x3(a1,a3,a4,b1,b3,b4,c1,c3,c4);

  b[0][2] =   Det3x3(b1,b2,b4,c1,c2,c4,d1,d2,d4);
  b[1][2] = - Det3x3(a1,a2,a4,c1,c2,c4,d1,d2,d4);
  b[2][2] =   Det3x3(a1,a2,a4,b1,b2,b4,d1,d2,d4);
  b[3][2] = - Det3x3(a1,a2,a4,b1,b2,b4,c1,c2,c4);

  b[0][3] = - Det3x3(b1,b2,b3,c1,c2,c3,d1,d2,d3);
  b[1][3] =   Det3x3(a1,a2,a3,c1,c2,c3,d1,d2,d3);
  b[2][3] = - Det3x3(a1,a2,a3,b1,b2,b3,d1,d2,d3);
  b[3][3] =   Det3x3(a1,a2,a3,b1,b2,b3,c1,c2,c3);
}

/*
	Calculate the inverse of a 4x4 matrix
    -1     1
   A  = ------- adjoint A
         det A
*/
int Inverse44(double a[4][4],double b[4][4])
{
   int i,j;
   double det;

   Adjoint44(a,b);
   det = Determinant44(a);

   if (ABS(det) < EPSILON) /* Singular */
		return(FALSE);

   for (i=0;i<4;i++)
      for (j=0;j<4;j++)
         b[i][j] = b[i][j] / det;

   return(TRUE);
}

void Identity44(double a[4][4])
{
   int i,j;

   for (i=0;i<4;i++)
      for (j=0;j<4;j++)
         a[i][j] = (double) (i == j);
}

/*
	c = a * b
*/
void Multiply44(double a[4][4],double b[4][4], double c[4][4])
{
	int i,j,k;

	for (i=0;i<4;i++) {
		for (j=0;j<4;j++) {
			c[i][j] = 0;
			for (k=0;k<4;k++) {
				c[i][j] += a[i][k] * b[k][j];
			}
		}
	}
}

/*
   Recursive definition of determinate using expansion by minors.
*/
double Determinant(double **a,int n)
{
   int i,j,j1,j2;
   double det = 0;
   double **m = NULL;

   if (n < 1) { /* Error */

   } else if (n == 1) { /* Shouldn't get used */
      det = a[0][0]