/* Output from p2c, the Pascal-to-C translator */
/* From input file "math.pas" */


#include <p2c/p2c.h>


#define MATH_G
#include "math.h"


Void vector_init(double i, double j, double k, vector *result)
{
  result->U1.i = i;
  result->U1.j = j;
  result->U1.k = k;
}


boolean vector_null_check(vector *v)
{
  boolean Result;

  Result = false;
  if (v->U1.i != 0.0)
    return Result;
  if (v->U1.j == 0.0) {
    if (v->U1.k == 0.0)
      return true;
  }
  return Result;
}


boolean vector_equal_check(vector *v1, vector *v2)
{
  boolean Result;

  Result = false;
  if (v1->U1.i != v2->U1.i)
    return Result;
  if (v1->U1.j == v2->U1.j) {
    if (v1->U1.k == v2->U1.k)
      return true;
  }
  return Result;
}


Void vector_negate(vector *v, vector *result)
{
  result->U1.i = -v->U1.i;
  result->U1.j = -v->U1.j;
  result->U1.k = -v->U1.k;
}


double vector_mag(vector *v)
{
  return sqrt(v->U1.i * v->U1.i + v->U1.j * v->U1.j + v->U1.k * v->U1.k);
}


Void vector_unit(vector *v, vector *result)
{
  double vmag;

  vmag = vector_mag(v);
  if (vmag == 0.0) {
    printf("ERROR: Attempt to unitize a NULL vector. Program aborted.\n");
    _Escape(0);
  }
  result->U1.i = v->U1.i / vmag;
  result->U1.j = v->U1.j / vmag;
  result->U1.k = v->U1.k / vmag;
}


Void vector_cross(vector *v1, vector *v2, vector *result)
{
  result->U1.i = v1->U1.j * v2->U1.k - v1->U1.k * v2->U1.j;
  result->U1.j = v1->U1.k * v2->U1.i - v1->U1.i * v2->U1.k;
  result->U1.k = v1->U1.i * v2->U1.j - v1->U1.j * v2->U1.i;
}


double vector_dot(vector *v1, vector *v2)
{
  return (v1->U1.i * v2->U1.i + v1->U1.j * v2->U1.j + v1->U1.k * v2->U1.k);
}


Void vector_add(vector *v1, vector *v2, vector *result)
{
  result->U1.i = v1->U1.i + v2->U1.i;
  result->U1.j = v1->U1.j + v2->U1.j;
  result->U1.k = v1->U1.k + v2->U1.k;
}


Void vector_subtract(vector *v1, vector *v2, vector *result)
{
  result->U1.i = v1->U1.i - v2->U1.i;
  result->U1.j = v1->U1.j - v2->U1.j;
  result->U1.k = v1->U1.k - v2->U1.k;
}


Void dir_cosines_init(vector *v1, vector *v2, vector *v3, dir_cosines *dc)
{
  dc->U1.row1 = *v1;
  dc->U1.row2 = *v2;
  dc->U1.row3 = *v3;
}


Void dir_cosines_calc(vector *ivec, vector *up, dir_cosines *dc)
{
  vector normal_ivec, normal_up;

  vector_unit(ivec, &normal_ivec);
  dc->U1.row1 = normal_ivec;
  vector_unit(up, &normal_up);
  vector_cross(&normal_up, &normal_ivec, &dc->U1.row2);
  vector_unit(&dc->U1.row2, &dc->U1.row2);
  vector_cross(&normal_ivec, &dc->U1.row2, &dc->U1.row3);
  dir_cosines_invert(dc, dc);
}


Void dir_cosines_invert(dir_cosines *dc, dir_cosines *result)
{
  *result = *dc;
  swap_(result->m[1], &result->m[0][1]);
  swap_(result->m[2], &result->m[0][2]);
  swap_(&result->m[2][1], &result->m[1][2]);
}


Void coordinate_transformation(vector *v, dir_cosines *dc, vector *result)
{
  long i, j;

  for (i = 0; i <= 2; i++) {
    result->a[i] = 0.0;
    for (j = 0; j <= 2; j++)
      result->a[i] += v->a[j] * dc->m[j][i];
  }
}


double modulo(double x, double y)
{
  double temp, DUMMY;

  temp = x / y;
  return (modf(temp, &DUMMY) * y);
}


double pos_modulo(double x, double y)
{
  double temp, DUMMY;

  temp = x / fabs(y);
  if (temp < 0.0)
    return ((1 + modf(temp, &DUMMY)) * fabs(y));
  else
    return (modf(temp, &DUMMY) * fabs(y));
}


double power(double x, double y)
{
  double Result, t, DUMMY;

  if (x == 0)
    return 0.0;
  if (x > 0 && y >= 0) {
    t = y * log(x);
    return exp(t);
  }
  if (x > 0 && y < 0) {
    t = -y * log(x);
    return (1 / exp(t));
  }
  if (x >= 0)
    return Result;
  if (modf(y, &DUMMY) == 0) {
    if (y >= 0) {
      t = y * log(-x);
      if (modulo(y, 2.0) == 0)
	return exp(t);
      else
	return (-exp(t));
    } else {
      t = -y * log(-x);
      if (modulo(-y, 2.0) == 0)
	return (1 / exp(t));
      else
	return (-1 / exp(t));
    }
  } else {
    t = sqrt(-2.0);
    return Result;
  }
}


double tan_(double x)
{
  return (sin(x) / cos(x));
}


double arcsin_(double x)
{
  x = bounds(-1.0, x, 1.0);
  return (2 * atan(x / (1 + sqrt(1 - x * x))));   /*ken's method...neat*/
}


double arccos_(double x)
{
  return (M_PI / 2 - arcsin_(x));
}


double arctan2(double sinx, double cosx)
{
  double temp;

  sinx = bounds(-1.0, sinx, 1.0);
  cosx = bounds(-1.0, cosx, 1.0);
  if (cosx == 0.0) {
    if (sinx == 0)
      return 0.0;
    else {
      if (sinx > 0)
	return (0.5 * M_PI);
      else
	return (1.5 * M_PI);
    }
  } else {
    temp = atan(sinx / cosx);
    if (cosx < 0.0)
      temp += M_PI;
    else {
      if (sinx < 0.0)
	temp += 2 * M_PI;
    }
    return temp;
  }
}


double sign(double a)
{
  if (a > 0.0)
    return 1.0;
  else {
    if (a < 0.0)
      return -1.0;
    else
      return 0.0;
  }
}


double bounds(double low, double value, double high)
{
  if (value < low || high < value) {
    if (value < low)
      return low;
    else
      return high;
  } else
    return value;
}


long ibounds(int low, int value, int high)
{
  if (value < low || high < value) {
    if (value < low)
      return low;
    else
      return high;
  } else
    return value;
}


double max(double v1, double v2)
{
  if (v1 >= v2)
    return v1;
  else
    return v2;
}


double min(double v1, double v2)
{
  if (v1 <= v2)
    return v1;
  else
    return v2;
}


long imax(int v1, int v2)
{
  if (v1 >= v2)
    return v1;
  else
    return v2;
}


long imin(int v1, int v2)
{
  if (v1 <= v2)
    return v1;
  else
    return v2;
}


Void swap_(double *a, double *b)
{
  double temp;

  temp = *a;
  *a = *b;
  *b = temp;
}


void _math_init()
{
  /*nothing to initialize*/
  static int _was_initialized = 0;
  if (_was_initialized++)
    return;
}
/* p2c: Note: Remember to call _math_init() in main program [215] */



/* End. */