#include <math.h>
#include "vector.h"
#include "alloc.h"

#include "_defs.h"
/*^*/

#include "file.h"
/*^*/

#include "point.h"
/*^*/

#ifndef _sf_vector_h

typedef struct{
    double dx,dz;
}vc2d;
/*^*/

typedef struct{
    double dx,dy,dz;
}vc3d;
/*^*/

#endif

/*------------------------------------------------------------*/
vc3d vec3d(pt3d* O, pt3d* A)
/*< build 3D vector >*/
{
    vc3d V;

    V.dx = A->x - O->x;
    V.dy = A->y - O->y;
    V.dz = A->z - O->z;

    return V;
}

/*------------------------------------------------------------*/
vc3d axa3d( int n)
/*< build 3D unit vector >*/
{
    vc3d V;
    V.dx=V.dy=V.dz=0.;

    if(n<1 || n>3) n=1;
    
    switch(n) {
	case 3:
	    V.dy=1;
	    break;		
	case 2:
	    V.dx=1;
	    break;	    
	case 1:
	    V.dz=1;
	    break;
    }
    return V;
}

/*------------------------------------------------------------*/
double scp3d(vc3d* U, vc3d* V)
/*< scalar product of 3D vectors >*/
{
    return U->dx*V->dx + U->dy*V->dy + U->dz*V->dz;
}

/*------------------------------------------------------------*/
vc3d vcp3d(vc3d* U, vc3d* V)
/*< vector product of 3D vectors >*/
{
    vc3d W;

    W.dx=(U->dy*V->dz) - (V->dy*U->dz);
    W.dy=(U->dz*V->dx) - (V->dz*U->dx);
    W.dz=(U->dx*V->dy) - (V->dx*U->dy);

    return W;
}

/*------------------------------------------------------------*/
double len3d(vc3d* V)
/*< 3D vector length >*/
{
    double l;
    l = sqrtf( V->dx*V->dx + 
	       V->dy*V->dy + 
	       V->dz*V->dz);
    return l;
}

/*------------------------------------------------------------*/
vc3d nor3d(vc3d* V)
/*< normalize 3D vector >*/
{
    vc3d W;
    double l;
    l = len3d(V);

    W.dx = V->dx / l;
    W.dy = V->dy / l;
    W.dz = V->dz / l;

    return W;
}

/*------------------------------------------------------------*/
double ang3d(vc3d* U, vc3d* V)
/*< angle between 3D vectors >*/
{
    double c,a;

    c = U->dx * V->dx +
	U->dy * V->dy +
	U->dz * V->dz;
    c/= len3d(U);
    c/= len3d(V);

    c = SF_SIG(c) * SF_MIN( 1., SF_ABS(c));

    a = 180*acosf(c)/SF_PI;

    return a;
}


/*------------------------------------------------------------*/
vc3d scl3d(vc3d* V, float s)
/*< scale a 3D vector >*/
{
     vc3d W;

     W.dx = V->dx * s;
     W.dy = V->dy * s;
     W.dz = V->dz * s;

     return W;
}