/*
    HWT_3D_mm: Acoustical 3-D Huygens wavefront/ray tracing for moving media.
    An extension for moving media by N. Zabotin, O.A. Godin and L.Ye. Zabotina
    (Copyright (C) 2013 University of Colorado Boulder) of an algorithm
    developed by P. Sava and S. Fomel as a part of Madagascar open project
    (Copyright (C) 2007 Colorado School of Mines).

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details (file COPYING.txt).

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/
#include "alloc.h"
#include "_bool.h"
#include "axa.h"
#include "_defs.h"
#include "error.h"
#include "file.h"
#include "getpar.h"
#include "point.h"
#include "simtab.h"
#include "vector.h"

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "math.h"

#include "HWT_3D_mm.h"

#ifndef _HWT_3D_mm_h

typedef struct{
    double x,y,z;       /* location */
    float  v,mx,my,mz;   /* value    */
}pt3d_cur;

typedef struct{
    double  v;   /* velocity */
    double  mx;   /* current */
    double  my;   /* current */
    double  mz;   /* current */
}vm;

typedef struct{
    double x,y,z;       /* location */
}point3d;

typedef struct{
    point3d P0,P1;       /* location */
}ray3d;


typedef struct{
    point3d V0,V1,V2;       /* location */
}triangle3d;

bool lg_surf;
point3d So;

#endif

#define MISSING -1
#define SMALL_NUM  0.00000001  // anything that avoids division overflow

static sf_axa az,ax,ay;
static sf_axa at,ag,ah;
static float ***vv;
static float ***curx;
static float ***cury;
static float ***curz;
static float Asurf,Bsurf,Csurf,dsurf;

/*------------------------------------------------------------*/
void hwt3d_init(float*** vv_in    /* velocity */,
                float*** curx_in    /* current */,
                float*** cury_in,
                float*** curz_in,
                sf_axis      az_in    /* z axis   */,
		sf_axis      ax_in    /* x axis   */,
		sf_axis      ay_in    /* y axis   */,
		sf_axis      at_in    /* t axis   */,
		sf_axis      ag_in    /* g axis   */,
		sf_axis      ah_in    /* h axis   */,
float Asurf_in,float Bsurf_in, float Csurf_in,float dsurf_in)
/*< initialize hwt3d >*/
{
    az = sf_nod(az_in);
    ax = sf_nod(ax_in);
    ay = sf_nod(ay_in);

    at = sf_nod(at_in);
    ag = sf_nod(ag_in);
    ah = sf_nod(ah_in);

    vv = vv_in;
    curx = curx_in;   cury = cury_in;    curz = curz_in;

    Asurf=Asurf_in;   Bsurf=Bsurf_in;  Csurf=Csurf_in;  dsurf=dsurf_in;
}

/*------------------------------------------------------------*/
vm hwt3d_getv(pt3d_cur P)
/*< get velocity from 3-D cube by linear interpolation >*/
{
    double  z, x, y;
    double rz,rx,ry;
    int    jz,jx,jy;
    double fz,fx,fy;
    vm v;
    int    kz,kx,ky;

    x = SF_MIN(SF_MAX(ax.o,P.x),ax.o+(ax.n-2)*ax.d);
    y = SF_MIN(SF_MAX(ay.o,P.y),ay.o+(ay.n-2)*ay.d);
    z = SF_MIN(SF_MAX(az.o,P.z),az.o+(az.n-2)*az.d);

    rx = (x-ax.o)/ax.d;
    jx = (int)(rx);
    fx =       rx-jx;
    kx = jx+1;

    ry = (y-ay.o)/ay.d;
    jy = (int)(ry);
    fy =       ry-jy;
    ky = jy+1;

    rz = (z-az.o)/az.d;
    jz = (int)(rz);
    fz =       rz-jz;
    kz = jz+1;

    if(ax.n==1) {jx=0; kx=0; fx=0;}
    if(ay.n==1) {jy=0; ky=0; fy=0;}
    if(az.n==1) {jz=0; kz=0; fz=0;}

    v.v = vv[ jy][ jx][ jz] * (1-fz)*(1-fx)*(1-fy) +
          vv[ jy][ jx][ kz] * (  fz)*(1-fx)*(1-fy) +
	  vv[ ky][ jx][ jz] * (1-fz)*(1-fx)*(  fy) +
	  vv[ ky][ jx][ kz] * (  fz)*(1-fx)*(  fy) +
	  vv[ jy][ kx][ jz] * (1-fz)*(  fx)*(1-fy) +
	  vv[ jy][ kx][ kz] * (  fz)*(  fx)*(1-fy) +
	  vv[ ky][ kx][ jz] * (1-fz)*(  fx)*(  fy) +
	  vv[ ky][ kx][ kz] * (  fz)*(  fx)*(  fy);

   v.mx = curx[ jy][ jx][ jz] * (1-fz)*(1-fx)*(1-fy) +
          curx[ jy][ jx][ kz] * (  fz)*(1-fx)*(1-fy) +
	  curx[ ky][ jx][ jz] * (1-fz)*(1-fx)*(  fy) +
	  curx[ ky][ jx][ kz] * (  fz)*(1-fx)*(  fy) +
	  curx[ jy][ kx][ jz] * (1-fz)*(  fx)*(1-fy) +
	  curx[ jy][ kx][ kz] * (  fz)*(  fx)*(1-fy) +
	  curx[ ky][ kx][ jz] * (1-fz)*(  fx)*(  fy) +
	  curx[ ky][ kx][ kz] * (  fz)*(  fx)*(  fy);

   v.my = cury[ jy][ jx][ jz] * (1-fz)*(1-fx)*(1-fy) +
          cury[ jy][ jx][ kz] * (  fz)*(1-fx)*(1-fy) +
	  cury[ ky][ jx][ jz] * (1-fz)*(1-fx)*(  fy) +
	  cury[ ky][ jx][ kz] * (  fz)*(1-fx)*(  fy) +
	  cury[ jy][ kx][ jz] * (1-fz)*(  fx)*(1-fy) +
	  cury[ jy][ kx][ kz] * (  fz)*(  fx)*(1-fy) +
	  cury[ ky][ kx][ jz] * (1-fz)*(  fx)*(  fy) +
	  cury[ ky][ kx][ kz] * (  fz)*(  fx)*(  fy);

   v.mz = curz[ jy][ jx][ jz] * (1-fz)*(1-fx)*(1-fy) +
          curz[ jy][ jx][ kz] * (  fz)*(1-fx)*(1-fy) +
	  curz[ ky][ jx][ jz] * (1-fz)*(1-fx)*(  fy) +
	  curz[ ky][ jx][ kz] * (  fz)*(1-fx)*(  fy) +
	  curz[ jy][ kx][ jz] * (1-fz)*(  fx)*(1-fy) +
	  curz[ jy][ kx][ kz] * (  fz)*(  fx)*(1-fy) +
	  curz[ ky][ kx][ jz] * (1-fz)*(  fx)*(  fy) +
	  curz[ ky][ kx][ kz] * (  fz)*(  fx)*(  fy);
    return(v);
}

/*------------------------------------------------------------*/
pt3d_cur get_point(vm v,pt3d_cur G)
/*< find velocity >*/
{
  pt3d_cur P;
  P.x = G.x;   P.y=G.y;  P.z=G.z;
  P.v = v.v;   P.mx=v.mx;  P.my=v.my;  P.mz=v.mz;
  return(P);
}

/*------------------------------------------------------------*/
pt3d_cur hwt3d_raytr(
		 pt3d_cur Tm,
		 pt3d_cur To,
		 float scaleray,FILE *stream,bool refl)
/*< ray tracing >*/
{
    pt3d_cur Tp;
    pt3d_cur Gm,Gp,Hm,Hp;
    vc3d TmTo;
    vm     v;

    double a1,a2,a3;
    vc3d   v1,v2,v3,  v4,v5;
    vc3d   qq,uu,ww; /* unit vectors */
    float ro;

    vc3d vv,uc;

    ro = To.v * at.d; /* ray step */

// wavefront normal
    vv.dx=To.x-Tm.x-at.d*Tm.mx;    // vector n (wave front normal)
    vv.dy=To.y-Tm.y-at.d*Tm.my;
    vv.dz=To.z-Tm.z-at.d*Tm.mz;
    TmTo=nor3d(&vv); // normalized vector vv

    /* axes unit vectors */
    v1 = axa3d(1);
    v2 = axa3d(2);
    v3 = axa3d(3);

    /* ray angle with axes*/
    a1 = ang3d(&TmTo, &v1); a1 = SF_ABS(a1);
    a2 = ang3d(&TmTo, &v2); a2 = SF_ABS(a2);
    a3 = ang3d(&TmTo, &v3); a3 = SF_ABS(a3);

    /* select reference unit vector as "most orthogonal" to an incomming ray */
    if(      SF_ABS(a1-90) <= SF_ABS(a2-90) &&
	     SF_ABS(a1-90) <= SF_ABS(a3-90) )
	ww=v1;
    else if( SF_ABS(a2-90) <= SF_ABS(a1-90) &&
	     SF_ABS(a2-90) <= SF_ABS(a3-90) )
	ww=v2;
    else
	ww=v3;

    /* build orthogonal wavefront (Gm,Gp,Hm,Hp) */

/*------------------------------------------------------------*/
    /* find Gm, Gp */
    uu   = vcp3d(&TmTo,&ww);       /* uu = TmTo x qq */
    qq   = nor3d(&uu);             /* qq = uu / |uu| */
    uu   = scl3d(&qq,ro*scaleray); /* uu = qq * ro */

    Gm   = tip3d_cur(&To,&uu);         /* Gm at tip of +uu from To */
    qq   = scl3d(&uu,-1);
    Gp   = tip3d_cur(&To,&qq);         /* Gp at tip of -uu from To */

    v=hwt3d_getv(Gm);
    Gm=get_point(v,Gm);

    v=hwt3d_getv(Gp);
    Gp=get_point(v,Gp);

/*------------------------------------------------------------*/
    uu = vec3d_cur(&Gm,&Gp);
    ww = nor3d(&uu);
/*------------------------------------------------------------*/
    /* find Hm, Hp */
    uu   = vcp3d(&TmTo,&ww);
    qq   = nor3d(&uu);
    uu   = scl3d(&qq,ro*scaleray);

    Hm   = tip3d_cur(&To,&uu);
    qq   = scl3d(&uu,-1);
    Hp   = tip3d_cur(&To,&qq);

    v=hwt3d_getv(Hm);
    Hm=get_point(v,Hm);
    v=hwt3d_getv(Hp);
    Hp=get_point(v,Hp);
    lg_surf=false;
/*------------------------------------------------------------*/
    /* execute HWT step from Tm & (Gm,Hm,To,Gp,Hp) to Tp */
    Tp=hwt3d_step(Tm,To,Gm,Gp,Hm,Hp,refl);
    return(Tp);
}

/*------------------------------------------------------------*/

pt3d_cur hwt3d_step(
		pt3d_cur Tm,
		pt3d_cur To,
		pt3d_cur Gm,
		pt3d_cur Gp,
		pt3d_cur Hm,
		pt3d_cur Hp,
		bool refl)
/*< one HWT time step >*/
{
    double gdx,gdy,gdz,gdr;
    double hdx,hdy,hdz,hdr;
    double ddx,ddy,ddz;
    double ax,bx,ay,by,az,bz;
    double a,b,c,del;

    double dx,dy,dz;
    double xo,yo,zo,pp;

    pt3d_cur Sm,Sp,S;
    vc3d TmTo,ToSm,ToSp;
    double am,ap;
    float ro;
    vm v;
/*------------------------------------------------------------*/

    ro = To.v * at.d;       /* ray step */

    gdr =(Gp.v-Gm.v)*at.d;
    gdx = Gp.x-Gm.x+ro*(Gp.mx/Gp.v-Gm.mx/Gm.v)+To.mx/To.v*gdr;   //Pm.mx/Pm.v
    gdy = Gp.y-Gm.y+ro*(Gp.my/Gp.v-Gm.my/Gm.v)+To.my/To.v*gdr;
    gdz = Gp.z-Gm.z+ro*(Gp.mz/Gp.v-Gm.mz/Gm.v)+To.mz/To.v*gdr;

    hdr =(Hp.v-Hm.v)*at.d;
    hdx = Hp.x-Hm.x+ro*(Hp.mx/Hp.v-Hm.mx/Hm.v)+To.mx/To.v*hdr;
    hdy = Hp.y-Hm.y+ro*(Hp.my/Hp.v-Hm.my/Hm.v)+To.my/To.v*hdr;
    hdz = Hp.z-Hm.z+ro*(Hp.mz/Hp.v-Hm.mz/Hm.v)+To.mz/To.v*hdr;


    /* find largest dd (avoid division by 0) */
    ddz = gdy*hdx - gdx*hdy;
    ddx = gdz*hdy - gdy*hdz;
    ddy = gdx*hdz - gdz*hdx;

    if( SF_ABS(ddz) >=SF_ABS(ddx) &&
	SF_ABS(ddz) > SF_ABS(ddy)) {

	ax  = gdr*hdy - hdr*gdy;
	bx  = gdy*hdz - gdz*hdy;
	ay  = gdr*hdx - gdx*hdr;
	by  = gdx*hdz - gdz*hdx;

	a = ddz*ddz + bx*bx + by*by;
	b =         -(ax*bx + ay*by);           b *= ro;
	c =           ax*ax + ay*ay - ddz*ddz;  c *= ro*ro;

	del = SF_MAX(b*b - a*c , 0.);
	del = sqrtf(del);

	dz = (-b+del)/a;  Sp.z = To.z + To.mz/To.v*ro + dz;
	dx = ro*ax-dz*bx; Sp.x = To.x + To.mx/To.v*ro + dx/  ddz;
	dy = ro*ay-dz*by; Sp.y = To.y + To.my/To.v*ro + dy/(-ddz);

	dz = (-b-del)/a;  Sm.z = To.z + To.mz/To.v*ro + dz;
	dx = ro*ax-dz*bx; Sm.x = To.x + To.mx/To.v*ro + dx/  ddz;
	dy = ro*ay-dz*by; Sm.y = To.y + To.my/To.v*ro + dy/(-ddz);

    } else if( SF_ABS(ddx) >=SF_ABS(ddy) &&
	       SF_ABS(ddx) > SF_ABS(ddz)) {

	ay  = gdr*hdz - hdr*gdz;
	by  = gdz*hdx - gdx*hdz;
	az  = gdr*hdy - gdy*hdr;
	bz  = gdy*hdx - gdx*hdy;

	a = ddx*ddx + by*by + bz*bz;
	b =         -(ay*by + az*bz);           b *= ro;
	c =           ay*ay + az*az - ddx*ddx;  c *= ro*ro;

	del = SF_MAX(b*b - a*c , 0.);
	del = sqrtf(del);

	dx = (-b+del)/a;  Sp.x = To.x + To.mx/To.v*ro + dx;
	dy = ro*ay-dx*by; Sp.y = To.y + To.my/To.v*ro + dy/  ddx;
	dz = ro*az-dx*bz; Sp.z = To.z + To.mz/To.v*ro + dz/(-ddx);

	dx = (-b-del)/a;  Sm.x = To.x + To.mx/To.v*ro + dx;
	dy = ro*ay-dx*by; Sm.y = To.y + To.my/To.v*ro + dy/  ddx;
	dz = ro*az-dx*bz; Sm.z = To.z + To.mz/To.v*ro + dz/(-ddx);

    } else {

	az  = gdr*hdx - hdr*gdx;
	bz  = gdx*hdy - gdy*hdx;
	ax  = gdr*hdz - gdz*hdr;
	bx  = gdz*hdy - gdy*hdz;

	a = ddy*ddy + bz*bz + bx*bx;
	b =         -(az*bz + ax*bx);           b *= ro;
	c =           az*az + ax*ax - ddy*ddy;  c *= ro*ro;

	del = SF_MAX(b*b - a*c , 0.);
	del = sqrtf(del);

	dy = (-b+del)/a;  Sp.y = To.y + To.my/To.v*ro + dy;
	dz = ro*az-dy*bz; Sp.z = To.z + To.mz/To.v*ro + dz/  ddy;
	dx = ro*ax-dy*bx; Sp.x = To.x + To.mx/To.v*ro + dx/(-ddy);

	dy = (-b-del)/a;  Sm.y = To.y + To.my/To.v*ro + dy;
	dz = ro*az-dy*bz; Sm.z = To.z + To.mz/To.v*ro + dz/  ddy;
	dx = ro*ax-dy*bx; Sm.x = To.x + To.mx/To.v*ro + dx/(-ddy);
    }

    TmTo = vec3d_cur( &Tm, &To); /* incomming ray */
    ToSm = vec3d_cur( &To, &Sm); /* candidate ray */
    ToSp = vec3d_cur( &To, &Sp); /* candidate ray */

    /* angle between ray segments */
    am = ang3d( &TmTo, &ToSm);
    ap = ang3d( &TmTo, &ToSp);

    /* select candidate point that moves forward */
    if(am<ap) S=Sm;
    else      S=Sp;

/*   Reflections */
// horizontal surface
    if ( refl) {
    if ( (Asurf*To.x+Bsurf*To.y+Csurf*To.z+dsurf)*(Asurf*S.x+Bsurf*S.y+Csurf*S.z+dsurf)<0 ){
       lg_surf=true;
       So.x=S.x;  So.y=S.y;  So.z=S.z;  //save the point before reflection
       xo=S.x;
       yo=S.y;
       zo=S.z;
       pp=(Asurf*xo+Bsurf*yo+Csurf*zo+dsurf)/(Asurf*Asurf+Bsurf*Bsurf+Csurf*Csurf);
       S.x=xo-2.0*Asurf*pp;
       S.y=yo-2.0*Bsurf*pp;
       S.z=zo-2.0*Csurf*pp;
    }
	}

/* ----------------------------------------------------------------- */
    v=hwt3d_getv(S);
    S=get_point(v,S);
    return S;
}

/*----------------------------------------------------------------*/
void wavefront(FILE *stream, int ah, int ag, pt3d_cur **wp, float oh,float og,float dh,float dg,float tt)
/*< wavefront text file >*/
{
float g,h;
int ih,ig;

for( ih=0; ih<ah; ih++) {
for( ig=0; ig<ag; ig++) {
   g=og+ig*dg;
   h=oh+ih*dh;
   fprintf(stream,"%8.4f %8.4f %10.5f %10.5f %10.5f %9.5f\n",g,h,wp[ih][ig].x,wp[ih][ig].y,wp[ih][ig].z,tt);
}}
}

/*----------------------------------------------------------------*/
point3d zinterp( pt3d_cur p1, point3d p0 )
/*< interpolation for reflected points to find x,y if z=0 >*/
{
point3d p;
p.z=0;
p.x=(p0.x*p1.z-p1.x*p0.z)/(p1.z-p0.z);
p.y=(p0.y*p1.z-p1.y*p0.z)/(p1.z-p0.z);
return p;
}

/*----------------------------------------------------------------*/
point3d yinterp( pt3d_cur p1, pt3d_cur p0, float Dpy )
/*< interpolatio for reflected points to find x,y if z=0 >*/
{
point3d p;
p.y=-Dpy;
p.x=p0.x*(p.y-p1.y)/(p0.y-p1.y)+p1.x*(p0.y-p.y)/(p0.y-p1.y);
p.z=p0.z*(p.y-p1.y)/(p0.y-p1.y)+p1.z*(p0.y-p.y)/(p0.y-p1.y);
return p;
}

/*--------------------------------------------------------------*/
pt3d_cur* pt3d_curalloc1( size_t n1)
/*< alloc point3d 1-D vector >*/
{
    pt3d_cur *ptr;
    ptr =  (pt3d_cur*) sf_alloc(n1,sizeof(pt3d_cur));
    return ptr;
}

/*--------------------------------------------------------------*/
pt3d_cur** pt3d_curalloc2( size_t n1,
		   size_t n2)
/*< alloc point3d 2-D vector >*/
{
    size_t i2;
    pt3d_cur **ptr;

    ptr = (pt3d_cur**) sf_alloc(n2,sizeof(pt3d_cur*));
    ptr[0] = pt3d_curalloc1(n1*n2);
    for (i2=1; i2<n2; i2++) {
	ptr[i2] = ptr[0]+i2*n1;
    }
    return ptr;
}

/*------------------------------------------------------------*/
vc3d vec3d_cur(pt3d_cur* O, pt3d_cur* 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 vec3d_uc(pt3d_cur* A)
/*< build 3D vector u/c >*/
{
    vc3d V;
    V.dx = A->mx/A->v;
    V.dy = A->my/A->v;
    V.dz = A->mz/A->v;
    return V;
}

/*-----------------------------------------------------------*/
pt3d_cur tip3d_cur(pt3d_cur* O, vc3d* V)
/*< tip of a 3D vector >*/
{
    pt3d_cur A;
    A.x = O->x + V->dx;
    A.y = O->y + V->dy;
    A.z = O->z + V->dz;
    A.v = 0;
    return A;
}


// 3D Time front
/*------------------------------------------------------------*/
point3d get_point3d(pt3d_cur G)
/*< transform triangle coordinates from pt3d_cur to point3d structure >*/
{
  point3d P;
  P.x = G.x;   P.y=G.y;  P.z=G.z;
  return(P);
}

/*------------------------------------------------------------*/
point3d crossp3d(point3d* U, point3d* V)
/*< vector product of 3D vectors - cross product >*/
{
    point3d W;
    W.x=(U->y*V->z) - (V->y*U->z);
    W.y=(U->z*V->x) - (V->z*U->x);
    W.z=(U->x*V->y) - (V->x*U->y);
    return W;
}

/*------------------------------------------------------------*/
point3d mult3d(point3d* U, float a)
/*< vector product of 3D vectors - cross product >*/
{
    point3d W;
    W.x=a*U->x;
    W.y=a*U->y;
    W.z=a*U->z;
    return W;
}

/*------------------------------------------------------------*/
point3d zeropoint(double x)
/*< fill 3d point >*/
{
    point3d V;
    V.x = x;
    V.y = x;
    V.z = x;
    return V;
}

/*------------------------------------------------------------*/
point3d dif3d(point3d* A, point3d* B)
/*< build 3D vector >*/
{
    point3d V;
    V.x = A->x - B->x;
    V.y = A->y - B->y;
    V.z = A->z - B->z;
    return V;
}

/*------------------------------------------------------------*/
point3d sum3d(point3d* A, point3d* B)
/*< build 3D vector >*/
{
    point3d V;
    V.x = A->x + B->x;
    V.y = A->y + B->y;
    V.z = A->z + B->z;
    return V;
}

/*------------------------------------------------------------*/
float dot3d(point3d* U, point3d* V)
/*< dot product (3D) which allows vector operations in arguments >*/
{
    return U->x*V->x + U->y*V->y + U->z*V->z;
}

/*------------------------------------------------------------*/
int intersect_RayTriangle( ray3d R, triangle3d T, point3d* Intsec)
/*<  intersect a straight line with a 3D triangle >*/
{
//    Input:  a straight line R, and a triangle T
//    Output: *I = intersection point (when it exists)
//    Return: -1 = triangle is degenerate (a segment or point)
//             0 = disjoint (no intersect)
//             1 = intersect in unique point Intsec
//             2 = are in the same plane

    point3d    u, v, n;             // triangle vectors
    point3d    dir, w0, w;          // straight line vectors
    float   r, a, b;             // params to calc straight line-plane intersect
    float    uu, uv, vv, wu, wv, D;
    float s, t;

    // get triangle edge vectors and plane normal
    u = dif3d(&T.V1,&T.V0);
    v = dif3d(&T.V2,&T.V0);

    n = crossp3d(&u,&v);             // cross product
    if (n.x == 0.0 && n.y==0.0 && n.z==0.0)            // triangle is degenerate
        return -1;                 // do not deal with this case

    dir = dif3d(&R.P1,&R.P0);       // straight line direction vector
    w0 = dif3d(&R.P0,&T.V0);
    a = -dot3d(&n,&w0);
    b = dot3d(&n,&dir);

    if (fabs(b) < SMALL_NUM) {     // straight line is parallel to triangle plane
        if (a == 0)                // straight line lies in triangle plane
            return 2;
        else return 0;             // straight line disjoint from plane
    }

    // get intersect point of the straight line with triangle plane
    r = a / b;
    if (r < 0.0)                   // straight line goes away from triangle
        return 0;                  // => no intersect
    // for a segment, also test if (r > 1.0) => no intersect

    w=mult3d(&dir,r);
    *Intsec = sum3d(&w,&R.P0);          // intersect point of straight line and plane

    // is I inside T?
    uu = dot3d(&u,&u);
    uv = dot3d(&u,&v);
    vv = dot3d(&v,&v);

    w=dif3d(Intsec,&T.V0);
    wu = dot3d(&w,&u);
    wv = dot3d(&w,&v);

    D = uv * uv - uu * vv;
    // get and test parametric coords
    s = (uv * wv - vv * wu) / D;
    if (s < 0.0 || s > 1.0)        // I is outside T
        return 0;
    t = (uv * wu - uu * wv) / D;
    if (t < 0.0 || (s + t) > 1.0)  // I is outside T
        return 0;

    return 1;                      // I is in T
}