GEOMETRY.C

Go to the documentation of this file.

// file GEOMETRY.CPP


#define MODULE_GEO

#include "render.h"


static unsigned char *AllocateShadowBuffer(void);
static void ApplyVertexEffect(long effect, vertex *v, long N, double time,
                              double vmin[], double vmax[], ivector Og);
static void FixUpMapEdge(GLfloat a[]);
static void NewPoint(double plane, vector p1, vector p2, vector p, 
                     vector2, vector2, vector2);



static void NewPoint(double dp, vector v1, vector v2, vector p, 
                               vector2 uv1, vector2 uv2, vector2 uv){
 double mu;
 mu=(v2[1]-v1[1]);
 if(fabs(mu) < 1.0)mu=0.0;
 else mu=(dp-v1[1])/mu;
 p[1]=dp;
 p[0]=v1[0]+mu*(v2[0]-v1[0]);
 p[2]=v1[2]+mu*(v2[2]-v1[2]);
 uv[0] = uv1[0]+mu*(uv2[0]-uv1[0]);
 uv[1] = uv1[1]+mu*(uv2[1]-uv1[1]);
}

static void  NewMapPoint(double dp, vector a1, vector a2, 
                         double u1, double u2, 
                         double v1, double v2,
                         double w1, double w2, 
                         vector r){
 double mu;
 mu=(a2[1]-a1[1]);
 if(fabs(mu) < 1.0)mu=0.0;
 else mu=(dp-a1[1])/mu;
 r[0] = u1 + mu*(u2 - u1);
 r[1] = v1 + mu*(v2 - v1);
 r[2] = w1 + mu*(w2 - w1);
}

short ClipObject(long O){
  vertex   *v;
  face     *f;
  double y1,y2,y3,dp;
  long i,i1,i2,i3,xv,xf,sz,vp,fp;
  vector2 uv1,uv2,uva,uvb;
  vector  a,b;
/* any faces that span the projection plane have to be turned into
   three faces with new vertices on the projection plane.
*/
//  dp=1.000;
  dp=FRONT_CLIPPING_PLANE;
  xv=0; xf=0;
  v=Object[O].Vbase;
  f=Object[O].Fbase;
  for(i=0;i<Object[O].NoFaces;i++){
     i1= -1;
     y1=v[f[i].V[0]].p[1];
     y2=v[f[i].V[1]].p[1];
     y3=v[f[i].V[2]].p[1];
     if     (y1 <= dp && y2 <=dp && y3 >  dp)i1 = 2;
     else if(y1 <= dp && y2 > dp && y3 <= dp)i1 = 1;
     else if(y1 <= dp && y2 > dp && y3 >  dp)i1 = 0;
     else if(y1 >  dp && y2 <=dp && y3 <= dp)i1 = 0;
     else if(y1 >  dp && y2 <=dp && y3 >  dp)i1 = 1;
     else if(y1 >  dp && y2 > dp && y3 <= dp)i1 = 2;
     if(i1 >= 0){ /* calculate number of extra faces and vertices */
       xv += 2; xf += 2;
     }
  }
  if(xv == 0)return OK;  /* no clipping necessary */
  //if(debug != NULL)fprintf(debug,"Clipping necessary extra vertices %ld faces %ld\n",xv,xf);
#if 0
fprintf(debug,"extra v=%ld face=%ld\n",xv,xf);
for(i=0;i<Object[O].NoFaces;i++){
int j;
fprintf(debug,"Face %ld vertices (%ld %ld %ld)\n",i,(f+i)->V[0],(f+i)->V[1],(f+i)->V[2]);
for(j=0;j<3;j++){
fprintf(debug,"(Vvert %8.1lf %8.1lf %8.1lf) (mapping  %lf %lf)\n",
(v+(f+i)->V[j])->p[0],(v+(f+i)->V[j])->p[1],(v+(f+i)->V[j])->p[2],
(f+i)->uu[j],(f+i)->vv[j]);
}
}
#endif
  /* resize the arrays and do it */
  sz=(Object[O].NoVertices+xv)*(long)sizeof(vertex);
  if((v=(vertex  *)X__Realloc(Object[O].Vbase,sz)) == NULL)goto X__ERROR;
  Object[O].Vbase=v;
  sz=(Object[O].NoFaces+xf)*(long)sizeof(face);
  if((f=(face  *)X__Realloc(Object[O].Fbase,sz)) == NULL)goto X__ERROR;
  Object[O].Fbase=f;
  vp=Object[O].NoVertices;
  fp=Object[O].NoFaces;
  for(i=0;i<Object[O].NoFaces;i++){
     i1= -1; i2= -1; i3= -1;
     y1=v[f[i].V[0]].p[1];
     y2=v[f[i].V[1]].p[1];
     y3=v[f[i].V[2]].p[1];
     if     (y1 <= dp && y2 <=dp && y3 >  dp){i1 = 2; i2=0; i3=1;}
     else if(y1 <= dp && y2 > dp && y3 <= dp){i1 = 1; i2=0; i3=2;}
     else if(y1 <= dp && y2 > dp && y3 >  dp){i1 = 0; i2=1; i3=2;}
     else if(y1 >  dp && y2 <=dp && y3 <= dp){i1 = 0; i2=1; i3=2;}
     else if(y1 >  dp && y2 <=dp && y3 >  dp){i1 = 1; i2=0; i3=2;}
     else if(y1 >  dp && y2 > dp && y3 <= dp){i1 = 2; i2=0; i3=1;}
     if(i1 >= 0){ /* calculate the new vertex positions on viewing plane */
       v[vp].x=v[vp].y=v[vp+1].x=v[vp+1].y=0.0; 
        /* calculate intersection point between points i1 - i2 and i1 - i3
           and place new vertices to interpolated position on viewing plane*/
       uv1[0]=v[f[i].V[i1]].u; uv1[1]=v[f[i].V[i1]].v;
       uv2[0]=v[f[i].V[i2]].u; uv2[1]=v[f[i].V[i2]].v;
       NewPoint(dp,v[f[i].V[i1]].p,v[f[i].V[i2]].p,v[vp].p,uv1,uv2,uva);
       v[vp].u=uva[0]; v[vp].v=uva[1];

       uv1[0]=v[f[i].V[i1]].u; uv1[1]=v[f[i].V[i1]].v;
       uv2[0]=v[f[i].V[i3]].u; uv2[1]=v[f[i].V[i3]].v;
       NewPoint(dp,v[f[i].V[i1]].p,v[f[i].V[i3]].p,v[vp+1].p,uv1,uv2,uvb);
       v[vp+1].u=uvb[0]; v[vp+1].v=uvb[1];

       NewMapPoint(dp,v[f[i].V[i1]].p,v[f[i].V[i2]].p,
                      f[i].uu[i1],f[i].uu[i2],
                      f[i].vv[i1],f[i].vv[i2],
                      f[i].ww[i1],f[i].ww[i2],a);
       NewMapPoint(dp,v[f[i].V[i1]].p,v[f[i].V[i3]].p,
                      f[i].uu[i1],f[i].uu[i3],
                      f[i].vv[i1],f[i].vv[i3],
                      f[i].ww[i1],f[i].ww[i3],b);
//fprintf(debug,"Splint face %ld  (%ld %ld %ld) between [%ld %ld] and [%ld %ld] \n",i,i1,i2,i3, f[i].V[i1],f[i].V[i2],
//f[i].V[i1],f[i].V[i3]); 
//fprintf(debug,"New mapping a %lf %lf\n",a[0],a[1]);
//fprintf(debug,"New mapping b %lf %lf\n",b[0],b[1]);

       f[fp].V[0]=vp;
       f[fp].V[1]=f[i].V[i2];
       f[fp].V[2]=f[i].V[i3];
       f[fp].uu[0]=a[0]; f[fp].vv[0]=a[1]; f[fp].ww[0]=a[2]; 
       f[fp].uu[1]=f[i].uu[i2]; f[fp].vv[1]=f[i].vv[i2]; f[fp].ww[1]=f[i].ww[i2]; 
       f[fp].uu[2]=f[i].uu[i3]; f[fp].vv[2]=f[i].vv[i3]; f[fp].ww[2]=f[i].ww[i3]; 
       
       f[fp+1].V[0]=vp;
       f[fp+1].V[1]=f[i].V[i3];
       f[fp+1].V[2]=vp+1;
       f[fp+1].uu[0]=a[0]; f[fp+1].vv[0]=a[1]; f[fp+1].ww[0]=a[2]; 
       f[fp+1].uu[1]=f[i].uu[i3]; f[fp+1].vv[1]=f[i].vv[i3]; f[fp+1].ww[1]=f[i].ww[i3]; 
       f[fp+1].uu[2]=b[0]; f[fp+1].vv[2]=b[1]; f[fp+1].ww[2]=b[2]; 

       f[i].V[i2]=vp;
       f[i].V[i3]=vp+1;
       f[i].uu[i2]=a[0]; f[i].vv[i2]=a[1]; f[i].ww[i2]=a[2]; 
       f[i].uu[i3]=b[0]; f[i].vv[i3]=b[1]; f[i].ww[i3]=b[2]; 
       f[fp+1].color[0]=f[fp].color[0]=f[i].color[0];
       f[fp+1].color[1]=f[fp].color[1]=f[i].color[1];
       f[fp+1].color[2]=f[fp].color[2]=f[i].color[2];
       f[fp+1].ffmap=f[fp].ffmap=f[i].ffmap;
       f[fp+1].ffmat=f[fp].ffmat=f[i].ffmat;
       f[fp+1].bFmc=f[fp].bFmc=f[i].bFmc;
       f[fp+1].bShiny=f[fp].bShiny=f[i].bShiny;
       f[fp+1].bSmooth=f[fp].bSmooth=f[i].bSmooth;
       f[fp+1].pn[0][0] = f[fp].pn[0][0] = f[i].pn[0][0];
       f[fp+1].pn[0][1] = f[fp].pn[0][1] = f[i].pn[0][1];
       f[fp+1].pn[0][2] = f[fp].pn[0][2] = f[i].pn[0][2];
       f[fp+1].pn[1][0] = f[fp].pn[1][0] = f[i].pn[1][0];
       f[fp+1].pn[1][1] = f[fp].pn[1][1] = f[i].pn[1][1];
       f[fp+1].pn[1][2] = f[fp].pn[1][2] = f[i].pn[1][2];
       f[fp+1].pn[2][0] = f[fp].pn[2][0] = f[i].pn[2][0];
       f[fp+1].pn[2][1] = f[fp].pn[2][1] = f[i].pn[2][1];
       f[fp+1].pn[2][2] = f[fp].pn[2][2] = f[i].pn[2][2];
       f[fp].c=f[fp+1].c=0.0;
       f[fp].n[0]=f[fp].n[1]=f[fp].n[2]=f[fp+1].n[0]
                 =f[fp+1].n[1]=f[fp+1].n[2]=0.0;
       vp += 2; fp += 2;
     }
  }
  Object[O].NoVertices = vp;
  Object[O].NoFaces = fp;
#if 0
fprintf(debug,"After v=%ld face=%ld\n",Object[O].NoVertices,Object[O].NoFaces);
for(i=0;i<Object[O].NoFaces;i++){
int j;
fprintf(debug,"Face %ld vertices (%ld %ld %ld)\n",i,(f+i)->V[0],(f+i)->V[1],(f+i)->V[2]);
for(j=0;j<3;j++){
fprintf(debug,"(Vvert %8.1lf %8.1lf %8.1lf) (mapping  %lf %lf)\n",
(v+(f+i)->V[j])->p[0],(v+(f+i)->V[j])->p[1],(v+(f+i)->V[j])->p[2],
(f+i)->uu[j],(f+i)->vv[j]);
}
} 
#endif
  return OK;
  X__ERROR:
  Render_Message(4,0,NULL);
  return FAIL;
}

long IShadow(long scanline, long v1x, long v1y, long v2x, long v2y,
             long *xi){
 double doff,ddxy;
 long dx,dy,off,dxy;
 if((v1y < scanline)  && (v2y < scanline))return 0;
 if((v1y > scanline)  && (v2y > scanline))return 0;
 dy = (v2y) - (v1y);
 if(dy == 0)return 0;
 dx = (v2x) - (v1x);
 off=scanline - v1y;
 doff=off;
 ddxy=dx*doff/dy;
 dxy=ddxy;
 *xi=(v1x)+dxy;
 return 1;
}

void ShadowFillLine(long scanline, long vx[3], long vy[3],
                      long *leftedge, long *rightedge){
 long x1,x2,x3,i1,i2,i3;
 i1=IShadow(scanline,vx[0],vy[0],vx[1],vy[1],&x1);
 i2=IShadow(scanline,vx[1],vy[1],vx[2],vy[2],&x2);
 i3=IShadow(scanline,vx[2],vy[2],vx[0],vy[0],&x3);
 if((i1 == 1) && (i2 == 1) && (i3 == 1))
 {*leftedge = min(min(x1,x2),x3);
  *rightedge= max(max(x1,x2),x3);}
 else if((i1 == 1) && (i2 == 1))
 {*leftedge = min(x1,x2);
  *rightedge= max(x1,x2);}
 else if((i2 == 1) && (i3 == 1))
 {*leftedge = min(x2,x3);
  *rightedge= max(x2,x3);}
 else if((i3 == 1) && (i1 == 1))
 {*leftedge = min(x3,x1);
  *rightedge= max(x3,x1);}
 else {
   Render_Message(5,0,NULL);
   return;
 }
 *leftedge  = max(0,min(XMAX-1,*leftedge));
 *rightedge = min(XMAX-1,max(0,*rightedge));
 return;
}

static unsigned char *AllocateShadowBuffer(void){
 long i,screensize;
 unsigned char *s,*screen;
 screensize = (long)XMAX * (long)YMAX;
 if((s = (unsigned char *)X__Malloc(screensize)) == NULL){
   Render_Message(14,0,NULL);
   return NULL;
 }
 screen=s; for(i=0;i<screensize;i++)*screen++ = 0;
 return s;
}

void ShadowGround(void){
 long i,j,k,l,m;
 long vi[3];
 long xx[3],yy[3],xxmin,xxmax,yymin,yymax;
 face    *f;
 vertex  *v;
 matl    *mat;
 unsigned char mattr;
 vector p,t1;
 double mu,nu,de,temp4;
 UCHAR  *screen;
 if((VGA_screen=AllocateShadowBuffer()) == NULL)return;
 screen=VGA_screen;
 for(k=0;k<Nlights;k++){
   if(Lights[k].type == NORMAL || Lights[k].type == SHADOWS){
     for(i=0;i<ObjectCount;i++){
       if(!Object[i].in_use)continue;
       if(!Object[i].cast_shadow)continue;
       v=Object[i].Vbase;
       f=Object[i].Fbase;
       mat=Object[i].Tbase;
       for(j=0;j<Object[i].NoFaces;j++){
         if((f+j)->ffmat >= 0 && (f+j)->ffmat < Object[i].NoMats){
           if(mat != NULL)mattr=(mat + f->ffmat)->transp;
           else mattr=0;
         } else mattr=0;
         if(mattr > 0)continue; /* glass no shadow */
         vi[0]=(f+j)->V[0];
         vi[1]=(f+j)->V[1];
         vi[2]=(f+j)->V[2];
         for(l=0;l<3;l++){ /* project each vertex to shadow plane */
           p[0]=(v+vi[l])->p[0];
           p[1]=(v+vi[l])->p[1];
           p[2]=(v+vi[l])->p[2];
           vecsub(Ground.p,Lights[k].p,t1);
           nu=dot(Ground.n,t1);
           vecsub(p,Lights[k].p,t1);
           de=dot(Ground.n,t1);
           if(fabs(de) < 1.e-5)goto FACEBREAK;
           mu=nu/de;
           if(mu < 1.0)goto FACEBREAK;
           vecsub(p,Lights[k].p,t1);
           vecscale(mu,t1,t1);
           vecsum(Lights[k].p,t1,t1);
           if(t1[1] > 1.000){
             temp4=(long)(scalex*t1[0]/(t1[1]));
             xx[l]=Xcentre+temp4;
             temp4=(long)(scaley*t1[2]/(t1[1]));
             yy[l]=Ycentre-temp4;
           }
           else goto FACEBREAK;
         }
         xxmin=min(xx[0],min(xx[1],xx[2]));
         xxmax=max(xx[0],max(xx[1],xx[2]));
         yymin=min(yy[0],min(yy[1],yy[2]));
         yymax=max(yy[0],max(yy[1],yy[2]));
         if((yymax > yymin) && (xxmax > xxmin))
         for(l=max(0,yymin);l<=min(yymax,YMAX-1);l++){
           ShadowFillLine(l,xx,yy,&xxmin,&xxmax);
           for(m=xxmin;m<=xxmax;m++) *(screen+l*(long)XMAX+m) = 1;
         }
         FACEBREAK: continue;
       }
     }
   }
 }
 return;
}

#define TOL 1.e-4;

double hitpoint(vector n, vector x, vector y, vector pb, vector pf,
                double * a, double * b){
 /* n        is normal to plane of map
    x & y    are vectors along the edges of the map  x= left y = down
    pb       is position vector to base point on map
    fp       is position vector to vertex to be mapped
    a & b    are output of mapping coordinates [0 - 1]
 */
 int k;
 vector v1,p;
 double mu,dm,detmax,det1,det2,det3,ve1,ve2,ve3,vp1,vp2,vp3,p1,p2,p3;
 VECSUB(pb,pf,v1)
 mu=DOT(n,v1);
 VECSCALE(mu,n,v1)
 VECSUM(pf,v1,p)
 ve1=x[0]; ve2=x[1]; ve3=x[2];
 vp1=y[0]; vp2=y[1]; vp3=y[2];
 p1=p[0]-pb[0];
 p2=p[1]-pb[1];
 p3=p[2]-pb[2];
 det1=ve1*vp2-vp1*ve2;
 det2=ve1*vp3-vp1*ve3;
 det3=ve2*vp3-vp2*ve3;
 k=0; detmax=TOL;
 if((dm=fabs(det1)) > detmax){k=1; detmax=dm;}
 if((dm=fabs(det2)) > detmax){k=2; detmax=dm;}
 if((dm=fabs(det3)) > detmax){k=3; detmax=dm;}
 if(k == 1){
   *a=( vp2*p1-vp1*p2)/det1;
   *b=(-ve2*p1+ve1*p2)/det1;
 }
 else if(k == 2){
   *a=( vp3*p1-vp1*p3)/det2;
   *b=(-ve3*p1+ve1*p3)/det2;
 }
 else if(k == 3){
   *a=( vp3*p2-vp2*p3)/det3;
   *b=(-ve3*p2+ve2*p3)/det3;
 }
 else {
   *a = *b = 0.0;
 }
 return mu;
}

BOOL hitpoint2(vector n, vector x, vector y, vector pb, vector pf,
                double * a, double * b){
 /* n        is normal to plane of map
    x & y    are vectors along the edges of the map  x= left y = down
    pb       is position vector to base point on map
    fp       is position vector to vertex to be mapped
    a & b    are output of mapping coordinates [0 - 1]
 */
 // an alternative hitpoint function that assumes that the point of
 // intersection pf lies in the plane !!
 vector c;
 double dd,xx,yy,xy,cx,cy;
 VECSUB(pb,pf,c)
 xx=DOT(x,x);
 yy=DOT(y,y);
 xy=DOT(x,y);
 cy=DOT(c,y);
 cx=DOT(c,x);
 dd=xx*yy-xy*xy;
 if(fabs(dd) < 1.e-10){
   *a = *b = 0.0;
   return FALSE;
 }
 dd = 1.0/dd;
 *a=(cx*yy-cy*xy)*dd;
 *b=(cy*xx-cx*xy)*dd;
 return TRUE;
}

#define RTD 0.159154943

#ifdef _SUNSTYLE
double hitcylinder(n,x,y,p0,p,a,b) vector n; vector x; vector y;
vector p0; vector p; double *a; double *b; {
#else
double hitcylinder(vector n, vector x, vector y, vector p0, vector p,
                double * a, double * b){
#endif
  double psi,theta,ddt;
  vector p1,q,p2;
  vecsub(p,p0,p1);
  psi =dot(p1,y)/dot(y,y);
  *b = psi;
  p2[0] = psi*y[0];     p2[1] = psi*y[1];      p2[2] = psi*y[2];
  q[0] = p1[0] - p2[0];  q[1] = p1[1] - p2[1];  q[2] = p1[2] - p2[2];
  normalize(q);
  ddt=DOT(q,n);
  if     (ddt >  1.0)ddt =  1.0;
  else if(ddt < -1.0)ddt = -1.0;
  theta=acos(ddt)*RTD;
  cross(q,n,p1);
  if(dot(p1,y) <= 0.0){
    if(theta <= 0.25) *a = (0.25 - theta);
    else  *a = (1.25 - theta);
  }
  else    *a = (0.25 + theta);
  return psi;
}

double hitsphere(vector y, vector dx, vector dy, vector p0, vector p,
                double * a, double * b){
  /* p0    is at centre of map
     dx    is vector in direction of seam
     dy    is directed to south pole
     p     is point of vertex to be assigned mapping coordinates
     a     is  phi  [0 - 1] (0 along x axis)
     b     is theta [0 - 1] (0  at north pole)
  */
  vector x,z,r;
  VECCOPY(dx,x)
  normalize(x);
  VECSCALE(-1.0,dy,z)
  normalize(z);
  VECSUB(p,p0,r)
  normalize(r);
  *b = acos(DOT(r,z))*0.3183;              /*  theta/pi                    */
  *a = atan2(DOT(r,y),DOT(r,x))*0.1592;    /*  phi/(2 pi)    0   < a < 0.5 */
  if(*a < 0.0) *a = 1.0 + *a;              /*                0.5 < a < 1.0 */
  return 0.0;
}

#define UpdateNoGlassBuffer(f,i,d,id,O) {                       \
  Zdepth[i]=d; Zbuffer[i]=id; Zobject[i]=O;                     \
}                                                               \

#if NGLASS == 1
#error BADGLASS1
#define UpdateGlassBuffer(f,i,d,id,O,mattr) {                   \
     if(mattr > 0){                                             \
       if(d < Zglassd[0][i]){                                   \
         Zglassd[0][i]=d; Zglass[0][i]=id; ZglassO[0][i]=O;     \
       }                                                        \
     }                                                          \
     else{                                                      \
       if(d <= Zglassd[0][i])Zglass[0][i] = -1;                 \
       Zdepth[i]=d; Zbuffer[i]=id; Zobject[i]=O;                \
     }                                                          \
}
#elif NGLASS == 2
#error BADGLASS2
#define UpdateGlassBuffer(f,i,d,id,O,mattr) {                   \
     if(mattr > 0){                                             \
       if(d < Zglassd[0][i]-600.0){                             \
         if(O != ZglassO[0][i] || id != Zglass[0][i]){          \
           Zglassd[1][i]=Zglassd[0][i];                         \
           ZglassO[1][i]=ZglassO[0][i];                         \
           Zglass[1][i]=Zglass[0][i];                           \
         }                                                      \
         Zglassd[0][i]=d; Zglass[0][i]=id; ZglassO[0][i]=O;     \
       }                                                        \
       else if(d > Zglassd[0][i]+600.0 &&                       \
               d < Zglassd[1][i]-600.0){                        \
         Zglassd[1][i]=d; Zglass[1][i]=id; ZglassO[1][i]=O;     \
       }                                                        \
     }                                                          \
     else{                                                      \
       if(d <= Zglassd[0][i]){                                  \
         Zglass[0][i] = -1;  Zglassd[0][i]=FARAWAY;             \
         Zglass[1][i] = -1;  Zglassd[1][i]=FARAWAY;             \
       }                                                        \
       else if(d < Zglassd[1][i]){                              \
         Zglass[1][i] = -1;                                     \
         Zglassd[1][i]=FARAWAY;                                 \
       }                                                        \
       Zdepth[i]=d; Zbuffer[i]=id; Zobject[i]=O;                \
     }                                                          \
}
#elif NGLASS == 4
#define UpdateGlassBuffer(f,i,d,id,O,mattr) {                   \
     if(mattr > 0){                                             \
       if(d < Zglassd[0][i]-600.0){                             \
         if(O != ZglassO[0][i] || id != Zglass[0][i]){          \
           Zglassd[3][i]=Zglassd[2][i];                         \
           ZglassO[3][i]=ZglassO[2][i];                         \
           Zglass[3][i]=Zglass[2][i];                           \
           Zglassd[2][i]=Zglassd[1][i];                         \
           ZglassO[2][i]=ZglassO[1][i];                         \
           Zglass[2][i]=Zglass[1][i];                           \
           Zglassd[1][i]=Zglassd[0][i];                         \
           ZglassO[1][i]=ZglassO[0][i];                         \
           Zglass[1][i]=Zglass[0][i];                           \
         }                                                      \
         Zglassd[0][i]=d; Zglass[0][i]=id; ZglassO[0][i]=O;     \
       }                                                        \
       else if(d > Zglassd[0][i]+600.0 &&                       \
               d < Zglassd[1][i]-600.0){                        \
         if(O != ZglassO[1][i] || id != Zglass[1][i]){          \
           Zglassd[3][i]=Zglassd[2][i];                         \
           ZglassO[3][i]=ZglassO[2][i];                         \
           Zglass[3][i]=Zglass[2][i];                           \
           Zglassd[2][i]=Zglassd[1][i];                         \
           ZglassO[2][i]=ZglassO[1][i];                         \
           Zglass[2][i]=Zglass[1][i];                           \
         }                                                      \
         Zglassd[1][i]=d; Zglass[1][i]=id; ZglassO[1][i]=O;     \
       }                                                        \
       else if(d > Zglassd[1][i]+600.0 &&                       \
               d < Zglassd[2][i]-600.0){                        \
         if(O != ZglassO[2][i] || id != Zglass[2][i]){          \
           Zglassd[3][i]=Zglassd[2][i];                         \
           ZglassO[3][i]=ZglassO[2][i];                         \
           Zglass[3][i]=Zglass[2][i];                           \
         }                                                      \
         Zglassd[2][i]=d; Zglass[2][i]=id; ZglassO[2][i]=O;     \
       }                                                        \
       else if(d > Zglassd[2][i]+600.0 &&                       \
               d < Zglassd[3][i]-600.0){                        \
         Zglassd[3][i]=d; Zglass[3][i]=id; ZglassO[3][i]=O;     \
       }                                                        \
     }                                                          \
     else{                                                      \
       if(d <= Zglassd[0][i]){                                  \
         Zglass[0][i] = -1;  Zglassd[0][i]=FARAWAY;             \
         Zglass[1][i] = -1;  Zglassd[1][i]=FARAWAY;             \
         Zglass[2][i] = -1;  Zglassd[2][i]=FARAWAY;             \
         Zglass[3][i] = -1;  Zglassd[3][i]=FARAWAY;             \
       }                                                        \
       else if(d < Zglassd[1][i]){                              \
         Zglass[1][i] = -1;  Zglassd[1][i]=FARAWAY;             \
         Zglass[2][i] = -1;  Zglassd[2][i]=FARAWAY;             \
         Zglass[3][i] = -1;  Zglassd[3][i]=FARAWAY;             \
       }                                                        \
       else if(d < Zglassd[2][i]){                              \
         Zglass[2][i] = -1;  Zglassd[2][i]=FARAWAY;             \
         Zglass[3][i] = -1;  Zglassd[3][i]=FARAWAY;             \
       }                                                        \
       else if(d < Zglassd[3][i]){                              \
         Zglass[3][i] = -1;  Zglassd[3][i]=FARAWAY;             \
       }                                                        \
       Zdepth[i]=d; Zbuffer[i]=id; Zobject[i]=O;                \
     }                                                          \
}
#else
#error BADGLASS0
#endif

void UpdateGlassScanBuffers(double xl, double xr, long id, long O,  // used with NOT tracing refractions
                       face  *f, double d2,
                       long *Zglass[], long *ZglassO[], double *Zglassd[],
                       long *Zbuffer, long *Zobject,
                       double *Zdepth, double *Zposn){
 long i,ixl,ixr,mapID;
 double depth,bL,bR,aL,aR,dL,dR;
 unsigned char mattr;
 if((mapID=f->ffmap) >= 0 && mapID < Object[O].NoMaps &&  // this face has a transparency map
     Object[O].Mbase[mapID].tp > 0)mattr=255;             // and we need to flag the glass buffers
 else if(f->ffmat < 0 || f->ffmat >= Object[O].NoMats)mattr=0;
 else   mattr=(Object[O].Tbase + f->ffmat)->transp;
/* xl += 0.0001; too small
   xr -= 0.0001;
   xl += 0.01;  too big
   xr -= 0.01;  */
 xl += 0.001;
 xr -= 0.001;
 ixl = (long)xl;
 ixr = (long)xr;
 if(ixl == ixr){
   dL = xl - (double)Xcentre; dL/=scalex;
   dR = xr - (double)Xcentre; dR/=scalex;
   bL = f->n[0]*dL + f->n[2]*d2 + f->n[1];
   aL = (f->c)/bL;
   bR = f->n[0]*dR + f->n[2]*d2 + f->n[1];
   aR = (f->c)/bR;
   depth=min(aL,aR);
   if(depth > 1.000 && depth <= Zdepth[ixl]){
     UpdateGlassBuffer(f,ixl,depth,id,O,mattr)
     Zposn[ixl]=(xr+xl)/2.0;
   }
   return;
 }
 else{
   dL = xl - (double)Xcentre; dL/=scalex;
   bL = f->n[0]*dL + f->n[2]*d2 + f->n[1];
   aL = (f->c)/bL;
   bR = f->n[0]*DR[ixl] + f->n[2]*d2 + f->n[1];
   aR = (f->c)/bR;
   depth=min(aL,aR);
   if(depth > 1.000 && depth <= Zdepth[ixl]){
     UpdateGlassBuffer(f,ixl,depth,id,O,mattr)
     Zposn[ixl]=xl;
   }
   bL = f->n[0]*DL[ixr] + f->n[2]*d2 + f->n[1];
   aL = (f->c)/bL;
   dR = xr - (double)Xcentre; dR/=scalex;
   bR = f->n[0]*dR + f->n[2]*d2 + f->n[1];
   aR = (f->c)/bR;
   depth=min(aL,aR);
   if(depth > 1.000 && depth <= Zdepth[ixr]){
     UpdateGlassBuffer(f,ixr,depth,id,O,mattr)
     Zposn[ixr]=xr;
   }
 }
 if(ixr-ixl  > 1)for(i=ixl+1;i<ixr;i++)
 {
   bL = f->n[0]*DL[i] + f->n[2]*d2 + f->n[1];
   aL = (f->c)/bL;
   bR = f->n[0]*DR[i] + f->n[2]*d2 + f->n[1];
   aR = (f->c)/bR;
   depth=min(aL,aR);
   if(depth > 1.000 && depth <= Zdepth[i]){
     UpdateGlassBuffer(f,i,depth,id,O,mattr)
     Zposn[i]=(double)i+0.5;
   }
 }
}

void UpdateNoGlassScanBuffers(double xl, double xr, long id, long O,   // used with tracing refractions
                       face  *f, double d2,
                       long *Zglass[], long *ZglassO[], double *Zglassd[],
                       long *Zbuffer, long *Zobject,
                       double *Zdepth, double *Zposn){
 long i,ixl,ixr;
 double depth,bL,bR,aL,aR,dL,dR;
 unsigned char mattr;
 if(f->ffmat < 0 || f->ffmat >= Object[O].NoMats)mattr=0;
 else   mattr=(Object[O].Tbase + f->ffmat)->transp;
/* xl += 0.0001; too small
   xr -= 0.0001;
   xl += 0.01;  too big
   xr -= 0.01;  */
 xl += 0.001;
 xr -= 0.001;
 ixl = (long)xl;
 ixr = (long)xr;
 if(ixl == ixr){
   dL = xl - (double)Xcentre; dL/=scalex;
   dR = xr - (double)Xcentre; dR/=scalex;
   bL = f->n[0]*dL + f->n[2]*d2 + f->n[1];
   aL = (f->c)/bL;
   bR = f->n[0]*dR + f->n[2]*d2 + f->n[1];
   aR = (f->c)/bR;
   depth=min(aL,aR);
   if(depth > 1.000 && depth <= Zdepth[ixl]){
     UpdateNoGlassBuffer(f,ixl,depth,id,O)
     Zposn[ixl]=(xr+xl)/2.0;
   }
   return;
 }
 else{
   dL = xl - (double)Xcentre; dL/=scalex;
   bL = f->n[0]*dL + f->n[2]*d2 + f->n[1];
   aL = (f->c)/bL;
   bR = f->n[0]*DR[ixl] + f->n[2]*d2 + f->n[1];
   aR = (f->c)/bR;
   depth=min(aL,aR);
   if(depth > 1.000 && depth <= Zdepth[ixl]){
     UpdateNoGlassBuffer(f,ixl,depth,id,O)
     Zposn[ixl]=xl;
   }
   bL = f->n[0]*DL[ixr] + f->n[2]*d2 + f->n[1];
   aL = (f->c)/bL;
   dR = xr - (double)Xcentre; dR/=scalex;
   bR = f->n[0]*dR + f->n[2]*d2 + f->n[1];
   aR = (f->c)/bR;
   depth=min(aL,aR);
   if(depth > 1.000 && depth <= Zdepth[ixr]){
     UpdateNoGlassBuffer(f,ixr,depth,id,O)
     Zposn[ixr]=xr;
   }
 }
 if(ixr-ixl  > 1)for(i=ixl+1;i<ixr;i++)
 {
   bL = f->n[0]*DL[i] + f->n[2]*d2 + f->n[1];
   aL = (f->c)/bL;
   bR = f->n[0]*DR[i] + f->n[2]*d2 + f->n[1];
   aR = (f->c)/bR;
   depth=min(aL,aR);
   if(depth > 1.000 && depth <= Zdepth[i]){
     UpdateNoGlassBuffer(f,i,depth,id,O)
     Zposn[i]=(double)i+0.5;
   }
 }
}

long Intersect(long scanline, vertex  *v1, vertex  *v2,
                double *xi){
 double ddx,ddy,doff,ddxy;
 if((v1->y < (double)scanline)  && (v2->y < (double)scanline))return 0;
 if((v1->y > (double)scanline)  && (v2->y > (double)scanline))return 0;
 ddy = (v2->y) - (v1->y);
 if(fabs(ddy) < 0.0001)return 0;
 ddx = (v2->x) - (v1->x);
 doff=(double)scanline - v1->y;
 ddxy=ddx*doff/ddy;
 *xi=(v1->x)+ddxy;
 return 1;
}

long Active(long scanline, vertex  *v, face  *f,
                      double *leftedge, double *rightedge){
 long v1,v2,v3
     ,i1,i2,i3;
 double x1,x2,x3;
 if(f->top > scanline)       return 0;
 if(f->bottom < scanline)    return 0;
 v1=f->V[0];
 v2=f->V[1];
 v3=f->V[2];
 i1=Intersect(scanline,(v+v1),(v+v2),&x1);
 i2=Intersect(scanline,(v+v2),(v+v3),&x2);
 i3=Intersect(scanline,(v+v3),(v+v1),&x3);
 if((i1 == 1) && (i2 == 1) && (i3 == 1))
 {*leftedge = min(min(x1,x2),x3);
  *rightedge= max(max(x1,x2),x3);}
 else if((i2 == 1) && (i3 == 1))
 {*leftedge = min(x2,x3);
  *rightedge= max(x2,x3);}
 else if((i3 == 1) && (i1 == 1))
 {*leftedge = min(x3,x1);
  *rightedge= max(x3,x1);}
 else if((i1 == 1) && (i2 == 1))
 {*leftedge = min(x1,x2);
  *rightedge= max(x1,x2);}
 else return 0;
 if(*leftedge >= (double)XMAX-0.002)return 0;
 else *leftedge  = max(0,*leftedge);
 if(*rightedge < (double)XMIN+0.002)return 0;
 *rightedge = min((double)XMAX-0.01,*rightedge);
 return 1;
}

void SetPerspectiveView(long StereoLeft){  
 double t1[4][4],t2[4][4],t3[4][4],t4[4][4];
 double lenz=50.0;
 if(ResolutionX > ResolutionY)
   scalex=scaley=(double)ResolutionX/2.0/21.22*lenz;
 else
   scalex=scaley=(double)ResolutionY/2.0/21.22*lenz;
 if(ResolutionX <= 320 && ResolutionY <= 200)
   scaley /= 1.2;     /* 3ds uses  *0.82 = /1.21 but this seems wrong */
 if(anti_alias == LOW)   {scalex *= 2;}
 if(anti_alias == MEDIUM){scalex *= 2; scaley *= 2;}
 if(anti_alias == HIGH)  {scalex *= 3; scaley *= 2;}
 if(anti_alias == ULTRA) {scalex *= 3; scaley *= 3;}
 CamTheta *= PI/180.0;
 CamPhi   *= PI/180.0;
 CamAlpha *= PI/180.0;
 scalex *= CamSx;  /* apply any more scaling or aspect ration */
 scaley *= CamSy;
 R_tram(t1,-ViewPoint[0],-ViewPoint[1],-ViewPoint[2]);
 R_rotz(t2,-CamPhi);
 R_m4by4(t2,t1,t3);
 R_rotx(t1,-CamAlpha);
 R_m4by4(t1,t3,t2);
 R_roty(t1,-CamTheta);   /* bank on axis */
 R_m4by4(t1,t2,ViewTransform);
 R_roty(t1,CamTheta);    /* for reflection maps get a transform back to */
 R_rotx(t2,CamAlpha);    /* co-ordinate system without camera,          */
 R_m4by4(t2,t1,t3);      /* but this is easier to compute               */
 R_rotz(t1,CamPhi);
 R_m4by4(t1,t3,AntiCamera);
 R_rotz(t1,CamPhi);      /* for sky environment map                     */
 R_rotx(t2,CamAlpha);
 // R_m4by4(t1,t2,CameraRotate);  // use this on its own for NO bank    */
 R_roty(t3,CamTheta);
 R_m4by4(t2,t3,t4);
 R_m4by4(t1,t4,CameraRotate);  // Camera Rotation matrix
 // The AntiCamera matrix ordered to accomodate the OpenFX to OpenGL 
 // coordinate differences. This does not include the Transpose required
 // when writing to the OpenGL transformation matrix.
 if(StereoRenderingON == YES){
   double tStereo[4][4],rStereo[4][4];
   double r,d;
   r=atan2(CamStereoSeparation/2.0,CamParallaxDepth);
   d=CamStereoSeparation/2.0*CamRuler;
   if(StereoLeft == 1){
     R_tram(tStereo,d,0.0,0.0);
     R_rotz(rStereo,r); 
   }
   else{
     R_tram(tStereo,-d,0.0,0.0);
     R_rotz(rStereo,-r); 
   }
   R_m4by4(rStereo,ViewTransform,t1);
   R_m4by4(tStereo,t1,ViewTransform);
   R_m4by4(rStereo,AntiCamera,t1);
   R_c4to4(t1,AntiCamera);
   R_m4by4(rStereo,CameraRotate,t1);
   R_c4to4(t1,CameraRotate);
 }
 AntiCameraGL[0][0]= AntiCamera[0][0];
 AntiCameraGL[0][1]= AntiCamera[0][2];
 AntiCameraGL[0][2]=-AntiCamera[0][1];
 AntiCameraGL[1][0]= AntiCamera[2][0];
 AntiCameraGL[1][1]= AntiCamera[2][2];
 AntiCameraGL[1][2]=-AntiCamera[2][1];
 AntiCameraGL[2][0]=-AntiCamera[1][0];
 AntiCameraGL[2][1]=-AntiCamera[1][2];
 AntiCameraGL[2][2]= AntiCamera[1][1];
 return;
}

/*  perform geometry transformations   */
void ApplyViewingTransformToPoint(vector vin, vector vout){
 double x1,y1,z1;
 R_m4by1(ViewTransform,vin[0],vin[1],vin[2],&x1,&y1,&z1);
 vout[0]=x1;    vout[1]=y1;    vout[2]=z1;
}

void BendNormal(vector nin, vector dn, vector nout){
 double x1,y1,z1,x2,y2,z2;
 R_m3by1(AntiCamera,nin[0],nin[1],nin[2],&x1,&y1,&z1);
 x2=x1+dn[0]; y2=y1+dn[1]; z2=z1+dn[2];
 R_m3by1(ViewTransform,x2,y2,z2,&nout[0],&nout[1],&nout[2]);
}

void MakeObjectTransformation(
             double fi, double theta, double alpha, // orientation
             short im, double ima,                  // internal rotations
             double sx, double sy, double sz,       // scales
             vector Centre,                         // centre point offset  
             vector Position,                       // Object Position in World
             double trpos[4][4],     // output the object transform (include view
             double trset[4][4],     // exclude the view transform 
             double trinv[4][4]){    // undos the whole transform (include view)
 double t1[4][4],t2[4][4],t3[4][4],t4[4][4];
 fi    *= PI/180.0;
 theta *= PI/180.0;
 alpha *= PI/180.0;
 R_tram(t2,-(double)Centre[0],-(double)Centre[1],-(double)Centre[2]);
 R_scal(t3,sx,sy,sz);
 R_m4by4(t3,t2,t1);
 if(im > 0){
   if(im == 1)R_rotz(t3,ima*PI/180.0);
   if(im == 2)R_rotx(t3,ima*PI/180.0);
   if(im == 3)R_roty(t3,ima*PI/180.0);
   R_m4by4(t3,t1,t2);
   R_roty(t1,theta);
   R_m4by4(t1,t2,t3);
 }
 else{
   R_roty(t2,theta);
   R_m4by4(t2,t1,t3);
 }
 R_rotx(t1,alpha);
 R_m4by4(t1,t3,t2);
 R_rotz(t1,fi);
 R_m4by4(t1,t2,t3);
 R_tram(t1,Position[0],Position[1],Position[2]);
 R_m4by4(t1,t3,trset);                 // excluding the viewing transform
 R_m4by4(ViewTransform,trset,trpos);   // including the viewing transform

 /* Get the inverse object transformation */

 R_tram(t1,-Position[0],-Position[1],-Position[2]);
 R_rotz(t2,-fi);
 R_m4by4(t2,t1,t3);
 R_rotx(t1,-alpha);
 R_m4by4(t1,t3,t2);
 if(im > 0){
   R_roty(t1,-theta);
   if(im == 1)R_rotz(t3,-ima*PI/180.0);
   if(im == 2)R_rotx(t3,-ima*PI/180.0);
   if(im == 3)R_roty(t3,-ima*PI/180.0);
   R_m4by4(t3,t1,t4);
 }
 else{
   R_roty(t4,-theta);
 }
 R_m4by4(t4,t2,t1);
 R_scal(t2,1.0/sx,1.0/sy,1.0/sz);
 R_m4by4(t2,t1,t3);
 R_tram(t2,(double)Centre[0],(double)Centre[1],(double)Centre[2]);
 R_m4by4(t2,t3,t1);
 R_m4by4(t1,AntiCamera,trinv);
 return;
}

void TransformMappingRectangle(
                    double  trpos[4][4], vector  P, vector  X, vector Y,
                    vector  p, vector  x, vector y, vector n){
  R_m4by1(trpos,P[0],P[1],P[2],&p[0],&p[1],&p[2]);
  R_m4by1(trpos,X[0],X[1],X[2],&x[0],&x[1],&x[2]);
  R_m4by1(trpos,Y[0],Y[1],Y[2],&y[0],&y[1],&y[2]);
  vecsub(x,p,x);
  vecsub(y,p,y);
  cross (x,y,n);
  if(normalize(n)){ /* bad axis just set up unit */
    x[0]=1.0; x[1]=0.0; x[2]=0.0;
    y[0]=0.0; y[1]=1.0; y[2]=0.0;
    cross(x,y,n);
  }
}


void R_rotz(double tr[4][4], double ang){
 long i,j;
 for(i=0;i<4;i++)
 for(j=0;j<4;j++)
 {
  tr[i][j]=0.0;
  if(i == j)tr[i][j]=1.0;
 }
  tr[0][0] =  cos(ang);
  tr[0][1] = -sin(ang);
  tr[1][0] =  sin(ang);
  tr[1][1] =  cos(ang);
  return;
}

#ifdef _SUNSTYLE
void R_rotx(tr,ang) double tr[4][4]; double ang; {
#else
void R_rotx(double tr[4][4], double ang){
#endif
 long i,j;
 for(i=0;i<4;i++)
 for(j=0;j<4;j++)
 {
  tr[i][j]=0.0;
  if(i == j)tr[i][j]=1.0;
 }
  tr[1][1] =  cos(ang);
  tr[1][2] = -sin(ang);
  tr[2][1] =  sin(ang);
  tr[2][2] =  cos(ang);
  return;
}

void R_roty(double tr[4][4], double ang){
 long i,j;
 for(i=0;i<4;i++)
 for(j=0;j<4;j++)
 {
  tr[i][j]=0.0;
  if(i == j)tr[i][j]=1.0;
 }
  tr[0][0] =  cos(ang);
  tr[0][2] = -sin(ang);
  tr[2][0] =  sin(ang);
  tr[2][2] =  cos(ang);
  return;
}

void R_tram(double t[4][4], double dx, double dy, double dz){
 long i,j;
 for(i=0;i<4;i++)
 for(j=0;j<4;j++)
 {
  t[i][j]=0.0;
  if(i == j)t[i][j]=1.0;
 }
 t[0][3]=dx;
 t[1][3]=dy;
 t[2][3]=dz;
 t[3][3]=1;
 return;
}

void R_scal(double t[4][4], double sx, double sy, double sz){
 long i,j;
 for(i=0;i<4;i++)
 for(j=0;j<4;j++)
 {
  t[i][j]=0.0;
  if(i == j)t[i][j]=1.0;
 }
 t[0][0]=sx;
 t[1][1]=sy;
 t[2][2]=sz;
 t[3][3]=1;
 return;
}

void R_m4by4(double t1[4][4], double t2[4][4], double tr[4][4]){
 long i,j,k;
 for(i=0;i<4;i++)
 for(j=0;j<4;j++)
 {
  tr[i][j]=0.0;
  for(k=0;k<4;k++)tr[i][j]=tr[i][j]+t1[i][k]*t2[k][j];
 }
 return;
}

void R_c4to4(double tin[4][4], double tout[4][4]){
 int i,j;
 for(i=0;i<4;i++)
 for(j=0;j<4;j++)
 tout[i][j]=tin[i][j];
 return;
}

void R_null_transform(double t[4][4]){
 short i,j;
 for(i=0;i<4;i++)
 for(j=0;j<4;j++)
 if(i == j)t[i][j]=1.0;
 else      t[i][j]=0.0;
}

void R_m4by1(double t4[4][4], double x, double y, double z,
           double *xx, double *yy, double *zz){
   *xx=t4[0][0]*x+t4[0][1]*y+t4[0][2]*z+t4[0][3];
   *yy=t4[1][0]*x+t4[1][1]*y+t4[1][2]*z+t4[1][3];
   *zz=t4[2][0]*x+t4[2][1]*y+t4[2][2]*z+t4[2][3];
 return;
}

void R_m3by1(double t4[4][4], double x, double y, double z,
           double *xx, double *yy, double *zz){
   *xx=t4[0][0]*x+t4[0][1]*y+t4[0][2]*z;
   *yy=t4[1][0]*x+t4[1][1]*y+t4[1][2]*z;
   *zz=t4[2][0]*x+t4[2][1]*y+t4[2][2]*z;
 return;
}


#ifdef _SUNSTYLE
void R_reflect(t,p,n) double t[4][4]; vector p; vector n; {
#else
void R_reflect(double t[4][4], vector p, vector n){
#endif
 /* p is point in reflection plane  */
 /* n is normal to reflection plane */
 long i,j;
 double dxy,phi,theta;
 double t1[4][4],t2[4][4],t3[4][4];
//{
//char temp[256];
//sprintf(temp,"n= %lf %lf %lf p=%lf %lf %lf",n[0],n[1],n[2],p[0],p[1],p[2]);
//MessageBox(NULL,temp,"Debug",MB_OK);
//}
 dxy=n[0]*n[0]+n[1]*n[1];
 R_tram(t1,-p[0],-p[1],-p[2]);
 if(fabs(n[2]) > 0.7071){
//MessageBox(NULL,"transform to Z","Debug",MB_OK);
  if(fabs(n[2]) > 0.998){ /* vertical so just reflect in XY plane and return */
// MessageBox(NULL,"Vertical","Debug",MB_OK);
     for(i=0;i<4;i++)for(j=0;j<4;j++){
       if(i == j)t2[i][j]=1.0;
       else      t2[i][j]=0.0;
     }
     t2[2][2]=-1.0;
     R_m4by4(t2,t1,t3);
     R_tram(t1,p[0],p[1],p[2]);
     R_m4by4(t1,t3,t);
     return;
   }
   dxy=sqrt(dxy);
   if     (fabs(n[0]) < 1.e-6 && n[1] > 0.0)phi =  1.570796;
   else if(fabs(n[0]) < 1.e-6 && n[1] < 0.0)phi = -1.570796;
   else phi = atan2(n[1],n[0]);
   if(n[2] > 0.0)theta =  1.570796-atan2(n[2],dxy);
   else          theta = -1.570796-atan2(n[2],dxy);
   R_rotz(t2, -phi);
   R_m4by4(t2,t1,t3);
   R_roty(t2, theta);
   R_m4by4(t2,t3,t1);
   for(i=0;i<4;i++)for(j=0;j<4;j++){
     if(i == j)t2[i][j]=1.0;
     else      t2[i][j]=0.0;
   }
   t2[2][2]=-1.0;  /* reflect in XY plane - perpend`r to Z axis */
   R_m4by4(t2,t1,t3);
   R_roty(t2,-theta);
   R_m4by4(t2,t3,t1);
   R_rotz(t2,phi);
   R_m4by4(t2,t1,t3);
   R_tram(t1,p[0],p[1],p[2]);
   R_m4by4(t1,t3,t);
 }
 else{
   dxy=sqrt(dxy);
   if     (fabs(n[0]) < 1.e-6 && n[1] > 0.0)phi =  1.570796;
   else if(fabs(n[0]) < 1.e-6 && n[1] < 0.0)phi = -1.570796;
   else phi = atan2(n[1],n[0]);
   theta = atan2(n[2],dxy);
   R_rotz(t2, -phi);
   R_m4by4(t2,t1,t3);
   R_roty(t2, -theta);
   R_m4by4(t2,t3,t1);
   for(i=0;i<4;i++)for(j=0;j<4;j++){
     if(i == j)t2[i][j]=1.0;
     else      t2[i][j]=0.0;
   }
   t2[0][0]=-1.0;  /* reflect in YZ plane - perpend`r to X axis */
   R_m4by4(t2,t1,t3);
   R_roty(t2,theta);
   R_m4by4(t2,t3,t1);
   R_rotz(t2,phi);
   R_m4by4(t2,t1,t3);
   R_tram(t1,p[0],p[1],p[2]);
   R_m4by4(t1,t3,t);
 }
 return;
}

void veccopy(vector v1, vector r){
 r[0] = v1[0];
 r[1] = v1[1];
 r[2] = v1[2];
}

void vecsub(vector v1, vector v2, vector r){
 r[0] = v1[0] - v2[0];
 r[1] = v1[1] - v2[1];
 r[2] = v1[2] - v2[2];
}

void vecsum(vector v1,vector v2, vector r){
 r[0] = v1[0] + v2[0];
 r[1] = v1[1] + v2[1];
 r[2] = v1[2] + v2[2];
}

void vecscale(double s, vector v, vector r){
 r[0] = s * v[0];
 r[1] = s * v[1];
 r[2] = s * v[2];
}

short normalize(vector v){
 double n,nn;
 n= (double)((v[0]*v[0]) + (v[1]*v[1]) + (v[2]*v[2]));
 if(n < 1.e-10)return 1;
 nn=sqrt(n);
 n = 1.0 / nn;
 vecscale(n,v,v);
 return 0;
}

double R_length(vector v){
 double n,nn;
 n= (double)((v[0]*v[0]) + (v[1]*v[1]) + (v[2]*v[2]));
 nn=sqrt(n);
 return (nn);
}

void cross(vector v1,vector v2, vector r){
 r[0] = (v1[1]*v2[2]) - (v2[1]*v1[2]);
 r[1] = (v1[2]*v2[0]) - (v1[0]*v2[2]);
 r[2] = (v1[0]*v2[1]) - (v2[0]*v1[1]);
}

double dot(vector v1, vector v2){
 double result;
 result  = v1[0] * v2[0];
 result += v1[1] * v2[1];
 result += v1[2] * v2[2];
 return result;
}

void AssignMappingCoordinatesToVertices(object *O){
 vertex *Vbase=O->Vbase;
 face   *f,*Fbase=O->Fbase;
 imap   *m,*Mbase=O->Mbase;
 matl   *Tbase=O->Tbase;
 axis   *a;
 long nfaces=O->NoFaces;
 long nvert=O->NoVertices;
 long nmaps=O->NoMaps;
 long nmats=O->NoMats;
 long i,j,mapid,matid;
 double alpha,beta,gamma,alphav[3],betav[3];
 if(nmaps > 0)for(i=0,m=Mbase;i<nmaps;i++,m++){
 //  m->map=m->mapin;
   vecsub(m->xin,m->pin,m->x);
   vecsub(m->yin,m->pin,m->y);
   cross(m->x,m->y,m->n);
   if(normalize(m->n)){ /* bad axis just set up unit */
     m->x[0]=1.0; m->x[1]=0.0; m->x[2]=0.0;
     m->y[0]=0.0; m->y[1]=1.0; m->y[2]=0.0;
     cross(m->x,m->y,m->n);
   }
 }
 for(i = -1; i<nmats; i++){ // for axes 
   if(i == -1)a = &(O->Abase);        // default axis for object
   else       a = &((Tbase+i)->ax);   // material axes
   vecsub(a->buin,a->boin,a->bu);
   vecsub(a->bvin,a->boin,a->bv);
   cross(a->bu,a->bv,a->bn);
   normalize(a->bn);
   a->tl=max(1.0,R_length(a->bu));
   a->tv=max(1.0,R_length(a->bv));
   vecscale(1.0/a->tl,a->bu,a->bu_norm);
   vecscale(1.0/a->tv,a->bv,a->bv_norm);
 }
 for(i=0,f=Fbase;i<nfaces;i++,f++){
   if(f->bFmc){;} // face has mapping coords from the designer - do nothing
   else if(f->ffmap >= 0){ // if mapped, assign surface coordinates from the map
     mapid=f->ffmap;
     if(mapid < nmaps){
       if((m=(Mbase+mapid))->map != MAP_BY_VERTEX){
         for(j=0;j<3;j++){
           if(m->map == CYLINDER){
              hitcylinder(m->n,m->x,m->y,m->p,
                    (Vbase+(f->V[j]))->p,&(alphav[j]),&(betav[j]));
           }
           else if(m->map == MAP_SPHERICAL){
             hitsphere(m->n,m->x,m->y,m->p
                    ,(Vbase+(f->V[j]))->p,&(alphav[j]),&(betav[j]));
           }
           else    /* assume planar map */
             hitpoint(m->n,m->x,m->y,m->p
                    ,(Vbase+(f->V[j]))->p,&(alphav[j]),&(betav[j]));
         }  
         if(m->map == CYLINDER || m->map == MAP_SPHERICAL)FixUpMapEdge(alphav);
         for(j=0;j<3;j++){
           if(m->map == CYLINDER)alphav[j] *= m->angle; // scale for multiple rotations
           f->uu[j]=alphav[j];
           f->vv[j]=betav[j];
           f->ww[j]=0.5;
          //if(debug != NULL)fprintf(debug,"Calculating  map coords for face %ld %ld (%lf %lf)\n",i,j,f->uu[j],f->vv[j]);
         }
       }
     }
   }
   else{ // if not mapped, assign texture coordinated from the shader axis
     for(j=0;j<3;j++){
       matid=f->ffmat;
       if(matid >= 0 && matid < nmats)a = &((Tbase+matid)->ax);
       else                           a = &(O->Abase);
       gamma=hitpoint(a->bn,
                      a->bu,
                      a->bv,
                      a->bo,
                     (Vbase+(f->V[j]))->p,&alpha,&beta);
       gamma /= max(1.0,a->tl);
       f->uu[j]=alpha;
       f->vv[j]=beta;
       f->ww[j]=gamma;
     }    
   }
 }
 // Reset to use mapping coords for all faces
 if(nmaps > 0)for(i=0,m=Mbase;i<nmaps;i++,m++){
   m->map = MAP_BY_VERTEX;
 }
}

#define MAPDELTA 0.33333

static void FixUpMapEdge(GLfloat a[]){
 if(a[0] < MAPDELTA){if(a[1] >= 1.0-MAPDELTA)a[1] -= 1.0; if(a[2] >= 1.0-MAPDELTA)a[2] -= 1.0; return;}
 if(a[1] < MAPDELTA){if(a[2] >= 1.0-MAPDELTA)a[2] -= 1.0; if(a[0] >= 1.0-MAPDELTA)a[0] -= 1.0; return;}
 if(a[2] < MAPDELTA){if(a[0] >= 1.0-MAPDELTA)a[0] -= 1.0; if(a[1] >= 1.0-MAPDELTA)a[1] -= 1.0; return;}
}