/*
 *   supportive subroutines...
 */
#include <math.h>
#include <stdio.h>
#undef min
#include "rtd.h"
#include "extern.h"
#define time(x) 0

/* returns dot product of two vectors */

double  dot (v1, v2)
    struct vector  *v1,
                   *v2;
{
    double  r;
    r = v1->x * v2->x + v1->y * v2->y + v1->z * v2->z;
    return (r);
}


/* multiplies a vector by a scalar */

scamult (scal, vect)
    double  scal;
    struct vector  *vect;
{
    vect->x *= scal;
    vect->y *= scal;
    vect->z *= scal;
    vect->l *= scal;
    vect->xzl *= scal; 
}


/* loads a vector with given x,y,z values */

mv (x, y, z, _vp) 
    double  x,
            y,
            z;
    struct vector  *_vp;
{
    _vp->x = x;
    _vp->y = y;
    _vp->z = z;
}


/* vector subtraction: v1 = v2 - v3 */

sv (v1, v2, v3) 
    struct vector  *v1,
                   *v2,
                   *v3;
{
    v1->x = v2->x - v3->x;
    v1->y = v2->y - v3->y;
    v1->z = v2->z - v3->z;
}



/* vector addition: v1 = v2 + v3 */

av (v1, v2, v3) 
    struct vector  *v1,
                   *v2,
                   *v3;
{
    v1->x = v2->x + v3->x;
    v1->y = v2->y + v3->y;
    v1->z = v2->z + v3->z;
}


/*  loads a transform matrix so that when vec is multiplied by it,
    the result is (1,0,0). it does so with out deforming space */

mt (vec, trans) 
    struct vector  *vec;
    struct mat *trans;
{
    vecl (vec);
    vexzl (vec);

    if (vec->xzl == 0.0) {
	trans->x.x = 0.0;
	trans->x.y = 1.0;
	trans->x.z = 0.0;
	trans->y.x = -1.0;
	trans->y.y = 0.0;
	trans->y.z = 0.0;
	trans->z.x = 0.0;
	trans->z.y = 0.0;
	trans->z.z = 1.0;
    }
    else {
	trans->x.x = (vec->x) / (vec->l);
	trans->x.y = (vec->y) / (vec->l);
	trans->x.z = (vec->z) / (vec->l);
	trans->y.x = -(vec->x) * (vec->y) / ((vec->l) * (vec->xzl));
	trans->y.y = (vec->xzl) / (vec->l);
	trans->y.z = -(vec->z) * (vec->y) / ((vec->l) * (vec->xzl));
	trans->z.x = -(vec->z) / (vec->xzl);
	trans->z.y = 0;
	trans->z.z = (vec->x) / (vec->xzl);
    }
}

/* multiply vector by matrix: vc1 = m*vc2 */

mt_vec (vc1, m, vc2)
    struct vector  *vc1,
                   *vc2;
    struct mat *m;
{
    double  x,
            y,
            z;

    x = dot (&(m->x), vc2);
    y = dot (&(m->y), vc2);
    z = dot (&(m->z), vc2);
    mv (x, y, z, vc1);
}



/* will update a file containing information on 
how we are coming along... useful for long runs*/

update (x, xmin, xmax)
    float   x,
            xmin,
            xmax;
{
    FILE * fff;
    static int  t;
    static double   size,
                    xm;
    int     tt,
            ttt,
            hr,
            min,
            sec;


    fprintf (stderr, "%d\n", (int) x);
    return;


    if (x == xmin) {
/*	  fff = fopen ("time.update", "w"); */
	fff = stderr;
	fprintf (fff, "starting...\n");
/*	  fclose (fff); */
	t = time (0);
	xm = xmin;
	size = (xmax - xmin);
    }
    else {
	x -= xm;
	tt = time (0) - t;
	ttt = (double) tt * size / x;
	tt = (double) tt * (size - x) / x;
	hr = tt / 3600;
	min = (tt % 3600) / 60;
	sec = tt % 60;
/*	  fff = fopen ("time.update", "w"); */
	fff = stderr;
	fprintf (fff, "%5.1f%% done.\n", 100.0 * x / size);
	fprintf (fff, "%d:%2d:%2d left.\n", hr, min, sec);
	hr = ttt / 3600;
	min = (ttt % 3600) / 60;
	sec = ttt % 60;
	fprintf (fff, "%d:%2d:%2d  total estimate.\n", hr, min, sec);
/*	  fclose (fff); */
    }
}


/* loads a vectors length... saves time to do this only when needed */

vecl (v)
    struct vector *v;
{
    v->l = sqrt ((v->x) * (v->x) + (v->y) * (v->y) + (v->z) * (v->z));
}


/* loads a vectors xz component. needed for mt */

vexzl (v)
    struct vector *v;
{
    v->xzl = sqrt ((v->x) * (v->x) + (v->z) * (v->z));
}