/***********************************************************************/
/* Open Visualization Data Explorer                                    */
/* (C) Copyright IBM Corp. 1989,1999                                   */
/* ALL RIGHTS RESERVED                                                 */
/* This code licensed under the                                        */
/*    "IBM PUBLIC LICENSE - Open Visualization Data Explorer"          */
/***********************************************************************/

#include <dxconfig.h>



#include <math.h>

/*
 * NB - not meant to be compiled directly
 * #define Type, c1,..., Neg3,... and then #include this file     
 */

#if INLINE
#    define STATIC static
#    if paxc
/*#        pragma inline (Neg,Add,Sub,Mul,Div,Min,Max,Lerp)*/
/*#        pragma inline (Dot,Cross,Length,Normalize)*/
/*#        pragma inline (Apply,Transpose,Invert,AdjointTranspose)*/
#    endif
#else
#    define STATIC
#endif

#ifdef Neg3
STATIC Point3 Neg3(Point3 v)
{
    v.c1 = -v.c1;
    v.c2 = -v.c2;
    v.c3 = -v.c3;
    return v;
}
#endif

#ifdef Add3
STATIC Point3 Add3(Point3 v, Point3 w)
{
    v.c1 += w.c1;
    v.c2 += w.c2;
    v.c3 += w.c3;
    return v;
}
#endif

#ifdef Sub3
STATIC Point3 Sub3(Point3 v, Point3 w)
{
    v.c1 -= w.c1;
    v.c2 -= w.c2;
    v.c3 -= w.c3;
    return v;
}
#endif

#ifdef Mul3
STATIC Point3 Mul3(Point3 v, double f)
{
    v.c1 *= f;
    v.c2 *= f;
    v.c3 *= f;
    return v;
}
#endif

#ifdef Div3
STATIC Point3 Div3(Point3 v, double f)
{
    v.c1 /= f;
    v.c2 /= f;
    v.c3 /= f;
    return v;
}
#endif

#ifdef Min3
STATIC Point3 Min3(Point3 v, Point3 w)
{
    v.c1 = v.c1<w.c1? v.c1 : w.c1;
    v.c2 = v.c2<w.c2? v.c2 : w.c2;
    v.c3 = v.c3<w.c3? v.c3 : w.c3;
    return v;
}
#endif

#ifdef Max3
STATIC Point3 Max3(Point3 v, Point3 w)
{
    v.c1 = v.c1>w.c1? v.c1 : w.c1;
    v.c2 = v.c2>w.c2? v.c2 : w.c2;
    v.c3 = v.c3>w.c3? v.c3 : w.c3;
    return v;
}
#endif

#ifdef Lerp3
STATIC Point3 Lerp3(double t, Point3 a, Point3 b)
{
    float s = 1-t;
    a.c1 = a.c1*t + b.c1*s;
    a.c2 = a.c2*t + b.c2*s;
    a.c3 = a.c3*t + b.c3*s;
    return a;
}
#endif

#ifdef Dot3
STATIC float Dot3(Point3 v, Point3 w)
{
    return v.x*w.x + v.y*w.y + v.z*w.z;
}
#endif

#ifdef Cross3
STATIC Point3 Cross3(Point3 v, Point3 w)
{
    Point3 r;
    r.c1 = v.c2*w.c3 - v.c3*w.c2;
    r.c2 = v.c3*w.c1 - v.c1*w.c3;
    r.c3 = v.c1*w.c2 - v.c2*w.c1;
    return r;
}
#endif

#ifdef Length3
STATIC double Length3(Point3 v)
{
    double v1 = v.c1;
    double v2 = v.c2;
    double v3 = v.c3;
    return sqrt(v1*v1 + v2*v2 + v3*v3);
}
#endif

#ifdef Normalize3
STATIC Point3 Normalize3(Point3 v)
{
    float d = Length3(v);
    if (d) {
	v.c1 /= d;
	v.c2 /= d;
	v.c3 /= d;
    }
    return v;
}
#endif

#ifdef Apply3
STATIC Point3 Apply3(Point3 p, Matrix3 t)
{
    Point3 q;
    q.c1 = t.A[0][0]*p.c1 + t.A[1][0]*p.c2 + t.A[2][0]*p.c3 + t.b[0];
    q.c2 = t.A[0][1]*p.c1 + t.A[1][1]*p.c2 + t.A[2][1]*p.c3 + t.b[1];
    q.c3 = t.A[0][2]*p.c1 + t.A[1][2]*p.c2 + t.A[2][2]*p.c3 + t.b[2];
    return q;
}
#endif

#ifdef Transpose3
STATIC Matrix3 Transpose3(Matrix3 t)
{
    Matrix3 s;
    s.A[0][0] = t.A[0][0];  s.A[0][1] = t.A[1][0];  s.A[0][2] = t.A[2][0];
    s.A[1][0] = t.A[0][1];  s.A[1][1] = t.A[1][1];  s.A[1][2] = t.A[2][1];
    s.A[2][0] = t.A[0][2];  s.A[2][1] = t.A[1][2];  s.A[2][2] = t.A[2][2];
    s.b[0] = s.b[1] = s.b[2] = 0;  /* XXX - ??? */
    return s;
}
#endif

#define COF(a,b,c,d) (t.A[a][c] * t.A[b][d] - t.A[a][d] * t.A[b][c])

#ifdef Invert3
STATIC Matrix3 Invert3(Matrix3 t)
{
    Matrix3 s;
    Point3 o;
    double inv;

    inv = t.A[0][0]*COF(1,2,1,2)+t.A[0][1]*COF(1,2,2,0)+t.A[0][2]*COF(1,2,0,1);
    inv = 1.0 / inv;

    s.A[0][0] = COF(1,2,1,2) * inv;
    s.A[1][0] = COF(1,2,2,0) * inv;
    s.A[2][0] = COF(1,2,0,1) * inv;
    s.A[0][1] = COF(2,0,1,2) * inv;
    s.A[1][1] = COF(2,0,2,0) * inv;
    s.A[2][1] = COF(2,0,0,1) * inv;
    s.A[0][2] = COF(0,1,1,2) * inv;
    s.A[1][2] = COF(0,1,2,0) * inv;
    s.A[2][2] = COF(0,1,0,1) * inv;

    s.b[0] = s.b[1] = s.b[2] = 0;
    *(Point3*)s.b = Apply3(Neg3(*(Point3*)t.b), s);
	
    return s;
}
#endif

#ifdef Determinant3
STATIC float Determinant3(Matrix3 t)
{
    return
	t.A[0][0]*COF(1,2,1,2)+t.A[0][1]*COF(1,2,2,0)+t.A[0][2]*COF(1,2,0,1);
}
#endif

#ifdef AdjointTranspose3
STATIC Matrix3 AdjointTranspose3(Matrix3 t)
{
    Matrix3 s;

    s.A[0][0] = COF(1,2,1,2);
    s.A[0][1] = COF(1,2,2,0);
    s.A[0][2] = COF(1,2,0,1);
    s.A[1][0] = COF(2,0,1,2);
    s.A[1][1] = COF(2,0,2,0);
    s.A[1][2] = COF(2,0,0,1);
    s.A[2][0] = COF(0,1,1,2);
    s.A[2][1] = COF(0,1,2,0);
    s.A[2][2] = COF(0,1,0,1);

    /* XXX - ??? - NB check users if you change this */
    s.b[0] = s.b[1] = s.b[2] = 0;

    return s;
}
#endif

#undef STATIC