MATERIAL.C

Go to the documentation of this file.

// file material.c  

#define MODULE_MATERIAL

/* Calculate the surface colour for given Pixel now that we know which     */
/* face it is on.  Intepolate the PHONG vector first. The ground is an     */
/* independent entity so that its caclulations can be optimised.           */

/* most shaders have a function to caclulate their value                   */

#include "render.h"

extern short mirror_ground;
extern double Rip_Freq;

#define FLOOR(x) ((x) >= 0.0 ? floor(x) : (0.0 - floor(0.0 - (x)) - 1.0))

static void Transform_to_Moziac(double, double ,double *, double *);
static void AttenuateGlassSurface(long , long , long ,
                                  UCHAR *, UCHAR *, UCHAR *,
                                  long *[], long *[], double *[], double *);
static void AttenuateGlass(long, long,
                    UCHAR *, UCHAR *, UCHAR *Blue,
                    long *[], long *[], double *[], double *);
static void PutBackPlanes(long i, long scanline, double d2,
                   UCHAR *Red, UCHAR *Green, UCHAR *Blue,
                   double *Zground, long *Zobject, double *Zdepth,
                   UCHAR *ZobjChar, long Wsky, double Rsky);
static void GetBackdropValue(long row, long col, vector p, vector d,
                             double *colour);

double spec_pow[16][max_spec+1],spec_gro[max_spec+1];

static double spec_val[16]={
       10.0, 15.0, 20.0, 25.0, 27.0, 30.0, 32.0, 35.0,
       40.0, 45.0, 50.0, 55.0, 60.0, 70.0, 80.0, 99.0};

void init_spec(void){
 long i,j; /* use LUT for specular highlight function to speed up rendering */
 double d;
 for(j=0;j<16;j++){
   spec_pow[j][max_spec] = 1.0;
   for(i=0;i<max_spec;i++){
     d = pow(((double)i/(double)max_spec),spec_val[j]);
     spec_pow[j][i]=d;
   }
 }
 for(i=0;i<max_spec;i++)spec_gro[i] = pow((1.0/max_spec*(double)i),200.0);
}


/* The functions for Bump and Reflection Mapping  */

/* Map->bumpX/Y is a scale factor to scale the change in normal            */
/* vector n is unity, max gradient is 255 -> scale factor is               */
/* 1/|x| * max_displ_of normal / 255 * height /100                         */
/* max_displ = tan(a) where a is the max angle the normal is bent          */
/* For Cylinder maps the X axis is defined by n X y which is normalized    */
/* so the appropriate scaling factor is different                          */

void BumpByCoordinates(imap *Map, short type, long obj, face *f, vector p, vector n){
 long  i,j,ialpha,ibeta,k,kl,kr,kt,kb,pixel,ix,iy;
 double alpha0,beta0,alpha,beta,a,b,red[3],cc0,cl0,cr0,ct0,cb0,
        ve1,ve2,ve3,vp1,vp2,vp3,p1,p2,p3,
        detmax,dm,det1,det2,det3,
        t1x,t2x,t1y,t2y;
 vertex   *V,*v0,*v1,*v2;
 static double dalpha_bump[3]={0.0,0.01,0.0},
                dbeta_bump[3]={0.0,0.0,0.01};
 vector vv1,vv2;
 unsigned char *pP,*pR;
 if(Map->bp < 0)return;
 j=Map->xmaxbump;
 if((V=Object[obj].Vbase) == NULL)return;
 v0=(V+(f->V[0]));
 v1=(V+(f->V[1]));
 v2=(V+(f->V[2]));
 ve1=v1->p[0] - v0->p[0]; ve2=v1->p[1] - v0->p[1]; ve3=v1->p[2] - v0->p[2];
 vp1=v2->p[0] - v0->p[0]; vp2=v2->p[1] - v0->p[1]; vp3=v2->p[2] - v0->p[2];
 p1=p[0]-v0->p[0];
 p2=p[1]-v0->p[1];
 p3=p[2]-v0->p[2];
 det1=ve1*vp2-vp1*ve2;
 det2=ve1*vp3-vp1*ve3;
 det3=ve2*vp3-vp2*ve3;
 k=0; detmax=1.e-3;  /* same as TOL in mirrors */
 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 == 0)return;
 else 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;
 }
 // Calculate UV map address of pixel
 t1x=f->uu[1] - f->uu[0]; t2x=f->uu[2] - f->uu[0];
 t1y=f->vv[1] - f->vv[0]; t2y=f->vv[2] - f->vv[0];
 alpha0 = f->uu[0] + a*t1x + b*t2x;
 beta0  = f->vv[0] + a*t1y + b*t2y;

 for(i=0;i<3;i++){
   alpha = alpha0+dalpha_bump[i];
   beta  = beta0+dbeta_bump[i];
   if(type == TILE){
     if(Map->moziac)Transform_to_Moziac(alpha,beta,&a,&b);
     else{
       a=(alpha-floor(alpha));
       b=(beta-floor(beta));
     }
     alpha=a*Map->xmaxbump;
     beta=b*Map->ymaxbump;
   }
   else{
     alpha=Map->xmaxbump*alpha;
     beta =Map->ymaxbump*beta;
   }
   ialpha = (long)alpha; alpha -= (double)ialpha;
   ibeta  = (long)beta;  beta  -= (double)ibeta;
   if(ialpha < 0 || ialpha >= Map->xmaxbump)return;
   if(ibeta  < 0 || ibeta  >= Map->ymaxbump)return;
   k =ibeta*j+ialpha;                       /* do the bi-linear interpolation */
   if(type == TILE){
     if(Map->moziac){
       if(ialpha <   1)kl=1;                     else kl=k-1;
       if(ialpha > j-2)kr=k-1;                   else kr=k+1;
       if(ibeta  <   1)kb=k+j;                   else kb=k-j;
       if(ibeta  > Map->ymaxbump-2)kt=k-j; else kt=k+j;
     }
     else{
       if(ialpha <   1)kl=k+j-1;                            else kl=k-1;
       if(ialpha > j-2)kr=ibeta*j;                          else kr=k+1;
       if(ibeta  <   1)kb=(Map->ymaxbump-1)*j+ialpha; else kb=k-j;
       if(ibeta  > Map->ymaxbump-2)kt=ialpha;         else kt=k+j;
     }
   }
   else{
     if(ialpha <   1)kl=k; else kl=k-1;
     if(ialpha > j-2)kr=k; else kr=k+1;
     if(ibeta  <   1)kb=k; else kb=k-j;
     if(ibeta  > Map->ymaxbump-2)kt=k; else kt=k+j;
   }

   pR=Map->p24Gbump;
   cc0 = (double)pR[k];
   cl0 = (double)pR[kl];
   cr0 = (double)pR[kr];
   ct0 = (double)pR[kt];
   cb0 = (double)pR[kb];
   if(alpha > 0.5){
     alpha = alpha-0.5;
     if(beta > 0.5){
       beta = beta-0.5;
       red[i] = (cr0-cc0)*alpha+(ct0-cc0)*beta+cc0;
     }
     else{
       beta = 0.5-beta;
       red[i] = (cr0-cc0)*alpha+(cb0-cc0)*beta+cc0;
     }
   }
   else{
     alpha = 0.5-alpha;
     if(beta > 0.5){
       beta = beta-0.5;
       red[i] = (cl0-cc0)*alpha+(ct0-cc0)*beta+cc0;
     }
     else{
       beta = 0.5-beta;
       red[i] = (cl0-cc0)*alpha+(cb0-cc0)*beta+cc0;
     }
   }
 }
 vecscale((red[1]-red[0])*Map->bumpx,Map->x,vv1);
 vecscale((red[2]-red[0])*Map->bumpy,Map->y,vv2);
 vecsum(n,vv1,n);
 vecsum(n,vv2,n);
 normalize(n);
 return;
}

void BoxMapIntersect(vector d, long *id, double *alpha, double *beta){
 // id 0=front 1=right 2=back 3=left 4=top 5=bottom
 double a,b,t=0.1;
 long i = -1;   // no hit
 normalize(d);
 if(d[1] >  t){ // try front
   a = d[0]/d[1]; b = d[2]/d[1];
   if(-1.0 <= a && a <= 1.0 && -1.0 <= b && b <= 1.0){i=0; goto HIT;}
 }
 if(d[0] >  t){ // try right
   a = d[1]/d[0]; b = d[2]/d[0]; a = -a;
   if(-1.0 <= a && a <= 1.0 && -1.0 <= b && b <= 1.0){i=1; goto HIT;}
 }
 if(d[1] < -t){ // try back
   a = -d[0]/d[1]; b = -d[2]/d[1]; a = -a;
   if(-1.0 <= a && a <= 1.0 && -1.0 <= b && b <= 1.0){i=2; goto HIT;}
 }
 if(d[0] < -t){ // try left
   a = -d[1]/d[0]; b = -d[2]/d[0];
   if(-1.0 <= a && a <= 1.0 && -1.0 <= b && b <= 1.0){i=3; goto HIT;}
 }
 if(d[2] >  t){ // try top
   a = d[0]/d[2]; b = d[1]/d[2];  b= -b;
   if(-1.0 <= a && a <= 1.0 && -1.0 <= b && b <= 1.0){i=4; goto HIT;}
 }
 if(d[2] < -t){ // try bottom
   a = -d[0]/d[2]; b = -d[1]/d[2];
   if(-1.0 <= a && a <= 1.0 && -1.0 <= b && b <= 1.0){i=5; goto HIT;}
 }
 HIT:
 *id = i;
 *alpha = 0.5*(a+1);
 *beta = 1.0 - 0.5*(b+1);
 return;
}

void MapMirror(double *R,double *G, double *B, imap *Map, vector p, vector n){
 /* p is input direction n is normal  d is reflection direction */
 /* map 4 images round the sphere and 2 up down                 */
 /* atan2 returns in range -pi to pi                            */
#define FOURDIV2PI 0.318309886     /* 2 images in 360 rotation */
#define TWODIVPI   0.318309886     /* 1 image in -90 to +90    */
 
 long j,pixel,k,ibeta,ialpha,kl,kr,kt,kb;
 double pdotn,alpha,beta,d12;
 double cc0,cc1,cc2,cl0,cl1,cl2,cr0,cr1,cr2,ct0,ct1,ct2,cb0,cb1,cb2;
 unsigned char *pP,*pR,*pG,*pB;
 vector d,dd;
 if(Map->rp < 0)return;
 pdotn=dot(p,n);
 vecscale(2*pdotn,n,d);
 vecsub(p,d,dd);   /* un camera the camera transform */
 R_m4by1(AntiCamera,dd[0],dd[1],dd[2],&d[0],&d[1],&d[2]);

 BoxMapIntersect(d,&j,&alpha,&beta); if(j < 0)return;

 alpha *= (double)Map->xmaxrefl;
 beta *= (double)Map->ymaxrefl;
 ialpha = (long)alpha; alpha -= (double)ialpha;
 ibeta=(long)beta; beta -= (double)ibeta;
 if(ialpha < 0 || ialpha >= Map->xmaxrefl)return;
 if(ibeta  < 0 || ibeta  >= Map->ymaxrefl)return;
 j=Map->xmaxrefl;
 k =ibeta*j+ialpha;

 if(ialpha <   1)kl=1;                 else kl=k-1;
 if(ialpha > j-2)kr=k-1;               else kr=k+1;
 if(ibeta  <   1)kb=k+j;               else kb=k-j;
 if(ibeta  > Map->ymaxrefl-2)kt=k-j;   else kt=k+j;

 pR=Map->p24Rrefl; pG=Map->p24Grefl; pB=Map->p24Brefl;
 cc0 = (double)pR[k];  cc1 = (double)pG[k];  cc2 = (double)pB[k];
 cl0 = (double)pR[kl]; cl1 = (double)pG[kl]; cl2 = (double)pB[kl];
 cr0 = (double)pR[kr]; cr1 = (double)pG[kr]; cr2 = (double)pB[kr];
 ct0 = (double)pR[kt]; ct1 = (double)pG[kt]; ct2 = (double)pB[kt];
 cb0 = (double)pR[kb]; cb1 = (double)pG[kb]; cb2 = (double)pB[kb];
 if(alpha > 0.5){
   alpha = alpha-0.5;
   if(beta > 0.5){
     beta = beta-0.5;
     *R   = (cr0-cc0)*alpha+(ct0-cc0)*beta+cc0;
     *G   = (cr1-cc1)*alpha+(ct1-cc1)*beta+cc1;
     *B   = (cr2-cc2)*alpha+(ct2-cc2)*beta+cc2;
   }
   else{
     beta = 0.5-beta;
     *R   = (cr0-cc0)*alpha+(cb0-cc0)*beta+cc0;
     *G   = (cr1-cc1)*alpha+(cb1-cc1)*beta+cc1;
     *B   = (cr2-cc2)*alpha+(cb2-cc2)*beta+cc2;
   }
 }
 else{
   alpha = 0.5-alpha;
   if(beta > 0.5){
     beta = beta-0.5;
     *R   = (cl0-cc0)*alpha+(ct0-cc0)*beta+cc0;
     *G   = (cl1-cc1)*alpha+(ct1-cc1)*beta+cc1;
     *B   = (cl2-cc2)*alpha+(ct2-cc2)*beta+cc2;
   }
   else{
     beta = 0.5-beta;
     *R   = (cl0-cc0)*alpha+(cb0-cc0)*beta+cc0;
     *G   = (cl1-cc1)*alpha+(cb1-cc1)*beta+cc1;
     *B   = (cl2-cc2)*alpha+(cb2-cc2)*beta+cc2;
   }
 }
}

void MapEnvironment(long col, long row, UCHAR *r, UCHAR *g, UCHAR *b, imap *Map){
 long j,pixel,k,ibeta,ialpha,kl,kr,kt,kb,w,kmin,kmax;
 double alpha,beta,R,G,B;
 double cc0,cc1,cc2,cl0,cl1,cl2,cr0,cr1,cr2,ct0,ct1,ct2,cb0,cb1,cb2;
 unsigned char *pP,*pR,*pG,*pB;
 vector d,dd;
 if(Map->pp < 0)return;
 d[0]=((double)col-(double)Xcentre)/scalex;
 d[1]=1.0;
 d[2]=((double)Ycentre-(double)row)/scaley;
 R_m3by1(CameraRotate,d[0],d[1],d[2],&dd[0],&dd[1],&dd[2]);
 BoxMapIntersect(dd,&j,&alpha,&beta); if(j < 0){
   *r = 0; *g=255; *b=0;  // for diagnosis
   return;
 }
 w=Map->xmax/6;   // 6 horizontal pictures assumed
 alpha = ((double)j+alpha)*(double)w;
 beta *= (double)Map->ymax;
 ialpha = (long)alpha; alpha -= (double)ialpha;
 ibeta=(long)beta; beta -= (double)ibeta;
 if(ialpha < 0)ialpha=0;
 if(ialpha >= Map->xmax)ialpha=Map->xmax-1;
 if(ibeta  < 0)ibeta=0;
 if(ibeta  >= Map->ymax)ibeta=Map->ymax-1;
 k = ibeta*Map->xmax+ialpha;
 kmin = j*w + 1;
 kmax = (j+1)*w - 2;
 if(ialpha <   kmin)kl=k;       else kl=k-1;
 if(ialpha > kmax)  kr=k;       else kr=k+1;
 if(ibeta  <   1)kb=k;          else kb=k-Map->xmax;
 if(ibeta  > Map->ymax-2)kt=k;  else kt=k+Map->xmax;
 pR=Map->p24R; pG=Map->p24G; pB=Map->p24B;
 // *r = pR[k];  *g = pG[k];  *b = pB[k]; return; // diagnostic
 cc0 = (double)pR[k];  cc1 = (double)pG[k];  cc2 = (double)pB[k];
 cl0 = (double)pR[kl]; cl1 = (double)pG[kl]; cl2 = (double)pB[kl];
 cr0 = (double)pR[kr]; cr1 = (double)pG[kr]; cr2 = (double)pB[kr];
 ct0 = (double)pR[kt]; ct1 = (double)pG[kt]; ct2 = (double)pB[kt];
 cb0 = (double)pR[kb]; cb1 = (double)pG[kb]; cb2 = (double)pB[kb];
 if(alpha > 0.5){
   alpha = alpha-0.5;
   if(beta > 0.5){
     beta = beta-0.5;
     R   = (cr0-cc0)*alpha+(ct0-cc0)*beta+cc0;
     G   = (cr1-cc1)*alpha+(ct1-cc1)*beta+cc1;
     B   = (cr2-cc2)*alpha+(ct2-cc2)*beta+cc2;
   }
   else{
     beta = 0.5-beta;
     R   = (cr0-cc0)*alpha+(cb0-cc0)*beta+cc0;
     G   = (cr1-cc1)*alpha+(cb1-cc1)*beta+cc1;
     B   = (cr2-cc2)*alpha+(cb2-cc2)*beta+cc2;
   }
 }
 else{
   alpha = 0.5-alpha;
   if(beta > 0.5){
     beta = beta-0.5;
     R   = (cl0-cc0)*alpha+(ct0-cc0)*beta+cc0;
     G   = (cl1-cc1)*alpha+(ct1-cc1)*beta+cc1;
     B   = (cl2-cc2)*alpha+(ct2-cc2)*beta+cc2;
   }
   else{
     beta = 0.5-beta;
     R   = (cl0-cc0)*alpha+(cb0-cc0)*beta+cc0;
     G   = (cl1-cc1)*alpha+(cb1-cc1)*beta+cc1;
     B   = (cl2-cc2)*alpha+(cb2-cc2)*beta+cc2;
   }
 }
 *r = (unsigned char)R;
 *g = (unsigned char)G;
 *b = (unsigned char)B;
 return;
}

void MapEnvironmentInDirection(vector dd, double *r, double *g, double *b, imap *Map){
 // this is like a reflection map
 long j,pixel,k,ibeta,ialpha,kl,kr,kt,kb,w,kmin,kmax;
 double alpha,beta,R,G,B;
 double cc0,cc1,cc2,cl0,cl1,cl2,cr0,cr1,cr2,ct0,ct1,ct2,cb0,cb1,cb2;
 unsigned char *pP,*pR,*pG,*pB;
 vector d;
 R_m4by1(AntiCamera,dd[0],dd[1],dd[2],&d[0],&d[1],&d[2]);
 if(Map->pp < 0)return; // Map inactive
 BoxMapIntersect(d,&j,&alpha,&beta); if(j < 0){
   *r = 0; *g=255; *b=0;  // for diagnosis
   return;
 }
 w=Map->xmax/6;   // 6 horizontal pictures assumed
 alpha = ((double)j+alpha)*(double)w;
 beta *= (double)Map->ymax;
 ialpha = (long)alpha; alpha -= (double)ialpha;
 ibeta=(long)beta; beta -= (double)ibeta;
 if(ialpha < 0)ialpha=0;
 if(ialpha >= Map->xmax)ialpha=Map->xmax-1;
 if(ibeta  < 0)ibeta=0;
 if(ibeta  >= Map->ymax)ibeta=Map->ymax-1;
 k = ibeta*Map->xmax+ialpha;
 kmin = j*w + 1;
 kmax = (j+1)*w - 2;
 if(ialpha <   kmin)kl=k;       else kl=k-1;
 if(ialpha > kmax)  kr=k;       else kr=k+1;
 if(ibeta  <   1)kb=k;          else kb=k-Map->xmax;
 if(ibeta  > Map->ymax-2)kt=k;  else kt=k+Map->xmax;
 pR=Map->p24R; pG=Map->p24G; pB=Map->p24B;
 // *r = pR[k];  *g = pG[k];  *b = pB[k]; return; // diagnostic
 cc0 = (double)pR[k];  cc1 = (double)pG[k];  cc2 = (double)pB[k];
 cl0 = (double)pR[kl]; cl1 = (double)pG[kl]; cl2 = (double)pB[kl];
 cr0 = (double)pR[kr]; cr1 = (double)pG[kr]; cr2 = (double)pB[kr];
 ct0 = (double)pR[kt]; ct1 = (double)pG[kt]; ct2 = (double)pB[kt];
 cb0 = (double)pR[kb]; cb1 = (double)pG[kb]; cb2 = (double)pB[kb];
 if(alpha > 0.5){
   alpha = alpha-0.5;
   if(beta > 0.5){
     beta = beta-0.5;
     R   = (cr0-cc0)*alpha+(ct0-cc0)*beta+cc0;
     G   = (cr1-cc1)*alpha+(ct1-cc1)*beta+cc1;
     B   = (cr2-cc2)*alpha+(ct2-cc2)*beta+cc2;
   }
   else{
     beta = 0.5-beta;
     R   = (cr0-cc0)*alpha+(cb0-cc0)*beta+cc0;
     G   = (cr1-cc1)*alpha+(cb1-cc1)*beta+cc1;
     B   = (cr2-cc2)*alpha+(cb2-cc2)*beta+cc2;
   }
 }
 else{
   alpha = 0.5-alpha;
   if(beta > 0.5){
     beta = beta-0.5;
     R   = (cl0-cc0)*alpha+(ct0-cc0)*beta+cc0;
     G   = (cl1-cc1)*alpha+(ct1-cc1)*beta+cc1;
     B   = (cl2-cc2)*alpha+(ct2-cc2)*beta+cc2;
   }
   else{
     beta = 0.5-beta;
     R   = (cl0-cc0)*alpha+(cb0-cc0)*beta+cc0;
     G   = (cl1-cc1)*alpha+(cb1-cc1)*beta+cc1;
     B   = (cl2-cc2)*alpha+(cb2-cc2)*beta+cc2;
   }
 }
 *r = R;
 *g = G;
 *b = B;
 return;
}

/* end the Bump and reflection Mapping functions  */

/* The functions for Image mapping and fast       */
/* anti-aliasing line texturing                   */

/* functions for the image mapping   - some old code is still here (#if 0) */

static void Transform_to_Moziac(double alpha, double beta,
                                double *a, double *b){
 long brkindx,brkindy;
 *a=FLOOR(alpha);  *b=FLOOR(beta);
 brkindx=((long)(*a) & 1L);
 brkindy=((long)(*b) & 1L);
 *a=(alpha - *a);  *b=(beta - *b);
 if(brkindx){
   if(brkindy){ *a = 1.0 - *a;  *b = 1.0 - *b; }
   else       { *a = 1.0 - *a; }
 }
 else if(brkindy){ *b = 1.0 - *b; }
 return;
}

/* Only used by Ground and Sky (when ray-tracing) */
long MapFromProjection(imap *Map, short type, short clear,
                        vector p , double *red
                                 , double *green
                                 , double *blue){
 /* return colour for map Map at point p */
 long  ialpha,ibeta,i,j,k,kl,kr,kt,kb,pixel;
 double alpha,beta,a,b,cc0,cc1,cc2,cl0,cl1,cl2,
           cr0,cr1,cr2,ct0,ct1,ct2,cb0,cb1,cb2;
 unsigned char *pP,*pR,*pG,*pB;
 /* checking for valid map must be done elsewhere */
 if(Map->pp < 0)return 0;
 j=Map->xmax;
 if(Map->map == CYLINDER){
   hitcylinder(Map->n,Map->x,Map->y,Map->p
              ,p,&alpha,&beta);
   alpha *= Map->angle;
 }
 else if(Map->map == MAP_SPHERICAL){
   hitsphere(Map->n,Map->x,Map->y,Map->p
              ,p,&alpha,&beta);
 }
 else    /* assume planar map */
   hitpoint(Map->n,Map->x,Map->y,Map->p
           ,p,&alpha,&beta);
 if(type == TILE){
   if(Map->moziac)Transform_to_Moziac(alpha,beta,&a,&b);
   else{
     a=(alpha-floor(alpha));
     b=(beta-floor(beta));
   }
   alpha=a*Map->xmax;
   beta=b*Map->ymax;
 }
 else{
   alpha=Map->xmax*alpha;
   beta =Map->ymax*beta;
 }
 ialpha = (long)alpha; alpha -= (double)ialpha;
 ibeta  = (long)beta;  beta  -= (double)ibeta;
 if(ialpha < 0 || ialpha >= Map->xmax)return 0;
 if(ibeta  < 0 || ibeta  >= Map->ymax)return 0;
 k =ibeta*j+ialpha;
 if(type == TILE){
   if(Map->moziac){
     if(ialpha <   1)kl=1;                 else kl=k-1;
     if(ialpha > j-2)kr=k-1;               else kr=k+1;
     if(ibeta  <   1)kb=k+j;               else kb=k-j;
     if(ibeta  > Map->ymax-2)kt=k-j; else kt=k+j;
   }
   else{
     if(ialpha <   1)kl=k+j-1;                        else kl=k-1;
     if(ialpha > j-2)kr=ibeta*j;                      else kr=k+1;
     if(ibeta  <   1)kb=(Map->ymax-1)*j+ialpha; else kb=k-j;
     if(ibeta  > Map->ymax-2)kt=ialpha;         else kt=k+j;
   }
 }
 else{
   if(ialpha <   1)kl=k; else kl=k-1;
   if(ialpha > j-2)kr=k; else kr=k+1;
   if(ibeta  <   1)kb=k; else kb=k-j;
   if(ibeta  > Map->ymax-2)kt=k; else kt=k+j;
 }
 pR=Map->p24R; pG=Map->p24G; pB=Map->p24B;
 if(clear == TRANSP){
   if(pR[k] == Map->kc[0] && pG[k] == Map->kc[1] &&  pB[k] == Map->kc[2])
     return 0;
 }
 cc0 = (double)pR[k];  cc1 = (double)pG[k];  cc2 = (double)pB[k];
 cl0 = (double)pR[kl]; cl1 = (double)pG[kl]; cl2 = (double)pB[kl];
 cr0 = (double)pR[kr]; cr1 = (double)pG[kr]; cr2 = (double)pB[kr];
 ct0 = (double)pR[kt]; ct1 = (double)pG[kt]; ct2 = (double)pB[kt];
 cb0 = (double)pR[kb]; cb1 = (double)pG[kb]; cb2 = (double)pB[kb];
 if(alpha > 0.5){
   alpha = alpha-0.5;
   if(beta > 0.5){
     beta = beta-0.5;
     *red   = (cr0-cc0)*alpha+(ct0-cc0)*beta+cc0;
     *green = (cr1-cc1)*alpha+(ct1-cc1)*beta+cc1;
     *blue  = (cr2-cc2)*alpha+(ct2-cc2)*beta+cc2;
   }
   else{
     beta = 0.5-beta;
     *red   = (cr0-cc0)*alpha+(cb0-cc0)*beta+cc0;
     *green = (cr1-cc1)*alpha+(cb1-cc1)*beta+cc1;
     *blue  = (cr2-cc2)*alpha+(cb2-cc2)*beta+cc2;
   }
 }
 else{
   alpha = 0.5-alpha;
   if(beta > 0.5){
     beta = beta-0.5;
     *red   = (cl0-cc0)*alpha+(ct0-cc0)*beta+cc0;
     *green = (cl1-cc1)*alpha+(ct1-cc1)*beta+cc1;
     *blue  = (cl2-cc2)*alpha+(ct2-cc2)*beta+cc2;
   }
   else{
     beta = 0.5-beta;
     *red   = (cl0-cc0)*alpha+(cb0-cc0)*beta+cc0;
     *green = (cl1-cc1)*alpha+(cb1-cc1)*beta+cc1;
     *blue  = (cl2-cc2)*alpha+(cb2-cc2)*beta+cc2;
   }
 }
 return 1;
}

long MapByCoordinates(imap *Map,
                        long obj, face *f,
                        short type, short clear,
                        vector p,
                        double *red,
                        double *green,
                        double *blue){
 /* return colour for map pointed to by Map at point p on face f of object # obj */
 long  ialpha,ibeta,i,j,k,kl,kr,kt,kb,pixel;
 double alpha,beta,a,b,cc0,cc1,cc2,cl0,cl1,cl2,
           cr0,cr1,cr2,ct0,ct1,ct2,cb0,cb1,cb2,
           ve1,ve2,ve3,vp1,vp2,vp3,p1,p2,p3,
           detmax,dm,det1,det2,det3,
           t1x,t2x,t1y,t2y;
 unsigned char *pP,*pR,*pG,*pB;
 vertex   *V,*v0,*v1,*v2;
 /* checking for valid map must be done elsewhere */
 if(Map->pp < 0)return 0;
 j=Map->xmax;
 if((V=Object[obj].Vbase) == NULL)return 0;
 v0=(V+(f->V[0]));
 v1=(V+(f->V[1]));
 v2=(V+(f->V[2]));
 ve1=v1->p[0] - v0->p[0]; ve2=v1->p[1] - v0->p[1]; ve3=v1->p[2] - v0->p[2];
 vp1=v2->p[0] - v0->p[0]; vp2=v2->p[1] - v0->p[1]; vp3=v2->p[2] - v0->p[2];
 p1=p[0]-v0->p[0];
 p2=p[1]-v0->p[1];
 p3=p[2]-v0->p[2];
 det1=ve1*vp2-vp1*ve2;
 det2=ve1*vp3-vp1*ve3;
 det3=ve2*vp3-vp2*ve3;
 k=0; detmax=1.e-3;  /* same as TOL in mirrors */
 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 == 0)return 0;
 else 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;
 }
 // Calculate UV map address of pixel
 t1x=f->uu[1] - f->uu[0]; t2x=f->uu[2] - f->uu[0];
 t1y=f->vv[1] - f->vv[0]; t2y=f->vv[2] - f->vv[0];
 alpha = f->uu[0] + a*t1x + b*t2x;
 beta  = f->vv[0] + a*t1y + b*t2y;
 if(type == TILE){
   if(Map->moziac)Transform_to_Moziac(alpha,beta,&a,&b);
   else{
     a=(alpha-floor(alpha));
     b=(beta-floor(beta));
   }
   alpha=a*Map->xmax;
   beta=b*Map->ymax;
 }
 else{
   alpha=Map->xmax*alpha;
   beta =Map->ymax*beta;
 }
 ialpha = (long)alpha; alpha -= (double)ialpha;
 ibeta  = (long)beta;  beta  -= (double)ibeta;
 if(ialpha < 0 || ialpha >= Map->xmax)return 0;
 if(ibeta  < 0 || ibeta  >= Map->ymax)return 0;
 k =ibeta*j+ialpha;

 if(type == TILE){
   if(Map->moziac){
     if(ialpha <   1)kl=1;                 else kl=k-1;
     if(ialpha > j-2)kr=k-1;               else kr=k+1;
     if(ibeta  <   1)kb=k+j;               else kb=k-j;
     if(ibeta  > Map->ymax-2)kt=k-j; else kt=k+j;
   }
   else{
     if(ialpha <   1)kl=k+j-1;                        else kl=k-1;
     if(ialpha > j-2)kr=ibeta*j;                      else kr=k+1;
     if(ibeta  <   1)kb=(Map->ymax-1)*j+ialpha; else kb=k-j;
     if(ibeta  > Map->ymax-2)kt=ialpha;         else kt=k+j;
   }
 }
 else{
   if(ialpha <   1)kl=k; else kl=k-1;
   if(ialpha > j-2)kr=k; else kr=k+1;
   if(ibeta  <   1)kb=k; else kb=k-j;
   if(ibeta  > Map->ymax-2)kt=k; else kt=k+j;
 }

 pR=Map->p24R; pG=Map->p24G; pB=Map->p24B;
 if(clear == TRANSP){
   if(pR[k] == Map->kc[0] && pG[k] == Map->kc[1] &&  pB[k] == Map->kc[2])
     return 0;
 }
 cc0 = (double)pR[k];  cc1 = (double)pG[k];  cc2 = (double)pB[k];
 cl0 = (double)pR[kl]; cl1 = (double)pG[kl]; cl2 = (double)pB[kl];
 cr0 = (double)pR[kr]; cr1 = (double)pG[kr]; cr2 = (double)pB[kr];
 ct0 = (double)pR[kt]; ct1 = (double)pG[kt]; ct2 = (double)pB[kt];
 cb0 = (double)pR[kb]; cb1 = (double)pG[kb]; cb2 = (double)pB[kb];
 if(alpha > 0.5){
   alpha = alpha-0.5;
   if(beta > 0.5){
     beta = beta-0.5;
     *red   = (cr0-cc0)*alpha+(ct0-cc0)*beta+cc0;
     *green = (cr1-cc1)*alpha+(ct1-cc1)*beta+cc1;
     *blue  = (cr2-cc2)*alpha+(ct2-cc2)*beta+cc2;
   }
   else{
     beta = 0.5-beta;
     *red   = (cr0-cc0)*alpha+(cb0-cc0)*beta+cc0;
     *green = (cr1-cc1)*alpha+(cb1-cc1)*beta+cc1;
     *blue  = (cr2-cc2)*alpha+(cb2-cc2)*beta+cc2;
   }
 }
 else{
   alpha = 0.5-alpha;
   if(beta > 0.5){
     beta = beta-0.5;
     *red   = (cl0-cc0)*alpha+(ct0-cc0)*beta+cc0;
     *green = (cl1-cc1)*alpha+(ct1-cc1)*beta+cc1;
     *blue  = (cl2-cc2)*alpha+(ct2-cc2)*beta+cc2;
   }
   else{
     beta = 0.5-beta;
     *red   = (cl0-cc0)*alpha+(cb0-cc0)*beta+cc0;
     *green = (cl1-cc1)*alpha+(cb1-cc1)*beta+cc1;
     *blue  = (cl2-cc2)*alpha+(cb2-cc2)*beta+cc2;
   }
 }
 return 1;
}

long MapTransByCoordinates(imap *Map, short type, long obj, face *f,
                      vector p , double *trans){
 long  ialpha,ibeta,i,j,k,kl,kr,kt,kb,pixel;
 double red,green,blue;
 double alpha,beta,a,b,cc0,cc1,cc2,cl0,cl1,cl2,
           cr0,cr1,cr2,ct0,ct1,ct2,cb0,cb1,cb2,
           ve1,ve2,ve3,vp1,vp2,vp3,p1,p2,p3,
           detmax,dm,det1,det2,det3,
           t1x,t2x,t1y,t2y;
 unsigned char *pP,*pR,*pG,*pB;
 vertex   *V,*v0,*v1,*v2;
 /* checking for valid map must be done elsewhere */
 j=Map->xmaxtran;
 if((V=Object[obj].Vbase) == NULL)return 0;
 v0=(V+(f->V[0]));
 v1=(V+(f->V[1]));
 v2=(V+(f->V[2]));
 ve1=v1->p[0] - v0->p[0]; ve2=v1->p[1] - v0->p[1]; ve3=v1->p[2] - v0->p[2];
 vp1=v2->p[0] - v0->p[0]; vp2=v2->p[1] - v0->p[1]; vp3=v2->p[2] - v0->p[2];
 p1=p[0]-v0->p[0];
 p2=p[1]-v0->p[1];
 p3=p[2]-v0->p[2];
 det1=ve1*vp2-vp1*ve2;
 det2=ve1*vp3-vp1*ve3;
 det3=ve2*vp3-vp2*ve3;
 k=0; detmax=1.e-3;  /* same as TOL in mirrors */
 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 == 0)return 0;
 else 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;
 }
 // Calculate UV map address of pixel
 t1x=f->uu[1] - f->uu[0]; t2x=f->uu[2] - f->uu[0];
 t1y=f->vv[1] - f->vv[0]; t2y=f->vv[2] - f->vv[0];
 alpha = f->uu[0] + a*t1x + b*t2x;
 beta  = f->vv[0] + a*t1y + b*t2y;

 if(type == TILE){
   if(Map->moziac)Transform_to_Moziac(alpha,beta,&a,&b);
   else{
     a=(alpha-floor(alpha));
     b=(beta-floor(beta));
   }
   alpha=a*Map->xmaxtran;
   beta=b*Map->ymaxtran;
 }
 else{
   alpha=Map->xmaxtran*alpha;
   beta =Map->ymaxtran*beta;
 }
 ialpha = (long)alpha; alpha -= (double)ialpha;
 ibeta  = (long)beta;  beta  -= (double)ibeta;

 if(ialpha < 0 || ialpha >= Map->xmaxtran)return 0;
 if(ibeta  < 0 || ibeta  >= Map->ymaxtran)return 0;
 k =ibeta*j+ialpha;
 if(type == TILE){
   if(Map->moziac){
     if(ialpha <   1)kl=1;                 else kl=k-1;
     if(ialpha > j-2)kr=k-1;               else kr=k+1;
     if(ibeta  <   1)kb=k+j;               else kb=k-j;
     if(ibeta  > Map->ymaxtran-2)kt=k-j; else kt=k+j;
   }
   else{
     if(ialpha <   1)kl=k+j-1;                        else kl=k-1;
     if(ialpha > j-2)kr=ibeta*j;                      else kr=k+1;
     if(ibeta  <   1)kb=(Map->ymaxtran-1)*j+ialpha; else kb=k-j;
     if(ibeta  > Map->ymaxtran-2)kt=ialpha;         else kt=k+j;
   }
 }
 else{
   if(ialpha <   1)kl=k; else kl=k-1;
   if(ialpha > j-2)kr=k; else kr=k+1;
   if(ibeta  <   1)kb=k; else kb=k-j;
   if(ibeta  > Map->ymaxtran-2)kt=k; else kt=k+j;
 }
 pR=Map->p24Rtran; pG=Map->p24Gtran; pB=Map->p24Btran;
 cc0 = (double)pR[k];  cc1 = (double)pG[k];  cc2 = (double)pB[k];
 cl0 = (double)pR[kl]; cl1 = (double)pG[kl]; cl2 = (double)pB[kl];
 cr0 = (double)pR[kr]; cr1 = (double)pG[kr]; cr2 = (double)pB[kr];
 ct0 = (double)pR[kt]; ct1 = (double)pG[kt]; ct2 = (double)pB[kt];
 cb0 = (double)pR[kb]; cb1 = (double)pG[kb]; cb2 = (double)pB[kb];
 if(alpha > 0.5){
   alpha = alpha-0.5;
   if(beta > 0.5){
     beta = beta-0.5;
     red   = (cr0-cc0)*alpha+(ct0-cc0)*beta+cc0;
     green = (cr1-cc1)*alpha+(ct1-cc1)*beta+cc1;
     blue  = (cr2-cc2)*alpha+(ct2-cc2)*beta+cc2;
   }
   else{
     beta = 0.5-beta;
     red   = (cr0-cc0)*alpha+(cb0-cc0)*beta+cc0;
     green = (cr1-cc1)*alpha+(cb1-cc1)*beta+cc1;
     blue  = (cr2-cc2)*alpha+(cb2-cc2)*beta+cc2;
   }
 }
 else{
   alpha = 0.5-alpha;
   if(beta > 0.5){
     beta = beta-0.5;
     red   = (cl0-cc0)*alpha+(ct0-cc0)*beta+cc0;
     green = (cl1-cc1)*alpha+(ct1-cc1)*beta+cc1;
     blue  = (cl2-cc2)*alpha+(ct2-cc2)*beta+cc2;
   }
   else{
     beta = 0.5-beta;
     red   = (cl0-cc0)*alpha+(cb0-cc0)*beta+cc0;
     green = (cl1-cc1)*alpha+(cb1-cc1)*beta+cc1;
     blue  = (cl2-cc2)*alpha+(cb2-cc2)*beta+cc2;
   }
 }
 *trans = (1.0- 1.3071895e-3*(red+green+blue)); /* scaled 0-1 Back is fully opaque white is transp */
 return 1;
}

long MapTextureByCoordinates(long obj,face *f, vector p,
     double *alpha, double *beta, double *gamma){
 /* return 2 D coords on face f of object # obj */
 long k;
 double a,b,
        ve1,ve2,ve3,vp1,vp2,vp3,p1,p2,p3,
        detmax,dm,det1,det2,det3,
        t1x,t2x,t1y,t2y,t1w,t2w;
 unsigned char *pP,*pR,*pG,*pB;
 vertex   *V,*v0,*v1,*v2;
 /* checking for valid map must be done elsewhere */
 if((V=Object[obj].Vbase) == NULL)return 0;
 v0=(V+(f->V[0])); 
 v1=(V+(f->V[1])); 
 v2=(V+(f->V[2])); 
 ve1=v1->p[0] - v0->p[0]; ve2=v1->p[1] - v0->p[1]; ve3=v1->p[2] - v0->p[2];
 vp1=v2->p[0] - v0->p[0]; vp2=v2->p[1] - v0->p[1]; vp3=v2->p[2] - v0->p[2];
 p1=p[0]-v0->p[0];
 p2=p[1]-v0->p[1];
 p3=p[2]-v0->p[2];
 det1=ve1*vp2-vp1*ve2;
 det2=ve1*vp3-vp1*ve3;
 det3=ve2*vp3-vp2*ve3;
 k=0; detmax=1.e-3;  /* same as TOL in mirrors */
 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 == 0)return 0;
 else 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;
 }
 // Calculate UV map address of pixel
 t1x=f->uu[1] - f->uu[0]; t2x=f->uu[2] - f->uu[0];
 t1y=f->vv[1] - f->vv[0]; t2y=f->vv[2] - f->vv[0];
 t1w=f->ww[1] - f->ww[0]; t2w=f->ww[2] - f->ww[0];
 *alpha = f->uu[0] + a*t1x + b*t2x;
 *beta  = f->vv[0] + a*t1y + b*t2y;
 *gamma = f->ww[0] + a*t1w + b*t2w;
// *gamma = 0.5;
 return 1;
}


long ShadedScanline(
                 long scanline ,double d2,
                 UCHAR *Red, UCHAR *Green, UCHAR *Blue,
                 double *Zground, long Wsky, double Rsky,
                 long *Zglass[], long *ZglassO[], double *Zglassd[],
                 long *Zbuffer, long *Zobject,
                 double *Zdepth, double *Zposn){
 face    *f;
 vertex  *v;
 imap    *Map;
 matl    *Mat;
 BOOL    bSmooth;
 long i,j,k,f1,obj;
 double det,alpha,beta;
 vector n,nn,dn,v0,v1,v2;
 UCHAR *ZobjChar;
 if(VGA_screen == NULL)ZobjChar=DummyShadow;
 else ZobjChar = &VGA_screen[scanline*XMAX]; /* shadow for this line */
 i=XMIN;
LOOP:
 f1=Zbuffer[i];
 if(f1 >= 0){
   obj=Zobject[i];
   v=Object[obj].Vbase;
   f=Object[obj].Fbase + f1;
   if(f->ffmap >= 0){                 /* mapped get map id */
     if((long)f->ffmap >= Object[obj].NoMaps)Map=NULL;
     else Map=(Object[obj].Mbase+(long)f->ffmap);
   }
   else Map=NULL;
   // Get the material ID or NULL
   if(f->ffmat < 0 || f->ffmat >= Object[obj].NoMats)Mat=NULL;
   else Mat=(Object[obj].Tbase+(long)f->ffmat);
 /* calculate the  alpha and beta for starting value */
   if(Mat != NULL)bSmooth=Mat->bSmooth;
   else           bSmooth=f->bSmooth;
   if(!bSmooth){   /* no smoothing required n is face normal */
     VECCOPY(f->n,n)
   }
   else{ /* interpolate normal over the face */
     det=(
         (v[f->V[2]].y - v[f->V[0]].y) * (v[f->V[1]].x - v[f->V[0]].x)
        -(v[f->V[2]].x - v[f->V[0]].x) * (v[f->V[1]].y - v[f->V[0]].y)
          );
     if(det == 0){
       VECCOPY(f->n,n)
     }
     else{
       VECCOPY(f->pn[0],v0)
       VECCOPY(f->pn[1],v1)
       VECCOPY(f->pn[2],v2)
       det=1.0/det;
       alpha= (double)(
        (v[f->V[2]].y - v[f->V[0]].y) * (Zposn[i] - v[f->V[0]].x)
       -(v[f->V[2]].x - v[f->V[0]].x) * ((double)scanline - v[f->V[0]].y)
                      )*det;
       beta= (double)(
        ((double)scanline - v[f->V[0]].y) * (v[f->V[1]].x - v[f->V[0]].x)
       -(Zposn[i] - v[f->V[0]].x) * (v[f->V[1]].y - v[f->V[0]].y)
                    )*det;
       n[0]=v0[0] + alpha*(v1[0] - v0[0])
                  +  beta*(v2[0] - v0[0]);
       n[1]=v0[1] + alpha*(v1[1] - v0[1])
                  +  beta*(v2[1] - v0[1]);
       n[2]=v0[2] + alpha*(v1[2] - v0[2])
                  +  beta*(v2[2] - v0[2]);
       normalize(n);
     }
   }
   GetRowPixelValue(scanline,i,Map,Mat,f,obj,n,Red,Green,Blue,
                      Zdepth,ZobjChar); /* pixel 1*/
   if(Zglass[0][i] > -1)AttenuateGlass(scanline,i,Red,Green,Blue,
                                       Zglass,ZglassO,Zglassd,Zposn);
   j=i+1;  /* interpolate normal over pixels showing same face */
   if(j < XMAX){
    while((Zbuffer[j] == f1) && (Zobject[j] == Zobject[i]))
                                {j++; if(j == XMAX)break; }   /* how many */
    if(j > i+1){ /* more than one pixel */
      if(j > i+2){ /* more than 2 pixels */
        if((det==0) || !(bSmooth)){
          dn[0]=0; dn[1]=0; dn[2]=0;
        }
        else{
          dn[0]=( (v[f->V[1]].y - v[f->V[2]].y) * v0[0]
                + (v[f->V[2]].y - v[f->V[0]].y) * v1[0]
                + (v[f->V[0]].y - v[f->V[1]].y) * v2[0]
                )*det;
          dn[1]=( (v[f->V[1]].y - v[f->V[2]].y) * v0[1]
                + (v[f->V[2]].y - v[f->V[0]].y) * v1[1]
                + (v[f->V[0]].y - v[f->V[1]].y) * v2[1]
                )*det;
          dn[2]=( (v[f->V[1]].y - v[f->V[2]].y) * v0[2]
                + (v[f->V[2]].y - v[f->V[0]].y) * v1[2]
                + (v[f->V[0]].y - v[f->V[1]].y) * v2[2]
                )*det;
        }
        for(k=i+1; k<j-1; k++){
          VECSUM(n,dn,n)
          VECCOPY(n,nn)
          normalize(nn);
          GetRowPixelValue(scanline,k,Map,Mat,f,obj,nn,Red,Green,Blue,
                             Zdepth,ZobjChar);
          if(Zglass[0][k] > -1)AttenuateGlass(scanline,k,Red,Green,Blue,
                                              Zglass,ZglassO,Zglassd,Zposn);
        }
      }
      /* do last pixel this must be separate - may not be on integer boundary */
      if((det==0) || !(bSmooth))VECCOPY(f->n,n)
      else{
        alpha = (double)(
         (v[f->V[2]].y - v[f->V[0]].y) * (Zposn[j-1] - v[f->V[0]].x)
        -(v[f->V[2]].x - v[f->V[0]].x) * ((double)scanline - v[f->V[0]].y)
                        )*det;
        beta = (double)(
         ((double)scanline - v[f->V[0]].y) * (v[f->V[1]].x - v[f->V[0]].x)
        -(Zposn[j-1] - v[f->V[0]].x) * (v[f->V[1]].y - v[f->V[0]].y)
                       )*det;
        n[0]=v0[0] + alpha*(v1[0] - v0[0])
                   +  beta*(v2[0] - v0[0]);
        n[1]=v0[1] + alpha*(v1[1] - v0[1])
                   +  beta*(v2[1] - v0[1]);
        n[2]=v0[2] + alpha*(v1[2] - v0[2])
                   +  beta*(v2[2] - v0[2]);
        normalize(n);
      }
      GetRowPixelValue(scanline,j-1,Map,Mat,f,obj,n,Red,Green,Blue,
                         Zdepth,ZobjChar);
      if(Zglass[0][j-1] > -1)AttenuateGlass(scanline,j-1,Red,Green,Blue,
                                            Zglass,ZglassO,Zglassd,Zposn);
    }
   }
   i=j;
 }
 else{  /* if no object get back plane values */
  PutBackPlanes(i,scanline,d2,Red,Green,Blue,
                Zground,Zobject,Zdepth,
                ZobjChar,Wsky,Rsky);
  if(Zglass[0][i] > -1)AttenuateGlass(scanline,i,Red,Green,Blue,
                                      Zglass,ZglassO,Zglassd,Zposn);
  i++;
 }
 if(i < XMAX)goto LOOP;
 return 0;
}

long FastAntiAliasShadedScanline(
                      long scanline ,long aaline, double d2,
                      UCHAR *Red, UCHAR *Green, UCHAR *Blue,
                      double *Zground, long Wsky, double Rsky,
                      long *Xbuffer, long *Xobject,
                      long *Zglass[], long *ZglassO[], double *Zglassd[],
                      long *Zbuffer, long *Zobject,
                      double *Zdepth, double *Zposn){
 face   *f;             /* super sample if faces are different  */
 vertex *v;
 imap   *Map;
 matl   *Mat;
 BOOL   bSmooth;
 long i,j,k,l,f1,obj,AlreadyCalculated;
 double det,alpha,beta;
 vector n,v0,v1,v2;
 UCHAR  *ZobjChar;
 if(VGA_screen == NULL)ZobjChar=DummyShadow;
 else ZobjChar = &VGA_screen[scanline*XMAX];
 i=0; l=i-1;
 for(k=0;k<XMAX;k+=aaXstep){
   for(j=0;j<aaXstep;j++){
     f1=Zbuffer[i]; obj=Zobject[i];
     if(f1 == Xbuffer[i] && obj == Xobject[i])AlreadyCalculated=1;
     else if(j == 0)AlreadyCalculated=0;
     else if(f1 == Zbuffer[i-1] && obj == Zobject[i-1]){
       AlreadyCalculated=1;
       Red[i] = Red[l]; Green[i]=Green[l]; Blue[i]=Blue[l];
     }
     else AlreadyCalculated=0;
     if(!AlreadyCalculated){
       if(f1 >= 0){
         v=Object[obj].Vbase;
         f=Object[obj].Fbase + f1;
         if(f->ffmap >= 0){ /* mapped */
           if((long)f->ffmap >= Object[obj].NoMaps)Map=NULL;
           else Map=(Object[obj].Mbase+(long)f->ffmap);
         }
         else Map=NULL;
         if(f->ffmat < 0 || f->ffmat >= Object[obj].NoMats)Mat=NULL;
         else Mat=(Object[obj].Tbase+(long)f->ffmat);
         if(Mat != NULL)bSmooth=Mat->bSmooth;
         else           bSmooth=f->bSmooth;
         if(!bSmooth)VECCOPY(f->n,n)
         else{
           det=(
             (v[f->V[2]].y - v[f->V[0]].y) * (v[f->V[1]].x - v[f->V[0]].x)
              -(v[f->V[2]].x - v[f->V[0]].x) * (v[f->V[1]].y - v[f->V[0]].y)
               );
           if(det == 0){ VECCOPY(f->n,n) }
           else{
             VECCOPY(f->pn[0],v0)
             VECCOPY(f->pn[1],v1)
             VECCOPY(f->pn[2],v2)
             det=1.0/det;
             alpha= (double)(
                     (v[f->V[2]].y - v[f->V[0]].y) *
                         (Zposn[i] - v[f->V[0]].x)
                    -(v[f->V[2]].x - v[f->V[0]].x) *
                 ((double)scanline - v[f->V[0]].y)
                            )*det;
             beta= (double)(
                  ((double)scanline - v[f->V[0]].y) *
                      (v[f->V[1]].x - v[f->V[0]].x)
                         -(Zposn[i] - v[f->V[0]].x) *
                      (v[f->V[1]].y - v[f->V[0]].y)
                           )*det;
             n[0]=v0[0] + alpha*(v1[0] - v0[0]) +  beta*(v2[0] - v0[0]);
             n[1]=v0[1] + alpha*(v1[1] - v0[1]) +  beta*(v2[1] - v0[1]);
             n[2]=v0[2] + alpha*(v1[2] - v0[2]) +  beta*(v2[2] - v0[2]);
             normalize(n);
           }
         }
         GetRowPixelValue(scanline,i,Map,Mat,f,obj,n,Red,Green,Blue,
                            Zdepth,ZobjChar);
         if(Zglass[0][i] > -1)AttenuateGlass(scanline,i,Red,Green,Blue,
                                             Zglass,ZglassO,Zglassd,Zposn);
       }
       else{
         PutBackPlanes(i,scanline,d2,Red,Green,Blue,
                       Zground,Zobject,Zdepth,ZobjChar,Wsky,Rsky);
         if(Zglass[0][i] > -1)AttenuateGlass(scanline,i,Red,Green,Blue,
                                             Zglass,ZglassO,Zglassd,Zposn);
       }
     }
     i++; l++;
   }
 }
 return 0;
}

/* end of Image mapping and fast AA line texturing     */

/* code for rendering specific elements such as ground */

void GetSkyLocation(void){
 double b;
 long i;
 for(i=XMIN;i<XMAX;i++){
   if(Sky.type == SKYMAPPED){
     b= (CamPhi+atan2(D1[i],1.000))/PI*180.0 * 4.0 ; /* 4 maps full rotation */
     while(b <  -180.0)b = b + 360.0;
     while(b >=  180.0)b = b - 360.0;
     ZskyBuffer[XMAX-1-i]=(short)((b+180.0)*skyXscale);
   }
   else{
     ZskyBuffer[XMAX-1-i]=0;
   }
 }
}

long PutBackIntersection(long scanline, double d2,
                         UCHAR *Alfa, double *Zground,
                         long *Zglass[], long *ZglassO[], double *Zglassd[],
                         long *Zbuffer, long *Zobject,
                         double *Zdepth, double *Zposn, double *Rsky){
 long i,j,Wsky;  /* background  setup */
 double b,c,mu;
 if(Sky.type == SKYGRADED){
   //*Rsky=fabs((CamAlpha+atan2(d2,1.000))/PI*2.0);  /* 0 Hor  1 +/- Zenith */
   *Rsky=fabs((double)(Ycentre-scanline))/(double)Ycentre;     
 }
 else if(Sky.type == SKYMAPPED){ /* CamAlpha is 0 on horizontal +90 look up*/
  Wsky=0; /* Wsky scanline V pos. in env.bru 0 is brush top*/
  b=90.0 - (CamAlpha+atan2(d2,1.000))/PI*180.0 * 2;
  while(b <   0.0)b += 90.0;
  while(b > 180.0)b -= 90.0;
  if(b > 0){          /* b is angle in range 0-180 for view to scanline */
    if(b > 90)Wsky=(short)((180-b)*skyYscale);      /* 0 is up 180 down */
    else      Wsky=(short)(b*skyYscale);            /* b > 90 picture is */
    Wsky = min(Wsky,skyBmax);                       /* inverted */
  }
 }
 for(i=XMIN;i<XMAX;i++){
   Zbuffer[i] = -1;
   Zdepth [i] = FARAWAY;
   Zposn  [i] = (double)i+0.5;
   for(j=0;j<NGLASS;j++){
     Zglassd[j][i] = FARAWAY;
     ZglassO[j][i] = -1;
     Zglass [j][i] = -1;  /* this will be non zero if a glass face is on top */
   }
   Zobject[i] = -2;
   Alfa   [i] =  0;
   if(R_Nground > 0){
     b=Ground.n[0]*D1[i]+Ground.n[1]+Ground.n[2]*d2; // Z buffer
//{
// vector dr;     //  if Z buffer is to contain actual depths
// dr[0]=D1[i]; dr[1]=1.0; dr[2]=d2;
// b=DOT(Ground.n,dr);
//}
     if(fabs(b) > 1.e-4){
       c=dot(Ground.n,Ground.p);
       mu=c/b;
       if(mu > 0.0){
         Zdepth[i] = mu;
         Zground[i]= mu;
         Zobject[i]= -1;
       }
     }
   }
 }
 return Wsky;
}

void GetGroundLight(vector p0, double *R, double *G ,double *B){
/* the R G B values may exceed 255.0 they must be clipped elsewhere */
 long j,brkindx;
 double c,cc,alpha,beta,Iab,IsR,IsG,IsB,Bump_Scale,IdR,IdG,IdB,xRGB[3]={0.0,0.0,0.0};
 vector p1,p2,ripple,gn,gnn,bv;
 double red,green,blue,map;
 IdR=IdG=IdB=0.0;
    if(Ground.texture > 0 && (Ground.texture < 4 || Ground.texture > 9)){
      switch(Ground.texture){
      case 1:  /* single wave  */
         calcripple(p0,Ground.p,wave_phase,Ground.lengthy,1.0,
         Ground.height,ripple); /* typical height 2.0 */
         Bump_Scale = - dot(p0,Ground.n)/R_length(p0);
         vecscale(Bump_Scale,ripple,ripple);
         vecsum(ripple,Ground.n,gn);
         normalize(gn);
         IsR=IsG=IsB=0.0;
         break;
      case 2: /* ripples using wave 2 */
// NEED FIXED because of hitpoint
//        waves2(p0,Ground.n,Ground.dx,Ground.dyy,Ground.p,gn,
//        Rip_Freq,wave_phase,Ground.height);  /* 0.2 is good for height */
        IsR=IsG=IsB=0.0;
        break;
      case 3: /* ocean multiple waves texture */
        Ocean(p0,Ground.n,gn);
        IsR=IsG=IsB=0.0;
        break;
      case 10: /* rough surface  ie bumpy */
        hitpoint(Ground.n,Ground.dx,Ground.dy,Ground.p,p0,&alpha,&beta);
        Dnoise2(gn,alpha*10.0,beta*10.0,0.0);                  /* scale bumps into */
        Bump_Scale = dot(p0,Ground.n)/R_length(p0);     /* distance*/
        Bump_Scale = min(0.5, Bump_Scale*2.0)*Ground.height*0.5;
        vecscale(Bump_Scale,gn,gn);
        vecsum(Ground.n,gn,gn);
        normalize(gn);
        break;
      case 16: /* external textures  this is not used */ {
        }
        break;
      default:
        veccopy(Ground.n,gn);
        break;
      }
    }
    else {
      veccopy(Ground.n,gn);
    }

    for(j=0;j<Nlights;j++)if(Lights[j].type != SHADOWS &&
                             Lights[j].type != DUMMYL){
      vecsub(Lights[j].p,p0,p1);
      cc=R_length(p1);
      vecscale(1.0/cc,p1,p2);
      Iab = Ground.ambient_light+(1.0-Ground.ambient_light)*fabs(DOT(p2,gn));
      if(Lights[j].dc){
        cc = DepthCueLight(j,cc);
        Iab *= cc;
      }
      if(Lights[j].type == SPOTLIGHT){
        if((alpha=DOT(p2,Lights[j].d)) < Lights[j].dot2)
          Iab=Ground.ambient_light;
        else if(alpha < Lights[j].dot1){
          Iab=Ground.ambient_light+(Iab-Ground.ambient_light)*
               (alpha-Lights[j].dot2)*Lights[j].falloff;
        }
        if(Iab > Ground.ambient_light && Lights[j].shadow_buffer != NULL)
            PointInShadow(p0,j,&Iab); /* for spotlights */
      }
      IdR += Iab*Lights[j].Ire;
      IdG += Iab*Lights[j].Igr;
      IdB += Iab*Lights[j].Ibl;
      if(Ground.shiny){
        vecsub(p1,p0,bv);
        normalize(bv);
        c=fabs(dot(bv,gn));
        if(c > 0.0){
          Iab = spec_gro[min(max_spec,(long)(c*max_spec))];
          if(Lights[j].dc)Iab *= cc;
          IsR += Iab*Lights[j].Ire;
          IsG += Iab*Lights[j].Igr;
          IsB += Iab*Lights[j].Ibl;
        }
      }
    }
    switch(Ground.texture){
    case 1:  /* waves have specular highlights */
    case 2:
    case 3:
      *R = min( ((double)Ground.color[0]*IdR + IsR*255), 255);
      *G = min( ((double)Ground.color[1]*IdG + IsG*255), 255);
      *B = min( ((double)Ground.color[2]*IdB + IsB*255), 255);
      break;
    case 4:  /* checkered ground */
      hitpoint2(Ground.n,Ground.dx,Ground.dy,Ground.p,p0,&alpha,&beta);
      brkindx = (long)FLOOR(alpha) + (long)FLOOR(beta);
      if(brkindx & 1){
        *R = (double)Ground.acolor[0]*IdR;
        *G = (double)Ground.acolor[1]*IdG;
        *B = (double)Ground.acolor[2]*IdB;
      }
      else{
        *R = (double)Ground.color[0]*IdR;
        *G = (double)Ground.color[1]*IdG;
        *B = (double)Ground.color[2]*IdB;
      }
      break;
    case 5: /* grid texture */
      hitpoint(Ground.n,Ground.dx,Ground.dy,Ground.p,p0,&alpha,&beta);
      alpha=(alpha-floor(alpha));
      beta=(beta-floor(beta));
      if(beta < 0.1 || alpha < 0.1){
        *R = (double)Ground.acolor[0]*IdR;
        *G = (double)Ground.acolor[1]*IdG;
        *B = (double)Ground.acolor[2]*IdB;
      }
      else{
        *R = (double)Ground.color[0]*IdR;
        *G = (double)Ground.color[1]*IdG;
        *B = (double)Ground.color[2]*IdB;
      }
      break;
    case 6: /* stripes texure */
      hitpoint(Ground.n,Ground.dx,Ground.dy,Ground.p,p0,&alpha,&beta);
      alpha=(alpha-floor(alpha));
      beta=(beta-floor(beta));
      if(alpha < 0.1){
        *R = (double)Ground.acolor[0]*IdR;
        *G = (double)Ground.acolor[1]*IdG;
        *B = (double)Ground.acolor[2]*IdB;
      }
      else{
        *R = (double)Ground.color[0]*IdR;
        *G = (double)Ground.color[1]*IdG;
        *B = (double)Ground.color[2]*IdB;
      }
      break;
    case 7: /* dot texture */
      hitpoint(Ground.n,Ground.dx,Ground.dyy,Ground.p,p0,&alpha,&beta);
      alpha=(alpha-floor(alpha))-0.5;
      beta=(beta-floor(beta))-0.5;
      if(beta*beta+alpha*alpha < Ground.range){
        *R = (double)Ground.acolor[0]*IdR;
        *G = (double)Ground.acolor[1]*IdG;
        *B = (double)Ground.acolor[2]*IdB;
      }
      else{
        *R = (double)Ground.color[0]*IdR;
        *G = (double)Ground.color[1]*IdG;
        *B = (double)Ground.color[2]*IdB;
      }
      break;
    case 8: /* cross texture */
      hitpoint(Ground.n,Ground.dx,Ground.dyy,Ground.p,p0,&alpha,&beta);
      alpha=fabs((alpha-floor(alpha))-0.5);
      beta=fabs((beta-floor(beta))-0.5);
      if((beta < 0.3 && alpha < 0.1) || (alpha < 0.3 && beta < 0.1)){
        *R = (double)Ground.acolor[0]*IdR;
        *G = (double)Ground.acolor[1]*IdG;
        *B = (double)Ground.acolor[2]*IdB;
      }
      else{
        *R = (double)Ground.color[0]*IdR;
        *G = (double)Ground.color[1]*IdG;
        *B = (double)Ground.color[2]*IdB;
      }
      break;
    case 9: /* brick texture */
      hitpoint(Ground.n,Ground.dx,Ground.dy,Ground.p,p0,&alpha,&beta);
      alpha=(alpha-floor(alpha));
      beta=(beta-floor(beta));
      if((beta > 0.95 || (beta > 0.45 && beta < 0.5)) ||
         (beta < 0.5 && ( (alpha > 0.45 && alpha < 0.5) || alpha > 0.95) ) ||
         (beta > 0.5 && ((alpha > 0.25 && alpha < 0.3) ||
         (alpha > 0.7 && alpha < 0.75)))
         ){
        *R = (double)Ground.acolor[0]*IdR;
        *G = (double)Ground.acolor[1]*IdG;
        *B = (double)Ground.acolor[2]*IdB;
      }
      else{
        *R = (double)Ground.color[0]*IdR;
        *G = (double)Ground.color[1]*IdG;
        *B = (double)Ground.color[2]*IdB;
      }
      break;
    case 10:
      *R = (double)Ground.color[0]*IdR;
      *G = (double)Ground.color[1]*IdG;
      *B = (double)Ground.color[2]*IdB;
      break;  /* bumpy */
    case 16:  /* external texture */
      *R = xRGB[0]*IdR; *G = xRGB[1]*IdG; *B = xRGB[2]*IdB;
      break;
    default:  /* ground is mapped by a brush               */
      map=0;  /* ground brush has been successfully loaded */

      if(Ground.mapped){
         map=MapFromProjection(&Ground.map,Ground.maptype,0,
                      p0,&red,&green,&blue);
      }
      if(map == 1){
        *R = (red*IdR);
        *G = (green*IdG);
        *B = (blue*IdB);
      }
      else{
        *R = (double)Ground.color[0]*IdR;
        *G = (double)Ground.color[1]*IdG;
        *B = (double)Ground.color[2]*IdB;
      }
    break;
    }
}

/* end of code for ground and other elements           */

double DepthCueLight(long j, double cc){
 double d;
 return 1.0/(Lights[j].dc_c + cc*(Lights[j].dc_l + cc*Lights[j].dc_q));
}

void GetPixelLighting(matl *Mat, vector p, vector n, double pdotn, double ddotn,
                      long obj, face *f, 
                      double *IdR, double *IdG, double *IdB,
                      double *IsR, double *IsG, double *IsB){
 long j,matsp;
 double cc,alpha,c,ldotn,falloff,Iab,amblight; 
 vector bv,pl,pr;
 BOOL bShiny;
 for(j=0;j<Nlights;j++)if(Lights[j].type != SHADOWS &&
                          Lights[j].type != DUMMYL){
   VECSUB(Lights[j].p,p,pl)
   cc=R_length(pl);
   VECSCALE(1.0/cc,pl,pr)
/*   ldotn=DOT(pr,f->n); this was used to try to eliminate fringe pixels
     but it appears to cause some faces to de-phong incorrectly */
   ldotn=DOT(pr,n);
   if((ldotn > 1.e-6)){
     amblight=Object[obj].ambient_light;
     //c = max(0.0,DOT(pr,n));
     c = ldotn;
     if(Lights[j].type == SPOTLIGHT){
/* note the < is because DOT is max when pr and n are aligned */
/* dot2 is the outer cone, dot1 is the inner cone             */
/* d of light is in OPPOSITE dirn. to the dirn of the light   */
/* pr is vector from point to light thus if pr.n = 1 light is */
/* pointing straight at the point                             */
       if((alpha=DOT(pr,Lights[j].d)) < Lights[j].dot2)falloff=0;
       else if(alpha < Lights[j].dot1){ /* in between cones */
         falloff=(alpha-Lights[j].dot2)*Lights[j].falloff;
       }
       else falloff=1.0; /* otherwise inside inner cone and no atten */
       c *= falloff;
       if(c > 0.0 && Lights[j].shadow_buffer != NULL &&
         Object[obj].show_shadow)PointInShadow(p,j,&c);
     }
     Iab = amblight+(1.0-amblight)*c;
   }
   else{  // Light backside
     Iab=Object[obj].ambient_light;
     falloff=1.0;
   }
   if(Lights[j].dc){
     cc = DepthCueLight(j,cc);
     Iab *= cc;
   }
   *IdR += Iab*Lights[j].Ire;
   *IdG += Iab*Lights[j].Igr;
   *IdB += Iab*Lights[j].Ibl;

   if(Mat != NULL){bShiny=Mat->bShiny; matsp=Mat->gloss;}
   else           {bShiny=f->bShiny;;  matsp=8;}
   if(bShiny){/* shiny  switch is on  */
     VECSUB(pl,p,bv)                /* calculate specular highlight */
     normalize(bv);
     c=DOT(bv,n);
     if(c > 0.0)Iab = spec_pow[matsp][min(max_spec,(long)(c*max_spec))];
     else Iab = 0;
     if(Lights[j].type == SPOTLIGHT)Iab *= falloff;
     if(Lights[j].dc)Iab *= cc;
     *IsR += Iab*Lights[j].Ire;
     *IsG += Iab*Lights[j].Igr;
     *IsB += Iab*Lights[j].Ibl;
   }
 }
}

static void PutBackPlanes(long i, long scanline, double d2,
                   UCHAR *Red, UCHAR *Green, UCHAR *Blue,
                   double *Zground, long *Zobject, double *Zdepth,
                   UCHAR *ZobjChar, long Wsky, double Rsky){
 vector p0;
 double red,green,blue;
 long ix;
   if(Zobject[i] == -1){   /* ground at this pixel */
    p0[0]=(((double)i-(double)Xcentre)/scalex*Zground[i]);
    p0[1]=min(FARAWAY,Zdepth[i]);
    p0[2]=(((double)Ycentre-(double)scanline)/scaley*Zground[i]);
    GetGroundLight(p0,&red,&green,&blue);
    Red  [i] = (unsigned char)min(255.0,red);
    Green[i] = (unsigned char)min(255.0,green);
    Blue [i] = (unsigned char)min(255.0,blue);
    if((ZobjChar[i] & 1) == 1){
      Red[i]   *= shadow_density;     /* do the shadows before reflections*/
      Green[i] *= shadow_density;
      Blue[i]  *= shadow_density;
    }
    if(mirror_ground == 1){
      Add_Ground_Mirror(scanline,i,p0,&Red[i],&Green[i],&Blue[i]);
    }
   }
   else if(Nsky > 0){
     /* if the map could not be loaded use this as well */
     if(Sky.type == SKYGRADED){
       Red  [i]=(unsigned char)(
                   ((double)Sky.Zenith[0]-(double)Sky.Horizon[0])*Rsky
                  + (double)Sky.Horizon[0]);
       Green[i]=(unsigned char)(
                   ((double)Sky.Zenith[1]-(double)Sky.Horizon[1])*Rsky
                  + (double)Sky.Horizon[1]);
       Blue [i]=(unsigned char)(
                   ((double)Sky.Zenith[2]-(double)Sky.Horizon[2])*Rsky
                  + (double)Sky.Horizon[2]);
     }
     else if(Sky.type <= SKYCOLOUR || !Sky.mapped){
       Red  [i]=Sky.Zenith[0];
       Green[i]=Sky.Zenith[1];
       Blue [i]=Sky.Zenith[2];
     }
     else if(Sky.type == SKYMAPPED){  /* Map the sky - do it directly */
       ix=Wsky*Sky.map.xmax+ZskyBuffer[i];
       Red[i]   = Sky.map.p24R[ix];
       Green[i] = Sky.map.p24G[ix];
       Blue[i]  = Sky.map.p24B[ix];
     }
     else if(Sky.type == BACKDROP){  /* Put in back drop */
       ix=((scanline*Sky.map.ymax)/YMAX)*Sky.map.xmax
         +(i*Sky.map.xmax)/XMAX;
       Red[i]   = Sky.map.p24R[ix];
       Green[i] = Sky.map.p24G[ix];
       Blue[i]  = Sky.map.p24B[ix];
     }
     else if(Sky.type == SKYCUBE){ // use environment map 
       MapEnvironment(i,scanline,&Red[i],&Green[i],&Blue[i], &(Sky.map));
     }
   }
   else{ Red[i]=0;   Green[i]=0;   Blue[i]=0; }
}

static void GetBackdropValue(long row, long col, vector p, vector d,
                             double *colour){
 /* this is used by the Tracer to get the background */
 vector dp,gp;
 double pdotn,ddotn,mu,d2,Rsky;
 long ix;
 if(R_Nground > 0 && fabs(ddotn=dot(Ground.n,d)) > 1.e-3){
   VECSUB(Ground.p,p,dp)
   pdotn=DOT(Ground.n,dp);
   mu=pdotn/ddotn;
   if(mu > 0.0){
     vecscale(mu,d,gp);
     vecsum(p,gp,gp);
     GetGroundLight(gp,&colour[0],&colour[1],&colour[2]);
   }
   else goto NEXTBLOCK;
 }
 else{
   NEXTBLOCK:
   if(Sky.type == SKYGRADED){
     d2=(double)(Ycentre-row);
     d2/=scaley;
     //Rsky=fabs((CamAlpha+atan2(d2,1.000))/PI*2.0);  // Replaced to match OpenGL
     Rsky=fabs((double)(Ycentre-row))/(double)Ycentre;     
     colour[0]=((double)Sky.Zenith[0]-(double)Sky.Horizon[0])*Rsky
               +(double)Sky.Horizon[0];
     colour[1]=((double)Sky.Zenith[1]-(double)Sky.Horizon[1])*Rsky
               +(double)Sky.Horizon[1];
     colour[2]=((double)Sky.Zenith[2]-(double)Sky.Horizon[2])*Rsky
               +(double)Sky.Horizon[2];
   }
   else if(Sky.type == SKYMAPPED && Sky.mapped){  // this looks like a cylindrical map
     ddotn=d[2];
     if(fabs(ddotn) > 1.e-2){
       vecsub(Sky.map.p,p,dp);
       pdotn=dp[2];
       mu=pdotn/ddotn;
       vecscale(fabs(mu),d,gp);
       vecsum(p,gp,gp);
       MapFromProjection(&Sky.map,TILE,EOPAQUE,gp,&colour[0],&colour[1],&colour[2]);
     }
   }
   else if(Sky.type == BACKDROP && Sky.mapped){
     ix=((row*Sky.map.ymax)/YMAX)*Sky.map.xmax
       +(col*Sky.map.xmax)/XMAX;
     colour[0] = Sky.map.p24R[ix];
     colour[1] = Sky.map.p24G[ix];
     colour[2] = Sky.map.p24B[ix];
   }
   else if(Sky.type == SKYCUBE && Sky.mapped){ // use environment map 
     unsigned char R,G,B;
     MapEnvironment(col,row,&R,&G,&B, &(Sky.map));
     colour[0]=(double)R; colour[1]=(double)G; colour[2]=(double)B;
   }
   else{
     colour[0] = (double)Sky.Zenith[0];
     colour[1] = (double)Sky.Zenith[1];
     colour[2] = (double)Sky.Zenith[2];
   }
 }
 return;
}

void GetMirrorValues(double *color, vector p, vector din, vector n,
                     int at_o, face *at_f,
                     long trace_reflection_depth,
                     long trace_refraction_depth,
                     long trace_refractive_flag){
  /* p    is point on mirror surface where ray hits   */
  /* din  is incident direction                       */
  vector d,dp,gp,di;
  double pdotn,ddotn,mu,ld;
  ddotn = dot(din,n);
  vecscale(2*ddotn,n,d);  /* d is reflected direction (same length as din) */
  vecsub(din,d,d);        /* d=2(din.n)n-din                               */
  if(R_Nground > 0){
    ddotn=dot(Ground.n,d);
    if(fabs(ddotn) > 1.e-2){
      vecsub(Ground.p,p,dp);
      pdotn=dot(Ground.n,dp);
      mu=pdotn/ddotn;
      if(mu > 0.0){
        vecscale(mu,d,di);
        vecsum(p,di,gp);
        GetGroundLight(gp,&color[0],&color[1],&color[2]);
        if(trace_reflections == 1){
          ld = (double)((di[0]*di[0]) + (di[1]*di[1]) + (di[2]*di[2]));
          if(fabs(ld) > 1.e-10){
            ld=sqrt(ld);
            VECSCALE((1.0/ld),di,di)
            trace_reflection_ray(p,di,ld,color,
                                 at_o,at_f,
                                 trace_reflection_depth,
                                 trace_refraction_depth,
                                 trace_refractive_flag);
          }
        }
        return;
      }
    }
  }
  color[0] = (double)Sky.Zenith[0];
  color[1] = (double)Sky.Zenith[1];
  color[2] = (double)Sky.Zenith[2];
  if(Sky.type == SKYMAPPED){
    ddotn=d[2];
    if(fabs(ddotn) > 1.e-2){
      vecsub(Sky.map.p,p,dp);
      pdotn=dp[2];
      mu=pdotn/ddotn;
      vecscale(fabs(mu),d,gp);
      vecsum(p,gp,gp);
      /* check for existance of sky map elsewhere */
      MapFromProjection(&Sky.map,TILE,EOPAQUE,gp,&color[0],&color[1],&color[2]);
    }
  }
  else if(Sky.type == SKYCUBE){
    MapEnvironmentInDirection(d,&color[0],&color[1],&color[2],&Sky.map);
  }
  if(trace_reflections == 1){
    ld= (double)((d[0]*d[0]) + (d[1]*d[1]) + (d[2]*d[2]));
    if(fabs(ld) > 1.e-10){
      ld=sqrt(ld);
      VECSCALE((1.0/ld),d,di)
      trace_reflection_ray(p,di,
                           (double)BIGP,color,
                           at_o,at_f,   
                           trace_reflection_depth,
                           trace_refraction_depth,
                           trace_refractive_flag);
    }
  }
  return;
}

static void AttenuateGlass(long row, long col,
                    UCHAR *Red, UCHAR *Green, UCHAR *Blue,
                    long *Zglass[], long *ZglassO[], double *Zglassd[],
                    double *Zposn){
 long i;                          /* if the face on top is glass */
 for(i=NGLASS-1;i>=0;i--){        /* add in its surface attribs  */
   if(Zglass[i][col] >= 0)AttenuateGlassSurface(i,row,col,Red,Green,Blue,
                                                Zglass,ZglassO,Zglassd,
                                                Zposn);
 }
}

static void AttenuateGlassSurface(long id, long row, long col,
                         UCHAR *Red, UCHAR *Green, UCHAR *Blue,
                         long *Zglass[], long *ZglassO[], double *Zglassd[],
                         double *Zposn){
 double tr,tr1,srf=0.0,srf1=1.0;
 double Iab,IdR,IdG,IdB,IsR,IsG,IsB,c,cc,det,alpha,beta,bend,bend1,falloff;
 double pdotn,ldotn;
 double colorb0=0,colorb1=0,colorb2=0,colorr0=0,colorr1=0,colorr2=0,
        colors[3],colorm[3];
 vector n,p,pr,pl,bv;
 long j,obj;
 short clear,type;
 face    *f;
 vertex  *v;
 imap    *Map;
 matl    *Mat;
 BOOL    bSmooth,xtext=FALSE;
 unsigned char mattr,matre;
 obj=ZglassO[id][col];
 f=Object[obj].Fbase+Zglass[id][col];
 v=Object[obj].Vbase;
 if(f->ffmap >= 0){
   if((long)f->ffmap >= Object[obj].NoMaps)Map=NULL;
   else Map=(Object[obj].Mbase+(long)f->ffmap);
 }
 else Map=NULL;
 if(f->ffmat < 0 || f->ffmat >= Object[obj].NoMats)Mat=NULL;
 else Mat=(Object[obj].Tbase+(long)f->ffmat);

 if(Mat != NULL)bSmooth=Mat->bSmooth;
 else           bSmooth=f->bSmooth;
 if(!bSmooth)VECCOPY(f->n,n)  /* phong = NO */
 else{
   det=(
       (v[f->V[2]].y - v[f->V[0]].y) * (v[f->V[1]].x - v[f->V[0]].x)
      -(v[f->V[2]].x - v[f->V[0]].x) * (v[f->V[1]].y - v[f->V[0]].y)
       );
   if(det == 0){
     VECCOPY(f->n,n)
   }
   else{
     det=1.0/det;
     alpha= (double)(
        (v[f->V[2]].y - v[f->V[0]].y) * (Zposn[col] - v[f->V[0]].x)
       -(v[f->V[2]].x - v[f->V[0]].x) * ((double)row - v[f->V[0]].y)
                    )*det;
     beta= (double)(
        ((double)row - v[f->V[0]].y) * (v[f->V[1]].x - v[f->V[0]].x)
       -(Zposn[col] - v[f->V[0]].x) * (v[f->V[1]].y - v[f->V[0]].y)
                   )*det;
     n[0]=f->pn[0][0] + alpha*(f->pn[1][0] - f->pn[0][0])
                     +  beta*(f->pn[2][0] - f->pn[0][0]);
     n[1]=f->pn[0][1] + alpha*(f->pn[1][1] - f->pn[0][1])
                     +  beta*(f->pn[2][1] - f->pn[0][1]);
     n[2]=f->pn[0][2] + alpha*(f->pn[1][2] - f->pn[0][2])
                     +  beta*(f->pn[2][2] - f->pn[0][2]);
     normalize(n);
   }
 }
 /* p is the point on the glass face that is seen at this pixel */
 p[0]=(((double)col-(double)Xcentre)/scalex*Zglassd[id][col]);
 p[1]=Zglassd[id][col];
 p[2]=(((double)Ycentre-(double)row)/scaley*Zglassd[id][col]);
 if((pdotn=DOT(p,f->n)) > 0){n[0] = -n[0]; n[1] = -n[1]; n[2] = -n[2];}
 IdR=IdG=IdB=0.0; IsR=IsG=IsB=0.0;
 if(Map != NULL && Map->bp >= 0){
   if(Map->bTiled)type=TILE;
   else           type=SINGLE;
   if(Map->map == MAP_BY_VERTEX)
     BumpByCoordinates(Map,type,obj,f,p,n);
   else{;}
 }
 if(Mat != NULL){
   VECCOPY((double)Mat->color,colors)
   tr=((double)(Mat->transp))*3.92e-3;
   if(Mat->bInternalShader)xtext=FALSE; else xtext=TRUE;
 }
 else{
   VECCOPY((double)f->color,colors) 
   tr=((double)(8))*3.92e-3;   // shoudn't happen as only materials can be glass
   xtext=FALSE;
 }
 if(xtext){ /* external shader */
   double alpha_channel=0.0;
   int s;
   double beta,gamma;
   MapTextureByCoordinates(obj,f,p,&alpha,&beta,&gamma);
   for(s=0;s<4;s++){
     if(Mat->shaderID[s] >= 0)
       RenderExternalTexture(Mat->shaderID[s],
           p,n,
           alpha,beta,gamma,
           (double)(Mat->params[0])*SHSCALE,  
           &alpha_channel,
           &tr,
           Mat->texco,
           obj,f,
           colors);
   }
 }
 if(Mat != NULL && Mat->bGlassEdge){       /* glass edge is on */
   if((bend=DOT(p,n)) < 0)bend = -bend;
   bend1=sqrt(bend/R_length(p));
   tr1  = 1.0-tr;
   tr *= bend1;                 /* to give glass some edge */
 }
 else{
   tr1  = 1.0-tr;
 }
 GetPixelLighting(Mat,p,n,pdotn,0.0,obj,f,
                  &IdR,&IdG,&IdB,&IsR,&IsG,&IsB);
 if(Map != NULL){
   if(Map->rp >= 0){
     MapMirror(&colorr0,&colorr1,&colorr2,Map,p,n);
     srf=Map->rd; srf1=(1.0-srf);
   }
   if(Map->pp < 0)goto MISSBRUSH;
   if(Map->bDecal)clear = TRANSP;
   else           clear = EOPAQUE;
   if(Map->map == MAP_BY_VERTEX){
     if(Map->bTiled)type = TILE;
     else           type = SINGLE;
     if(MapByCoordinates(Map,obj,f,type,clear,p,&colorb0,&colorb1,&colorb2)
        == 0)goto MISSBRUSH1;
   }
   else{;}
   if(!(Map->bShaded)){       /* unshaded map required */
     IdR=IdG=IdB=1.0; IsR=IsG=IsB=0.0;
   }
   tr1=Map->pd; tr=(1.0-tr1);    /* tr is 0-1 % of transparency */
   *(colors  ) = colorb0;
   *(colors+1) = colorb1;
   *(colors+2) = colorb2;
   MISSBRUSH:
   if(Map->tp >= 0){ /* transparency map - white is fully opaque */
     if(Map->bTiled)type = TILE;
     else           type = SINGLE;
     if(Map->map == MAP_BY_VERTEX){
       if(MapTransByCoordinates(Map,type,obj,f,p,&tr) != 0)tr1=(1.0-tr);
     }
     else{;}
   }
   MISSBRUSH1:
   IdR=(double)Red[col]*tr   + (colors[0]*IdR+IsR*255)*tr1;
   IdG=(double)Green[col]*tr + (colors[1]*IdG+IsG*255)*tr1;
   IdB=(double)Blue[col]*tr  + (colors[2]*IdB+IsB*255)*tr1;
   IdR=IdR*srf1 + colorr0*srf;
   IdG=IdG*srf1 + colorr1*srf;
   IdB=IdB*srf1 + colorr2*srf;
 }
 else{
   IdR=(double)Red[col]*tr   + colors[0]*tr1*IdR+IsR*255;
   IdG=(double)Green[col]*tr + colors[1]*tr1*IdG+IsG*255;
   IdB=(double)Blue[col]*tr  + colors[2]*tr1*IdB+IsB*255;
   if(Mat != NULL)matre=Mat->refl; else matre=0;
   if(matre > 0){
     GetMirrorValues(colorm,p,p,n,0,NULL,0,0,0);
     srf = (double)(matre)/255.0; srf1=(1.0-srf);
     IdR=IdR*srf1 + colorm[0]*srf;
     IdG=IdG*srf1 + colorm[1]*srf;
     IdB=IdB*srf1 + colorm[2]*srf;
   }
 }
 Red[col]  =(unsigned char)min(255.0,IdR);
 Green[col]=(unsigned char)min(255.0,IdG);
 Blue[col] =(unsigned char)min(255.0,IdB);
}

void GetPixelValue(vector p, vector dinput, // p is point on surface  dinput is input direction
                   vector n, double pdotn, 
                   imap *Map, matl *Mat, face *f, long obj,
                   long row, long col, BOOL first_trace,
                   long trace_relection_depth, long trace_refraction_depth, long refractive_flag,
                   double *Red, double *Green, double *Blue){
 BOOL     xtext,fullsat=FALSE;
 vector   pr,bv,pl,d;
 long     j,flag_mirror=0,t,flag_trans=0,refx=1000;
 short    type,textID,clear;
 double   alpha,beta,gamma;
 double   colors[3],colorr[3],colorm[3],colort[3];
 double   Iab,IdR,IdG,IdB,IsR,IsG,IsB,c,cc,
          ldotn,amblight,falloff;
 double   sf,sf1,tr,tr1,rf,rf1;
 unsigned char mattr,matre;
 
 IdR=IdG=IdB=0.0; IsR=IsG=IsB=0.0;
 colorr[0]=colorr[1]=colorr[2]=0.0;
 colorm[0]=colorm[1]=colorm[2]=0.0;
 colort[0]=colort[1]=colort[2]=0.0;

 if(Mat != NULL){
   VECCOPY((double)Mat->color,colors)
   if(Mat->bInternalShader){textID=Mat->internalShaderID; xtext=FALSE;}
   else                    {textID=0; xtext=TRUE;}
   tr=((double)(Mat->transp))*3.92e-3;
   refx=Mat->refractive_index;
 }
 else{
   VECCOPY((double)f->color,colors)
   textID=0; 
   xtext=FALSE;
   tr=0;
   refx=1000;
 }
 sf=1.0;
 /* n is to be normal to surface at that point p */
 if(Map != NULL && Map->bp >= 0){
   if(Map->bTiled)type=TILE;
   else           type=SINGLE;
   if(Map->map == MAP_BY_VERTEX)
     BumpByCoordinates(Map,type,obj,f,p,n);
   else{;}
 }
 if(textID > 0 || xtext){
   MapTextureByCoordinates(obj,f,p,&alpha,&beta,&gamma);
   if(xtext){ /* external shader */
     double alpha_channel=0.0;
     int s; 
     for(s=0; s<4;s++){
       if(Mat->shaderID[s] >= 0)
         fullsat=RenderExternalTexture(Mat->shaderID[s],
           p,n,
           alpha,beta,gamma,
           (double)(Mat->params[0])*SHSCALE,  
           &alpha_channel,
           &tr,
           Mat->texco,
           obj,f,
           colors);
     }
   }
   else fullsat=GetInternalTexture(Mat,alpha,beta,gamma,textID,p,n,obj,f,colors,
                     colort,&flag_trans,tr,&tr1,
                     trace_relection_depth,trace_refraction_depth,
                     refractive_flag);
 }
 GetPixelLighting(Mat,p,n,pdotn,0.0,obj,f,
                  &IdR,&IdG,&IdB,&IsR,&IsG,&IsB);
 colorm[0]=colorm[1]=colorm[2]=0.0;
 if(Map != NULL){
   if(Map->rp >= 0){
     MapMirror(&colorr[0],&colorr[1],&colorr[2],Map,p,n);
     flag_mirror=1; rf=Map->rd;
   }
   if(Map->tp >= 0){ /* transparency map - white is fully transparent */
     if(Map->bTiled)type = TILE;
     else           type = SINGLE;
     if(MapTransByCoordinates(Map,type,obj,f,p,&tr) != 0)tr1=(1.0-tr);
     if(trace_refractions && tr > 0.0){
       c=(double)((dinput[0]*dinput[0]) + (dinput[1]*dinput[1]) + (dinput[2]*dinput[2]));
       c=sqrt(c);
       VECSCALE((1.0/c),dinput,d)
       if(first_trace && refx == 1000){
         GetBackdropValue(row,col,p,d,colort);
       }
       else {
         colort[0]=*Red; colort[1]=*Green; colort[2]=*Blue;
       }
       if(refx == 1000)trace_transmission_ray(p,d,n,obj,f,(double)BIGP,colort, 
         trace_relection_depth,trace_refraction_depth,refractive_flag);
       else trace_refraction_ray(p,d,n,obj,f,(double)BIGP,colort,refx,
         trace_relection_depth,trace_refraction_depth,refractive_flag);
       flag_trans=1;
     }
   }
   if(Map->pp < 0)goto ALGOTEXT;
   if(Map->bDecal)clear = TRANSP;
   else           clear = EOPAQUE;
   if(Map->map == MAP_BY_VERTEX){
     if(Map->bTiled)type = TILE;
     else           type = SINGLE;
     if(MapByCoordinates(Map,obj,f,type,clear,p,
                         &colorm[0],&colorm[1],&colorm[2])
        == 0)goto ALGOTEXT;
     flag_trans=2;
   }
   if(Map->bShaded){       /* shaded map required */
     colorm[0] = colorm[0]*IdR+IsR*255;        /* put spec after mix  */
     colorm[1] = colorm[1]*IdG+IsG*255;
     colorm[2] = colorm[2]*IdB+IsB*255;
   }
   sf=(1.0-Map->pd);          /* map "pd" is 0.0-1.0  or pp */
   if(Map->pp == 100){        /* make sure only map is used */
     colors[0]=colors[1]=colors[2]=0.0;
     goto SKIPALGO;
   }
 }
 /* algorithmic shaders */
 ALGOTEXT:
 if(fullsat){IdR=IdG=IdB=1.0; IsR=IsG=IsB=0.0;}

 SKIPALGO:
 /* combine the shader and map with appropriate illimination */
 sf1=1.0-sf;
 colors[0] = colors[0]*IdR*sf + colorm[0]*sf1;
 colors[1] = colors[1]*IdG*sf + colorm[1]*sf1;
 colors[2] = colors[2]*IdB*sf + colorm[2]*sf1;
 /* colors now is the surface value including diffuse lighting */
 if(trace_shadows && Object[obj].show_shadow &&
    (j=check_for_shadow(obj,p,f->n)) != 0){
   IsR = IsG = IsB = 0.0;
   if(j < 0){
     colors[0] *= shadow_density;     /* do the shadows */
     colors[1] *= shadow_density;
     colors[2] *= shadow_density;
   }
   else{ /* j layers of glass */
     colors[0] *= glass_shadow_density;
     colors[1] *= glass_shadow_density;
     colors[2] *= glass_shadow_density;
   }
 }
 if(Mat != NULL)mattr=Mat->transp;
 else           mattr=0;
 if(trace_refractions && mattr > 0 && flag_trans == 0){ 
   double bend,bend1;
   BOOL glass_edge;
   if(Mat != NULL)glass_edge=Mat->bGlassEdge;
   else           glass_edge=FALSE;
     if(glass_edge){       /* glass edge is on */
     if((bend=DOT(p,n)) < 0)bend = -bend;
     bend1=sqrt(bend/R_length(p));
     tr1  = 1.0-tr;
     tr *= bend1;                 /* to give glass some edge */
   }
   else{
     tr1  = 1.0-tr;
   }
   c=(double)((dinput[0]*dinput[0]) + (dinput[1]*dinput[1]) + (dinput[2]*dinput[2]));
   c=sqrt(c);
   VECSCALE((1.0/c),dinput,d)
   if(refx == 1000){ // This not refractive glass
     GetBackdropValue(row,col,p,d,colort);
     trace_transmission_ray(p,d,n,obj,f,(double)BIGP,colort,
         trace_relection_depth,trace_refraction_depth,refractive_flag);
  }
  else  trace_refraction_ray(p,d,n,obj,f,(double)BIGP,colort,refx,
                        trace_relection_depth,trace_refraction_depth,refractive_flag);
   flag_trans=1;
 }
 if(flag_trans == 1){  /* mix any transmissivity */
   colors[0] = colort[0]*tr + colors[0]*tr1;
   colors[1] = colort[1]*tr + colors[1]*tr1;
   colors[2] = colort[2]*tr + colors[2]*tr1;
 }
 colors[0] += IsR*255; /* mix the specular highlights */
 colors[1] += IsG*255;
 colors[2] += IsB*255;

 if(Mat != NULL)matre=Mat->refl;
 else           matre=0;
 if(matre > 0 && ! flag_mirror){    /* mirror added to finish */
   GetMirrorValues(colorr,p,dinput,n,obj,f,
      trace_relection_depth,trace_refraction_depth,refractive_flag);
   rf = (double)(matre)/255.0;
   flag_mirror=1;
 }

 if(flag_mirror){
   rf1=1.0-rf;
   *Red  = (unsigned char)min(255,colorr[0]*rf+colors[0]*rf1);
   *Green= (unsigned char)min(255,colorr[1]*rf+colors[1]*rf1);
   *Blue = (unsigned char)min(255,colorr[2]*rf+colors[2]*rf1);
 }
 else{
   *Red  = (unsigned char)min(255,colors[0]);
   *Green= (unsigned char)min(255,colors[1]);
   *Blue = (unsigned char)min(255,colors[2]);
 }
}


void GetRowPixelValue(long row, long col, imap *Map, matl *Mat,
                   face *f, long obj, vector nn,  /* obj is object id */
                   UCHAR *Red, UCHAR *Green, UCHAR *Blue,
                   double *Zdepth,UCHAR *ZobjChar){
 double r,g,b,pdotn;
 vector p,n;
 /* p is the point on the face that is seen at this pixel */
 p[0]=(((double)col-(double)Xcentre)/scalex*Zdepth[col]);
 p[1]=Zdepth[col];
 p[2]=(((double)Ycentre-(double)row)/scaley*Zdepth[col]);
 if((pdotn=DOT(p,f->n)) < 0)VECCOPY(nn,n)
 else     {n[0] = -nn[0]; n[1] = -nn[1]; n[2] = -nn[2];}
 if(trace_refractions){  // if transparent must restart transparent ray because ray does not match Z buffer
   double maptr,mattr;
   if(Map != NULL)maptr=Map->tp;
   else           maptr=0;
   if(Mat != NULL)mattr=Mat->transp;
   else           mattr=0;
   if(mattr > 0.0 || maptr > 0.0){ // Transparent so start ray back at origin
     int fc;
     vector d; 
     d[0]=(((double)col-(double)Xcentre+0.5)/scalex);
     d[1]=1.0;
     d[2]=(((double)Ycentre-(double)row)/scaley);
     normalize(d);
     p[0]=p[1]=p[2]=0.0;
     if(!trace_starting_ray(d,p,&obj,&fc,&f))return;   // p is output
     p[0]=p[1]=p[2]=0.0;   // Get Surface needs (0,0,0) to start 
     GetSurfaceValue(obj,f,Object[obj].Vbase,p,d,&r,&g,&b,0,0,0);
     goto ENDP;
   } 
 }
 GetPixelValue(p,p,n,pdotn,Map,Mat,f,obj,row,col,TRUE,0,0,0,&r,&g,&b);
 ENDP: 
 Red[col]  = (unsigned char)min(255,r);
 Green[col]= (unsigned char)min(255,g);
 Blue[col] = (unsigned char)min(255,b);
}

void GetSurfaceValue(int obj, face *f, vertex *v,
                     vector p, vector dr, double *r, double *g, double *b,
                     long trace_reflection_depth, long trace_refraction_depth,
                     long refractive_flag){
/* used by ray tracing and reflection               */
/* p is origin point of ray, dr is direction of ray */
/* face ID is known - return RGB                    */
 imap *Map;
 matl *Mat;
 vector d,n,pr,pl,dp,bv,dl,nn;
 double ddotn,mu,pdotn,alpha,beta;
 BOOL bSmooth;
 ddotn=DOT(f->n,dr);
 if(fabs(ddotn) > 1.e-3){/* vector will hit face */
   VECSUB(v[f->V[0]].p,p,dp)
   pdotn=DOT(f->n,dp);
   mu=pdotn/ddotn;
   VECSCALE(mu,dr,dp)          /* dp is vector from ray origin      */
   VECSUM(p,dp,d)              /* d is the point on visible face    */
   if(f->ffmat < 0 || f->ffmat >= Object[obj].NoMats)Mat=NULL;
   else Mat=(Object[obj].Tbase+(long)f->ffmat);
   if(Mat != NULL)bSmooth=Mat->bSmooth;
   else           bSmooth=f->bSmooth;
   if(!bSmooth)veccopy(f->n,n);
   else{                       /* phong smoooth */
     VECSUB(v[f->V[1]].p,v[f->V[0]].p,pr)
     VECSUB(v[f->V[2]].p,v[f->V[0]].p,pl)
     hitpoint(f->n,pr,pl,v[f->V[0]].p,d,&alpha,&beta);
     VECSUB(f->pn[1],f->pn[0],pr)
     VECSUB(f->pn[2],f->pn[0],pl)
     VECSCALE(alpha,pr,pr)
     VECSCALE(beta,pl,pl)
     VECSUM(f->pn[0],pr,n)
     VECSUM(n,pl,n)
     normalize(n);
   }
   if(f->ffmap >= 0){ /* choose map */
      if((long)f->ffmap >= Object[obj].NoMaps)Map=NULL;
      else Map=(Object[obj].Mbase+(long)f->ffmap);
   }
   else Map=NULL;
   // added to allow for double sided faces in reflection and trace
   if(DOT(dp,n) < 0.0)VECCOPY(n,nn)
   else    {nn[0] = -n[0]; nn[1] = -n[1]; nn[2] = -n[2];}
   // otherwise just use "n" in arguments below
   GetPixelValue(d,dr,nn,-1.0,Map,Mat,f,obj,
                 0,0,FALSE,        // NOT first hit so no backdrop to get   
                 trace_reflection_depth,trace_refraction_depth,refractive_flag,
                 r,g,b); // returned values.
 }
}