/* $Header: /server/users/students/chow/mci/netrestoreW/sparcsrc/RCS/rmath.c,v 1.1 1995/07/09 16:50:41 chow Exp chow $ */
/* Output from p2c, the Pascal-to-C translator */
/* From input file "math.pas" */


#include "p2c/p2c.h"


#define MATH_G
#include "rmath.h"
#include <math.h>


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


boolean vector_null_check(v)
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(v1, v2)
vector *v1, *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(v, result)
vector *v, *result;
{
  result->U1.i = -v->U1.i;
  result->U1.j = -v->U1.j;
  result->U1.k = -v->U1.k;
}


double vector_mag(v)
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(v, result)
vector *v, *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(v1, v2, result)
vector *v1, *v2, *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(v1, v2)
vector *v1, *v2;
{
  return (v1->U1.i * v2->U1.i + v1->U1.j * v2->U1.j + v1->U1.k * v2->U1.k);
}


Void vector_add(v1, v2, result)
vector *v1, *v2, *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(v1, v2, result)
vector *v1, *v2, *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(v1, v2, v3, dc)
vector *v1, *v2, *v3;
dir_cosines *dc;
{
  dc->U1.row1 = *v1;
  dc->U1.row2 = *v2;
  dc->U1.row3 = *v3;
}


Void dir_cosines_calc(ivec, up, dc)
vector *ivec, *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(dc, result)
dir_cosines *dc, *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(v, dc, result)
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(x, y)
double x, y;
{
  double temp, DUMMY;

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


double pos_modulo(x, y)
double x, 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(x, y)
double x, 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_(x)
double x;
{
  return (sin(x) / cos(x));
}


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


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


double arctan2(sinx, cosx)
double sinx, 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(a)
double a;
{
  if (a > 0.0)
    return 1.0;
  else {
    if (a < 0.0)
      return -1.0;
    else
      return 0.0;
  }
}


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


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


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


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


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


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


Void swap_(a, b)
double *a, *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. */