/*
 * mountain.c
 *
 * Peter Janssen
 *
 */

#include "geom.h"
#include "mountain.h"
#include "allocmatrix.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>


static void MountainEvolve();
static void MountainSetupLevelList();
static void MountainSetupCache();
static void MountainSetupCalcList();
int MountainBoundsIntersect();
int CheckGrid();
int IntersectTriangle();

#define    x2grid(m, x) (((x) - m->D2Bounds[0][0]) * m->InvXSize)
#define    y2grid(m, y) (((y) - m->D2Bounds[0][1]) * m->InvYSize)

#define    grid2x(m, x) ((x) * m->XSize + m->D2Bounds[0][0])
#define    grid2y(m, y) ((y) * m->YSize + m->D2Bounds[0][1])

#define OutOf2DBounds(p, b)                                 \
               ((p).x < (b)[0][0] || (p).x > (b)[1][0] ||   \
                (p).y < (b)[0][1] || (p).y > (b)[1][1] )

static Methods *iMountainMethods = NULL; 
static char mountainName[] = "mountain";
unsigned long MountainTests, MountainHits;
static int WhichCache[3][3] = { {TRUE, TRUE, FALSE},
                                {FALSE, FALSE, FALSE},
                                {FALSE, TRUE, TRUE}};


int MyPower(a, b)
int   a, b;
{
int   result, i;

   for (result = 1, i = 1; i <= b; i++) result *= a;
   return result;
}

/* Random Generator */

#define rand_m  (unsigned long)2147483647
#define rand_q  (unsigned long)127773

#define rand_a (unsigned int)16807
#define rand_r (unsigned int)2836

/*
 * F(z) = (az)%m
 *      = az-m(az/m)
 *
 * F(z)  = G(z)+mT(z)
 * G(z)  = a(z%q)- r(z/q)
 * T(z)  = (z/q) - (az/m)
 *
 * F(z)  = a(z%q)- rz/q+ m((z/q) - a(z/m))
 *       = a(z%q)- rz/q+ m(z/q) - az
*/

static long rand_seed;

static long rand_number()
{
long    lo, hi, test;

    hi   = rand_seed/rand_q;
    lo   = rand_seed%rand_q;

    test = rand_a*lo - rand_r*hi;

    if (test > 0)
        rand_seed = test;
    else
        rand_seed = test + rand_m;

    return rand_seed;
}

double GAUSS(s, x, y, level)
int  s, x, y, level;
{
   long    x1, x2;
   double  r, v1, v2;
   static double inv_rand_m;
   inv_rand_m = 1.0 / (double) rand_m;

   rand_seed = s + (y<<15) + x + (level<<8);
   do {
       x1 = rand_number();
       x2 = rand_number();
       v1 = (((double)x1 * inv_rand_m) * 2) - 1;
       v2 = (((double)x2 * inv_rand_m) * 2) - 1;
       r = v1*v1 + v2*v2;
   } while (r >= 1.0);
   r = sqrt(-2.0 * log(r) / r);
   return v1*r;
}


/* 
 * Create MountainObject and return Reference to it
 */
Mountain *MountainCreate(p1, p2, number, dimension, evolvelevel,
                         iterationlevel, scale, startgrid)
Vec2d            p1, p2;
int              number, evolvelevel, iterationlevel;
Float            dimension, scale;
StartGridStruct *startgrid;
{
        Mountain *mountain;
        Vec2d     size;
        int       gridmax;

        mountain = (Mountain *)share_malloc(sizeof(Mountain));
        size.u = fabs(p1.u - p2.u);
        size.v = fabs(p1.v - p2.v);
        if (equal(size.u, 0.) || equal(size.v, 0.)) {
            RLerror(RL_ADVISE, "Degenerated startgrid.\n");
            return (Mountain *) NULL;
        }
        mountain->D2Bounds[LOW][X] = min(p1.u, p2.u);
        mountain->D2Bounds[HIGH][X] = max(p1.u, p2.u);
        mountain->D2Bounds[LOW][Y] = min(p1.v, p2.v);
        mountain->D2Bounds[HIGH][Y] = max(p1.v, p2.v);

        if (number < 0 || number > 65565) {
            RLerror(RL_WARN, "Illegal MountainNumber : %d", number);
            number = 1000;
        }
        if (dimension < 0.5 || dimension > 3.5) {
            RLerror(RL_WARN, "Illegal FractalDimension : %f", dimension);
            dimension = 2.5;
        }
        if (iterationlevel < 0) {
            RLerror(RL_WARN, "Illegal IterationLevel : %d", iterationlevel);
            iterationlevel = 2;
        }
        if (evolvelevel < 0 || evolvelevel > iterationlevel) {
            RLerror(RL_WARN, "Illegal EvolveLevel : %d", evolvelevel);
            evolvelevel = 0;
        }
        if (scale < 0) {
            RLerror(RL_WARN, "Illegal ScaleFactor : %f", scale);
            scale = 1.0;
        }

        gridmax  = MyPower(2, iterationlevel) * (startgrid->Size - 1) + 1;
        mountain->XSize    = size.u / gridmax;
        mountain->YSize    = size.v / gridmax;
        mountain->InvXSize = 1.0 / mountain->XSize;
        mountain->InvYSize = 1.0 / mountain->YSize;

        mountain->Number           = number;
        mountain->IterationLevel   = iterationlevel;
        mountain->EvolveLevel      = evolvelevel;
        mountain->FractalDimension = dimension;
        mountain->Scale            = scale;
        mountain->StartGridSize    = startgrid->Size;
        mountain->StartGridData    = startgrid->Data;
        mountain->GridMax          = gridmax;


        MountainEvolve(mountain);
        MountainSetupLevelList(mountain);
        MountainSetupCache(mountain);
        MountainSetupCalcList(mountain);

        return mountain;
}

/*
 * Calculating MountainData up till iterationlevel 'EvolveLevel'
 */
static void MountainEvolve(m)
Mountain *m;
{
   Float      Scale;
   Float    **OldMountain, **NewMountain;
   int        OldSize, NewSize;
   int        XFar, XShort, YFar, YShort;
   int        i, x, y;

   Scale         = m->Scale;
   OldMountain   = m->StartGridData;
   OldSize       = m->StartGridSize;

   for (i = 1; i <= m->EvolveLevel; i++) {
       NewSize = 2 * OldSize - 1;
       Scale = Scale * pow((OldSize / (Float) NewSize), m->FractalDimension);
       NewMountain = AllocMatrix(Float, NewSize+1, NewSize+1);
       for (x = 0; x <= NewSize; x++)
           for (y = 0; y <= NewSize; y++) {
               XShort = (x+1)>>1;
               XFar   = (x|1) - XShort;
               YShort = (y+1)>>1;
               YFar   = (y|1) - YShort;
               NewMountain[x][y] =
                   (OldMountain[XFar][YFar] +
                    3 * (OldMountain[XFar][YShort] +
                         OldMountain[XShort][YFar] +
                         3 * OldMountain[XShort][YShort])) * 0.0625 +
                    GAUSS(m->Number, x, y, i) * Scale;
           }
       FreeMatrix(OldMountain, OldSize+1);
       OldMountain = NewMountain; OldSize = NewSize;
   }
   m->EvolvedData = OldMountain;
   m->EvolvedSize = OldSize;
}

/*
 * Calculate how many cells on every iterationlevel
 */
static void MountainSetupLevelList(m)
Mountain  *m;
{
int i, Ilevel, Elevel, Dlevel;

    Ilevel = m->IterationLevel;
    Elevel = m->EvolveLevel;

    Dlevel = Ilevel - Elevel + 1;
    m->MaxSize = (int *) share_malloc(Dlevel * sizeof(int));
    m->MaxSize[0] = MyPower(2, Elevel) * (m->StartGridSize - 1) +  1;
    for (i = 1; i < Dlevel; i++)
         m->MaxSize[i] = 2 * m->MaxSize[i - 1] - 1;
}

/*
 * Allocate memory for cache on every iterationlevel
 */
static void MountainSetupCache(m)
Mountain *m;
{
int   i,j, Ilevel, Elevel, Dlevel;

    Ilevel = m->IterationLevel;
    Elevel = m->EvolveLevel;

    if (Dlevel = Ilevel - Elevel) {
        m->Cache1 = (HeightCache **) share_malloc(Dlevel * sizeof(HeightCache *));
        m->Cache2 = (HeightCache **) share_malloc(Dlevel * sizeof(HeightCache *));
        for (i = 0; i < Dlevel; i++)  {
             m->Cache1[i] = (HeightCache *) share_malloc((2 * m->MaxSize[i + 1] + 1) * sizeof(HeightCache));
             m->Cache2[i] = (HeightCache *) share_malloc((2 * m->MaxSize[i + 1] + 1) * sizeof(HeightCache));
             for (j = 0; j <= 2 * m->MaxSize[i + 1]; j++) {
                  m->Cache1[i][j].x = -1;
                  m->Cache2[i][j].x = -1;
             }
        }
    }
}


/* 
 * Setup list to keep track of which cells are needed
 * on the different iterationlevels
 */
static void MountainSetupCalcList(m)
Mountain   *m;
{
int       i, j;
CalcList  *Chasingcl, *Cl;

    Chasingcl   = (CalcList *) share_malloc(sizeof(CalcList));
    for (i = 0; i <= 2; i++)
         for (j = 0; j <= 2; j++)
              Chasingcl->ToDo[i][j] = FALSE;
    m->ToCalcList       = Chasingcl;
    Chasingcl->Previous = NULL;
    for (i = m->IterationLevel; i > m->EvolveLevel + 1; i--) {
        Cl              = (CalcList *) share_malloc(sizeof(CalcList));
        Chasingcl->Next = Cl;
        Cl->Previous    = Chasingcl;
        Chasingcl       = Cl;
    }
    Chasingcl->Next = NULL;
}


/*
 * Intersect ray with mountain.
 */
int MountainIntersect(mountain, ray, mindist, maxdist)
Mountain   *mountain;
Ray        *ray;
Float       mindist, *maxdist;
{
        Vector  currpos;
        Float   offset, currmaxdist;
        int     x, y, sum;
        int     hit=FALSE;
        Float   Abs_t_X, Delta_t_X;
        int     incX, boundX;
        Float   Abs_t_Y, Delta_t_Y;
        int     incY, boundY;
        Float   Delta_w_X_z, Delta_w_Y_z;
        Float   wx_z, wy_z, currpos_z;
        int     triangleorder;
        Float   Rx, Ry, Rsum;
        

        MountainTests++;

        /*
         * Find out where the projected ray hits the grid.
         */
        VecAddScaled(ray->pos, mindist, ray->dir, &currpos);
        if (OutOf2DBounds(currpos, mountain->D2Bounds)) {
             offset = *maxdist;
             if (!MountainBoundsIntersect(ray, mountain->D2Bounds, mindist, &offset))
                        return FALSE;
             else {
                  VecAddScaled(ray->pos, offset, ray->dir, &currpos);
                  if (equal(currpos.x, mountain->D2Bounds[LOW][X]) ||
                      equal(currpos.y, mountain->D2Bounds[LOW][Y]))
                            triangleorder = 1;
                  else triangleorder = -1;
              }
        } else {
             offset = mindist;
             Rx = x2grid(mountain, currpos.x);
             Ry = y2grid(mountain, currpos.y);
             Rsum = Rx + Ry;
             sum = (int) Rx + (int) Ry;
             if ((Rsum - sum) <= 0.5) triangleorder = 1;
             else triangleorder = -1;
        }

        x = x2grid(mountain, currpos.x);
        if (x == mountain->GridMax) x--;

        if (fabs(ray->dir.x) < EPSILON) {
                Abs_t_X   = FAR_AWAY;
                Delta_t_X = 0;
        } else if (ray->dir.x < 0.) {
                Abs_t_X     = offset + (grid2x(mountain, x) - currpos.x) / ray->dir.x;
                Delta_t_X   = mountain->XSize / (-ray->dir.x);
                incX = -1;  boundX = -1;
        } else {
                Abs_t_X     = offset + (grid2x(mountain, x+1) - currpos.x) / ray->dir.x;
                Delta_t_X   = mountain->XSize / ray->dir.x;
                incX = 1;  boundX = mountain->GridMax;
        }

        y = y2grid(mountain, currpos.y);
        if (y == mountain->GridMax) y--;

        if (fabs(ray->dir.y) < EPSILON) {
                Abs_t_Y   = FAR_AWAY;
                Delta_t_Y = 0;
         } else if (ray->dir.y < 0.) {
                Abs_t_Y     = offset + (grid2y(mountain, y) - currpos.y) / ray->dir.y;
                Delta_t_Y   = mountain->YSize / (-ray->dir.y);
                incY = -1;  boundY = -1;
        } else {
                Abs_t_Y     = offset + (grid2y(mountain, y+1) - currpos.y) / ray->dir.y;
                Delta_t_Y   = mountain->YSize / ray->dir.y;
                incY = 1;  boundY = mountain->GridMax;
        }


        Delta_w_X_z = Delta_t_X * ray->dir.z;
        Delta_w_Y_z = Delta_t_Y * ray->dir.z;

        wx_z = ray->pos.z + Abs_t_X * ray->dir.z;
        wy_z = ray->pos.z + Abs_t_Y * ray->dir.z;

        currpos_z   = currpos.z;

        /*
         * Traverse grid with DDA-algorithm
         */
        while(TRUE) {
            if (Abs_t_X < Abs_t_Y) {
                currmaxdist = min(Abs_t_X, *maxdist);
                if (CheckGrid(mountain, x, y, ray, offset, &currmaxdist, triangleorder, currpos_z, wx_z)) {
                    *maxdist = currmaxdist;
                    return TRUE;
                }
                x += incX;
                if (x == boundX) break;
                triangleorder = incX;
                offset = currmaxdist;
                currpos_z  = wx_z;
                Abs_t_X   += Delta_t_X;
                wx_z      += Delta_w_X_z;
            } else {
                currmaxdist = min(Abs_t_Y, *maxdist);
                if (CheckGrid(mountain, x, y, ray, offset, &currmaxdist, triangleorder, currpos_z, wy_z)) {
                    *maxdist = currmaxdist;
                    return TRUE;
                }
                y += incY;
                if (y == boundY) break;
                triangleorder = incY;
                offset = currmaxdist;
                currpos_z   = wy_z;
                Abs_t_Y   += Delta_t_Y;
                wy_z      += Delta_w_Y_z;
           }
        }
        return FALSE;
}


int MountainBoundsIntersect(ray, D2Bounds, mindist, maxdist)
Ray     *ray;
Float    D2Bounds[2][2];
Float    mindist;
Float   *maxdist;
{

    Float    t, tmin, tmax;
    Float    dir, pos;


    tmax = *maxdist;
    tmin =  mindist;

    dir = ray->dir.x;
    pos = ray->pos.x;

    if (dir < 0) {
        t = (D2Bounds[LOW][X] - pos) / dir;
        if (t < tmin) return FALSE;
        if (t <= tmax) tmax = t;
        t = (D2Bounds[HIGH][X] - pos) / dir;
        if (t >= tmin) {
            if (t > tmax) return FALSE;
            tmin = t;
        }
    } else if (dir > 0) {
        t = (D2Bounds[HIGH][X] - pos) / dir;
        if (t < tmin) return FALSE;
        if (t <= tmax) tmax = t;
        t = (D2Bounds[LOW][X] - pos) / dir;
        if (t >= tmin) {
            if (t > tmax) return FALSE;
            tmin = t;
        }
    } else if (pos < D2Bounds[LOW][X] ||
               pos > D2Bounds[HIGH][X]) return FALSE;

    dir = ray->dir.y;
    pos = ray->pos.y;

    if (dir < 0) {
        t = (D2Bounds[LOW][Y] - pos) / dir;
        if (t < tmin) return FALSE;
        if (t <= tmax) tmax = t;
        t = (D2Bounds[HIGH][Y] - pos) / dir;
        if (t >= tmin) {
            if (t > tmax) return FALSE;
            tmin = t;
        }
    } else if (dir > 0) {
        t = (D2Bounds[HIGH][Y] - pos) / dir;
        if (t < tmin) return FALSE;
        if (t <= tmax) tmax = t;
        t = (D2Bounds[LOW][Y] - pos) / dir;
        if (t >= tmin) {
            if (t > tmax) return FALSE;
            tmin = t;
        }
    } else if (pos < D2Bounds[LOW][Y] ||
               pos > D2Bounds[HIGH][Y]) return FALSE;


    if (tmin == mindist) {
        if (tmax < *maxdist) {
            *maxdist = tmax;
            return TRUE;
        }
    } else {
        if (tmin < *maxdist) {
            *maxdist = tmin;
            return TRUE;
        }
    }
    return FALSE;
}

#define CHECKCACHE(cache1, cache2, Rx, Ry, calclist, xx, yy)    \
    for (ii = 0; ii <= 1; ii++, Rx++, Ry-=2) {                  \
        for (jj = 0; jj <= 1; jj++, Ry++) {                     \
            if ((cache1->x == Rx) && (cache1->y == Ry))         \
                calclist->Data[xx+ii][yy+jj] = cache1->z;       \
            else if ((cache2->x == Rx) && (cache2->y == Ry))    \
                calclist->Data[xx+ii][yy+jj] = cache2->z;       \
            else calclist->ToDo[xx+ii][yy+jj] = TRUE;           \
            cache1++;                                           \
            cache2--;                                           \
        }                                                       \
        cache1--;                                               \
        cache2+=3;                                              \
    }


/*
 * Calculate heightvalues at four corners of grid
 * and check if ray possibly intersects the two triangles
 * defined by the four points
 */
int CheckGrid(m, x, y, ray, offset, maxdist, triangleorder, z_near, z_far)
Mountain    *m;
int          x, y;
Ray         *ray;
Float        offset;
Float       *maxdist;
int          triangleorder;
Float        z_near, z_far;
{

   HeightCache *Cache1, *Cache2;
   int          NX, NY, TX, TY, TTX, TTY, XP, YP;
   int          level, i, j, xx, yy, ii, jj;
   int          XShort, XFar, YShort, YFar;
   CalcList    *TCL, *TCLP;
   Vector       point1, point2, point3, edge1, edge2, edge3;
   Float        scale;

   static CalcList  EvolveLevelData;



   level = m->IterationLevel - m->EvolveLevel - 1;
   if (level >= 0) {
       TCL = m->ToCalcList;
       TCL->x = x;
       TCL->y = y;
       for (i = 0; i <= 1; i++)
            for (j = 0; j<= 1; j++)
                 TCL->ToDo[i][j] = FALSE;
       /* 
        * Check Cache if cell is available 
        */
       Cache1 = &m->Cache1[level][x+y];
       Cache2 = &m->Cache2[level][m->MaxSize[level + 1] + x - y];
       xx = x;  yy = y;
       CHECKCACHE(Cache1, Cache2, xx, yy, TCL, 0, 0);
       xx= x; yy=y;
       level--;
       /*  
        * Calculate the coordinates of the cells necsarry for calculating
        * the top cell, check cache if they are available 
        */
       for (TCLP = TCL, TCL = TCL->Next; TCL; TCLP = TCL, TCL = TCL->Next,
            level--, xx = NX, yy = NY) {
              NX = TCL->x = xx>>1;
              NY = TCL->y = yy>>1;
              for (i = 0; i <= 2; i++)
                   for (j = 0; j <= 2; j++)
                        TCL->ToDo[i][j] = FALSE;
              for (i = 0; i <= 2; i++)
                  for (j = 0; j <= 2; j++)
                      if (TCLP->ToDo[i][j]) {
                          TX = (xx+i)>>1;  TY = (yy+j)>>1;
                          Cache1 = &m->Cache1[level][TX+TY];
                          Cache2 = &m->Cache2[level][m->MaxSize[level + 1]+TX-TY];
                          TTX = TX - NX;  TTY = TY - NY;
                          CHECKCACHE(Cache1, Cache2, TX, TY, TCL, TTX, TTY);
                      }
       }
       xx = xx>>1;  yy = yy>>1;
       for (i = 0; i <= min(2, m->MaxSize[0] - xx); i++)
           for (j = 0; j <= min(2, m->MaxSize[0] - yy); j++)
               EvolveLevelData.Data[i][j] = m->EvolvedData[xx + i][yy + j];
       level = 0;
       XP    = xx; YP    = yy;
       /*
        * Calculate heightvalues of corners of cell who were
        * not available in cache 
        */
       for (TCL = TCLP, TCLP = &EvolveLevelData; TCL; TCLP = TCL, TCL = TCL->Previous,
            level++, XP = xx, YP = yy) {
              xx = TCL->x;  yy = TCL->y;
              Cache1 = &m->Cache1[level][xx + yy];
              Cache2 = &m->Cache2[level][m->MaxSize[level + 1] + xx - yy];
              scale = m->Scale * pow((m->StartGridSize/ (Float) m->MaxSize[level + 1]), m->FractalDimension);
              for (i = 0; i<=2; i++) {
                  for (j = 0; j <= 2; j++) {
                      if (TCL->ToDo[i][j]) {
                          XShort = (xx + 1 + i)>>1;
                          XFar   = ((xx + i)|1) - XShort;
                          YShort = (yy + 1 + j)>>1;
                          YFar   = ((yy + j)|1) - YShort;
                          XShort -= XP;
                          XFar -= XP;
                          YShort -= YP;
                          YFar -= YP;
                          TCL->Data[i][j] = (3*(TCLP->Data[XShort][YFar]
                                                + TCLP->Data[XFar][YShort]
                                                + 3 * TCLP->Data[XShort][YShort])
                                             + TCLP->Data[XFar][YFar]) * 0.0625 +
                                             GAUSS(m->Number, xx + i, yy + j, level + m->EvolveLevel + 1) * scale;

                          if (WhichCache[i][j]) {
                              Cache1->x = xx + i;
                              Cache1->y = yy + j;
                              Cache1->z = TCL->Data[i][j];
                          } else {
                              Cache2->x = xx + i;
                              Cache2->y = yy + j;
                              Cache2->z = TCL->Data[i][j];
                          }
                      }
                      Cache1++;
                      Cache2--;
                  }
                  Cache1 -= 2;
                  Cache2 += 4;
              }
          }
   } else {
       for (i = 0; i <= 1; i++)
            for (j = 0; j <= 1; j++)
                 m->ToCalcList->Data[i][j] = m->EvolvedData[x+i][y+j];
   }
   TCL = m->ToCalcList;
   /*
    * Check if ray crosses cell at right height 
    */
   if ((min(z_near, z_far) > max(max(max(TCL->Data[0][0], TCL->Data[0][1]), TCL->Data[1][0]), TCL->Data[1][1])) ||
       (max(z_near, z_far) < min(min(min(TCL->Data[0][0], TCL->Data[0][1]), TCL->Data[1][0]), TCL->Data[1][1]))) return FALSE;
   /*
    * Test for ray-triangle intersection 
    */
   for (i = 1; i <=2 ; i++, triangleorder += 2)
       if (triangleorder == 1) {
           point1.x = (x * m->XSize) + m->D2Bounds[LOW][X];
           point1.y = (y * m->YSize) + m->D2Bounds[LOW][Y];
           point1.z = TCL->Data[0][0];
           point2.x = ((x+1) * m->XSize) + m->D2Bounds[LOW][X];
           point2.y = (y * m->YSize) + m->D2Bounds[LOW][Y];
           point2.z = TCL->Data[1][0];
           point3.x = (x * m->XSize) + m->D2Bounds[LOW][X];
           point3.y = ((y+1) * m->YSize) + m->D2Bounds[LOW][Y];
           point3.z = TCL->Data[0][1];
           if (IntersectTriangle(m, ray, &point1, &point2, &point3, offset, maxdist)) return TRUE;
       } else {
           point1.x = ((x+1) * m->XSize) + m->D2Bounds[LOW][X];
           point1.y = ((y+1) * m->YSize) + m->D2Bounds[LOW][Y];
           point1.z = TCL->Data[1][1];
           point2.x = (x * m->XSize) + m->D2Bounds[LOW][X];
           point2.y = ((y+1) * m->YSize) + m->D2Bounds[LOW][Y];
           point2.z = TCL->Data[0][1];
           point3.x = ((x+1) * m->XSize) + m->D2Bounds[LOW][X];
           point3.y = (y * m->YSize) + m->D2Bounds[LOW][Y];
           point3.z = TCL->Data[1][0];
           if (IntersectTriangle(m, ray, &point1, &point2, &point3, offset, maxdist)) return TRUE;
       }
   return FALSE;
}

/* 
 * Check if ray intersects triangle 
 */

int IntersectTriangle(m, ray, point1, point2, point3, mindist, maxdist)
Mountain   *m;
Ray        *ray;
Vector     *point1, *point2, *point3;
Float       mindist, *maxdist;
{
    Vector  edge1, edge2, edge3;
    Float   k, s;
    Vector  normal, absnorm;
    Vector  pos, dir;
    Float   qi1, qi2, b0, b1, b2;


    VecSub(*point2, *point1, &edge1);
    VecSub(*point3, *point2, &edge2);
    VecSub(*point1, *point3, &edge3);

    /* 
     * Find plane normal. 
     */
    VecCross(&edge1, &edge2, &normal);

    pos = ray->pos;
    dir = ray->dir;
    /*
     * Plane intersection.
     */
    k = dotp(&normal, &dir);
    if (fabs(k) < EPSILON)
        return FALSE;
    s = (dotp(&normal, point1) - dotp(&normal, &pos)) / k;
    if (s < mindist || s > *maxdist)
                return FALSE;

    /*
     * Find "dominant" part of normal vector.
     */
    absnorm.x = fabs(normal.x);
    absnorm.y = fabs(normal.y);
    absnorm.z = fabs(normal.z);

    if (absnorm.x > absnorm.y && absnorm.x > absnorm.z) {
        qi1 = pos.y + s * dir.y;
        qi2 = pos.z + s * dir.z;
        b0 = (edge2.y * (qi2 - point2->z) -
              edge2.z * (qi1 - point2->y)) / normal.x;
        if (b0 < 0. || b0 > 1.) return FALSE;
        b1 = (edge3.y * (qi2 - point3->z) -
              edge3.z * (qi1 - point3->y)) / normal.x;
        if (b1 < 0. || b1 > 1.) return FALSE;
        b2 = (edge1.y * (qi2 - point1->z) -
              edge1.z * (qi1 - point1->y)) / normal.x;
        if (b2 < 0. || b2 > 1.) return FALSE;
    } else if (absnorm.y > absnorm.z) {
        qi1 = pos.x + s * dir.x;
        qi2 = pos.z + s * dir.z;
        b0 = (edge2.z * (qi1 - point2->x) -
              edge2.x * (qi2 - point2->z)) / normal.y;
        if (b0 < 0. || b0 > 1.) return FALSE;
        b1 = (edge3.z * (qi1 - point3->x) -
              edge3.x * (qi2 - point3->z)) / normal.y;
        if (b1 < 0. || b1 > 1.) return FALSE;
        b2 = (edge1.z * (qi1 - point1->x) -
              edge1.x * (qi2 - point1->z)) / normal.y;
        if (b2 < 0. || b2 > 1.) return FALSE;
    } else {
        qi1 = pos.x + s * dir.x;
        qi2 = pos.y + s * dir.y;
        b0 = (edge2.x * (qi2 - point2->y) -
              edge2.y * (qi1 - point2->x)) / normal.z;
        if (b0 < 0. || b0 > 1.) return FALSE;
        b1 = (edge3.x * (qi2 - point3->y) -
              edge3.y * (qi1 - point3->x)) / normal.z;
        if (b1 < 0. || b1 > 1.) return FALSE;
        b2 = (edge1.x * (qi2 - point1->y) -
              edge1.y * (qi1 - point1->x)) / normal.z;
        if (b2 < 0. || b2 > 1.) return FALSE;
    }
    MountainHits++;
    *maxdist = s;
    m->Normal = normal;
    VecAddScaled(ray->pos, s, ray->dir, &m->IntersectPos);

    return TRUE;
}

#define VecEqual(a,b) ((equal((a).x, (b).x)) && (equal((a).y, (b).y)) && (equal((a).z, (b).z)))

/*
 * Calculate normal on triangle in pos
 */
int MountainNormal(m, pos, nrm, gnrm)
Mountain *m;
Vector   *pos, *nrm, *gnrm;
{
    if (!VecEqual(*pos, m->IntersectPos))
        RLerror(RL_ABORT, "Can't compute mountain normal in point other than the last intersection.\n");
        *gnrm = m->Normal;

        VecNormalize(gnrm);
        *nrm = *gnrm;
        return FALSE;
}


char *MountainName()
{
        return mountainName;
}

void MountainStats(tests, hits)
unsigned long *tests, *hits;
{
        *tests = MountainTests;
        *hits = MountainHits;
}

void MountainUV(m, pos, norm, uv, dpdu, dpdv)
Mountain   *m;
Vector     *pos, *norm, *dpdu, *dpdv;
Vec2d      *uv;
{
        uv->u = pos->x;
        uv->v = pos->y;
        if (dpdu) {
              dpdu->x = 1.;
              dpdu->y = dpdu->z = 0.;
              dpdv->x = dpdv->z = 0.;
              dpdv->y = 1.;
        }
}

void MountainBounds(m, bounds)
Mountain *m;
Float     bounds[2][3];
{
        bounds[LOW][X] = m->D2Bounds[LOW][X];
        bounds[LOW][Y] = m->D2Bounds[LOW][Y];
        bounds[HIGH][X] = m->D2Bounds[HIGH][X];
        bounds[HIGH][Y] = m->D2Bounds[HIGH][Y];
        bounds[LOW][Z] = -FAR_AWAY;
        bounds[HIGH][Z] = FAR_AWAY;
}

void MountainMethodRegister(meth)
UserMethodType  meth;
{
        if (iMountainMethods)
            iMountainMethods->user = meth;
}


Methods *MountainMethods()
{
        if (iMountainMethods == (Methods *)NULL) {
                iMountainMethods = MethodsCreate();
                iMountainMethods->create = (GeomCreateFunc *)MountainCreate;
                iMountainMethods->methods = MountainMethods;
                iMountainMethods->name = MountainName;
                iMountainMethods->intersect = MountainIntersect;
                iMountainMethods->normal = MountainNormal;
                iMountainMethods->uv = MountainUV;
                iMountainMethods->bounds = MountainBounds;
                iMountainMethods->stats = MountainStats;
                iMountainMethods->checkbounds = FALSE;
                iMountainMethods->closed = FALSE;
        }
        return iMountainMethods;
}