REFLECT.C

Go to the documentation of this file.

// file reflect.c

#define MODULE_MIRRORS

/* The functions in this file get surface colours for an reflection vector and */
/* the generation of the ground reflection. ground reflection is done by using */
/* a full Z screen buffer to record the identity of the faces that would be    */
/* seen in the ground mirrors. Then when a ground pixel needs to get a         */
/* reflection it is calculated from the reflection ray given that we know the  */
/* id of the face involved. This is quicker that calculating a reflection image*/
/* and using that as a map - I think.                                          */
/* They also implement the Spotlight shadows - which is still NOT satisfactory */

#include "render.h"

static long IntersectMirror(long scanline, double v1x, double v1y,
                     double v2x, double v2y,  double *xi);
static short ShadowMirrorLine(long scanline, double vx[3], double vy[3],
                      double *leftedge, double *rightedge,
                      double xmin, double xmax);
static void UpdateShadowZbuffer(float *id,float *id2,
                                long l, double xl, double xr,
                                double d, double d2, vector en,
                                double ec, double SBC);

extern double spec_pow[16][max_spec+1],Rip_Freq;

short mirror_ground=0;

//#define US unsigned short
#define US UNSIGNED
#define UC unsigned char

static UC     *grf_o;
static US     *grf_f;
static float  *grf_z;
static long   shadow_buffer_size2;

#ifdef _SUNSTYLE
static long IntersectMirror(scanline,v1x,v1y,v2x,v2y,xi)
long scanline; double v1x; double v1y; double v2x; double v2y; double *xi; {
#else
static long IntersectMirror(long scanline, double v1x, double v1y,
                     double v2x, double v2y,  double *xi){
#endif
 double ddx,ddy,doff,ddxy;    /* scanline function for mirror */
 if((v1y < scanline)  && (v2y < scanline))return 0;
 if((v1y > scanline)  && (v2y > scanline))return 0;
 ddy = (v2y) - (v1y);
 if(fabs(ddy) < 0.0001)return 0;
 ddx = (v2x) - (v1x);
 doff=(double)scanline - v1y;
 ddxy=ddx*doff/ddy;
 *xi=(v1x)+ddxy;
 return 1;
}

#ifdef _SUNSTYLE
static short ShadowMirrorLine(scanline,vx,vy,leftedge,rightedge,xmin,xmax)
long scanline; double vx[3]; double vy[3]; double *leftedge;
double *rightedge; double xmin; double xmax; {
#else
static short ShadowMirrorLine(long scanline, double vx[3], double vy[3],
                      double *leftedge, double *rightedge,
                      double xmin, double xmax){
#endif
 long i1,i2,i3;                /* scanline function for mirror */
 double x1,x2,x3;
 i1=IntersectMirror(scanline,vx[0],vy[0],vx[1],vy[1],&x1);
 i2=IntersectMirror(scanline,vx[1],vy[1],vx[2],vy[2],&x2);
 i3=IntersectMirror(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((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 >= xmax)return 0;
 else *leftedge  = max(0,*leftedge);
 if(*rightedge < xmin)return 0;
 *rightedge = min(xmax-0.01,*rightedge);
 return 1;
}

static void UpdateZbuffer(long l, double xl, double xr, long id, long O,
                          double d2, vector en, double ec,
                          long Mxmax, long Mxcentre, double Mscalex){
 long i,ixl,ixr;               /* scanline function for mirror */
 double depth,bL,bR,aL,aR,dL,dR;
 xl += 0.01;
 xr -= 0.01;
 ixl = (long)xl;
 ixr = (long)xr;
 if(ixl == ixr){
   dL = xl - (double)Mxcentre; dL/=Mscalex;
   dR = xr - (double)Mxcentre; dR/=Mscalex;
   bL = en[0]*dL + en[2]*d2 + en[1];
   if(fabs(bL) < 1.e-30)aL=BIGP;
   else aL = ec/bL;
   bR = en[0]*dR + en[2]*d2 + en[1];
   if(fabs(bR) < 1.e-30)aR=BIGP;
   else aR = ec/bR;
   depth=min(aL,aR);
   if(depth > 1.000 && depth <= *(grf_z+l*(long)Mxmax+ixl)){
       *(grf_z+l*(long)Mxmax+ixl) = depth;
       *(grf_o+l*(long)Mxmax+ixl) = (UC)O;
       *(grf_f+l*(long)Mxmax+ixl) = (US)id;
   }
 }
 else{
   dL = xl - (double)Mxcentre; dL/=Mscalex;
   bL = en[0]*dL + en[2]*d2 + en[1];
   if(fabs(bL) < 1.e-30)aL=BIGP;
   else aL = (ec)/bL;
   bR = en[0]*DR[ixl] + en[2]*d2 + en[1];
   if(fabs(bR) < 1.e-30)aR=BIGP;
   else aR = ec/bR;
   depth=min(aL,aR);
   if(depth > 1.000 && depth <= *(grf_z+l*(long)Mxmax+ixl)){
     *(grf_z+l*(long)Mxmax+ixl) = depth;
     *(grf_o+l*(long)Mxmax+ixl) = (UC)O;
     *(grf_f+l*(long)Mxmax+ixl) = (US)id;
   }
   bL = en[0]*DL[ixr] + en[2]*d2 + en[1];
   if(fabs(bL) < 1.e-30)aL=BIGP;
   else aL = ec/bL;
   dR = xr - (double)Mxcentre; dR/=Mscalex;
   bR = en[0]*dR + en[2]*d2 + en[1];
   if(fabs(bR) < 1.e-30)aR=BIGP;
   else aR = (ec)/bR;
   depth=min(aR,aL);
   if(depth > 1.000 && depth <= *(grf_z+l*(long)Mxmax+ixr)){
     *(grf_z+l*(long)Mxmax+ixr) = depth;
     *(grf_o+l*(long)Mxmax+ixr) = (UC)O;
     *(grf_f+l*(long)Mxmax+ixr) = (US)id;
   }
   if(ixr-ixl  > 1)for(i=ixl+1;i<ixr;i++){
     bL = en[0]*DL[i] + en[2]*d2 + en[1];
     if(fabs(bL) < 1.e-30)aL=BIGP;
     else aL = ec/bL;
     bR = en[0]*DR[i] + en[2]*d2 + en[1];
     if(fabs(bL) < 1.e-30)aL=BIGP;
     else aR = ec/bR;
     depth=min(aL,aR);
     if(depth > 1.00 && depth <= *(grf_z+l*(long)Mxmax+i)){
       *(grf_z+l*(long)Mxmax+i) = depth;
       *(grf_o+l*(long)Mxmax+i) = (UC)O;
       *(grf_f+l*(long)Mxmax+i) = (US)id;
     }
   }
 }
}

void Free_Mirror_On_Ground(void){
 if(mirror_ground == 0)return;
 X__Free(grf_o); /* z buffer is X__Freed after it is used */
 X__Free(grf_f);
 mirror_ground=0;
}

short Mirror_In_Ground(void){       
 long i,j,l,fbc;
 long vi[3],screensize,Mxmax,Mymax,Mxcentre,Mycentre;
 face    *f;
 vertex  *v;
 vector p,pp[3],en;
 double dp,temp4,d2,ec,xx[3],yy[3],xxmin,xxmax,yymin,yymax,Mscalex,Mscaley;
 double itr[4][4],tr[4][4],t1[4][4],t2[4][4],t3[4][4],x,y,z,xg,yg,zg;
 Mscalex=scalex; Mscaley=scaley;
 if(full_mirror){
   Mxmax=XMAX; Mymax=YMAX;
 }
 else{
   Mxmax=ResolutionX; Mymax=ResolutionY;
   if(anti_alias == LOW)   {Mscalex /= 2;}
   if(anti_alias == MEDIUM){Mscalex /= 2; Mscaley /= 2;}
   if(anti_alias == HIGH)  {Mscalex /= 3; Mscaley /= 2;}
   if(anti_alias == ULTRA) {Mscalex /= 3; Mscaley /= 3;}
 }
 Mxcentre=Mxmax/2; Mycentre=Mymax/2;
 screensize = Mxmax*Mymax;
 if((grf_o = (UC *)X__Malloc(screensize*sizeof(UC))) == NULL){
   Render_Message(1,0,NULL); return FAIL;
 }
 if((grf_f = (US *)X__Malloc(screensize*sizeof(US))) == NULL){
   X__Free(grf_o);
   Render_Message(1,0,NULL); return FAIL;
 }
 if((grf_z = (float *)X__Malloc(screensize*sizeof(float))) == NULL){
   X__Free(grf_o);
   X__Free(grf_f);
   Render_Message(1,0,NULL); return FAIL;
 }
 mirror_ground=1;
 for(i=0;i<screensize;i++){
   *(grf_o + i) =  255;
   *(grf_f + i) =  65535;
   *(grf_z + i) =  BIGP;
 }
 R_roty(t1,CamTheta);   /* Get needed transformations          */
 R_rotx(t2,CamAlpha);
 R_m4by4(t2,t1,t3);
 R_rotz(t1,CamPhi);
 R_m4by4(t1,t3,t2);
 R_tram(t1,ViewPoint[0],ViewPoint[1],ViewPoint[2]);
 R_m4by4(t1,t2,itr);   /* ITR is inverse viewing transformation */
 R_tram(t2,-ViewPoint[0],-ViewPoint[1],-ViewPoint[2]);
 R_rotz(t3,-CamPhi);
 R_m4by4(t3,t2,t1);
 R_rotx(t2,-CamAlpha);
 R_m4by4(t2,t1,t3);
 R_roty(t1,-CamTheta);
 R_m4by4(t1,t3,tr);    /* TR is the viewing transformation      */
 R_m4by1(itr,Ground.p[0],Ground.p[1],Ground.p[2],&xg,&yg,&zg);
 for(i=0;i<Mxmax;i++) {
   DL[i]=(double)(i - Mxcentre);        DL[i]/=Mscalex;
   DR[i]=(double)(i+1-Mxcentre);        DR[i]/=Mscalex;
   D1[i]=(double)(i - Mxcentre) + 0.5;  D1[i]/=Mscalex;
 }
 for(i=0;i<ObjectCount;i++){
   if(!Object[i].in_use)continue;
   v=Object[i].Vbase;
   f=Object[i].Fbase;
   for(j=0;j<Object[i].NoFaces;j++){
     vi[0]=(f+j)->V[0];
     vi[1]=(f+j)->V[1];
     vi[2]=(f+j)->V[2];
     fbc=0;  /* count vertices below ground */
     for(l=0;l<3;l++){
       R_m4by1(itr,(v+vi[l])->p[0],(v+vi[l])->p[1],(v+vi[l])->p[2],&x,&y,&z);
       dp=z-zg;                    /* delta point and mirror plane */
       if(dp < 0.0){fbc++; z = zg;} /*  goto FACEBREAK;  point below plane*/
       else   z = zg-dp;
       R_m4by1(tr,x,y,z,&p[0],&p[1],&p[2]);
       veccopy(p,pp[l]);
       if(p[1] >= 1.000){  /* vertex in front of projection plane */
         temp4=(Mscalex*(p[0])/(p[1]));
         xx[l]=(double)Mxcentre + temp4;
         temp4=(Mscaley*(p[2])/(p[1]));
         yy[l]=(double)Mycentre - temp4;
       }
       else goto FACEBREAK;
     }
     if(fbc == 3)goto FACEBREAK; /* all vertices below ground */
     xxmin=min(xx[1],xx[2]);     xxmin=min(xx[0],xxmin);
     xxmax=max(xx[1],xx[2]);     xxmax=max(xx[0],xxmax);
     yymin=min(yy[1],yy[2]);     yymin=min(yy[0],yymin);
     yymax=max(yy[1],yy[2]);     yymax=max(yy[0],yymax);
     if((yymax > yymin) && (xxmax > xxmin)){
       VECSUB(pp[2],pp[0],pp[2])
       VECSUB(pp[1],pp[0],pp[1])
       cross(pp[1],pp[2],en);
       if(normalize(en) == 0){ /* only for good normals */
         ec=DOT(pp[0],en);
         if(ec > 0.0){
           VECSCALE((double)(-1.0),en,en)
           ec = -ec;
         }
         for(l=max(0,(long)yymin);l<=min((long)yymax,Mymax-1);l++){
           if(ShadowMirrorLine(l,xx,yy,&xxmin,&xxmax,
                               0.0,(double)Mxmax)){
             d2=(double)(Mycentre-l)/Mscaley;
             UpdateZbuffer(l,xxmin,xxmax,j,i,d2,en,ec,Mxmax,Mxcentre,Mscalex);
           }
         }
       }
     }
     FACEBREAK: continue;
   }
 }
 X__Free(grf_z); grf_z=NULL;
 return 1;
}

void Add_Ground_Mirror(long scanline, long col, vector p,
       unsigned char *R, unsigned char *G, unsigned char *B){
/* p is vector to point on ground where reflection is */
 register long addr;
 int obj;
 face    *f;
 vertex  *v;
 double r1,pdotn,r=0.0,g=0.0,b=0.0;
 vector d;

 if(full_mirror)addr=XMAX*scanline+col;
 else addr=ResolutionX*(scanline/aaYstep)+col/aaXstep;

 obj=(int)(*(grf_o + addr));
 if( obj < 255){
   f = (Object + obj)->Fbase + (long)(*(grf_f + addr));
   v = (Object + obj)->Vbase;
   pdotn=DOT(p,Ground.n);
   VECSCALE(2*pdotn,Ground.n,d)
   VECSUB(p,d,d)                 /* vector in the reflected direction */
   GetSurfaceValue(obj,f,v,p,d,&r,&g,&b,0,0,0);
   r1=(1.0 - Ground.refl);
   *R = (unsigned char)
        ((double)*R * r1 + r*Ground.refl);
   *G = (unsigned char)
        ((double)*G * r1 + g*Ground.refl);
   *B = (unsigned char)
        ((double)*B * r1 + b*Ground.refl);
 }
 return;
}

void Free_Lights_Shadows(void){
 long i;
 for(i=0;i<Nlights;i++){
   if(Lights[i].shadow_buffer != NULL){
     X__Free(Lights[i].shadow_buffer);
     Lights[i].shadow_buffer=NULL;
   }
   if(Lights[i].object_buffer != NULL){
     X__Free(Lights[i].object_buffer);
     Lights[i].object_buffer=NULL;
   }
   if(Lights[i].face_buffer != NULL){
     X__Free(Lights[i].face_buffer);
     Lights[i].face_buffer=NULL;
   }
   if(Lights[i].atmospheric_shadow_buffer != NULL){
     X__Free(Lights[i].atmospheric_shadow_buffer);
     Lights[i].atmospheric_shadow_buffer=NULL;
   }
 }
}

void SetLightShadow(long id){  /* for light "id" */
 double itr[4][4],tr[4][4],t1[4][4],t2[4][4],t3[4][4],d,x,y,z,
        xx[3],yy[3],d2,ec,xxmin,xxmax,yymin,yymax;
 long i,j,l,vi[3];
 face   *f;
 vertex *v;
 vector en,pp[3],dr;
 float  d2f,*depth2=NULL;
 if(shadow_buffer_size == 0)return;
 if((Lights[id].shadow_buffer = (float *)X__Malloc((long)shadow_buffer_size*
     (long)shadow_buffer_size*(long)sizeof(float))) == NULL){
   Render_Message(2,0,NULL);
   return;
 }
 if((depth2 = (float *)X__Malloc((long)shadow_buffer_size*
     (long)shadow_buffer_size*(long)sizeof(float))) == NULL){
   X__Free(Lights[id].shadow_buffer);
   Lights[id].shadow_buffer=NULL;
   Render_Message(2,0,NULL);
   return;
 }
 else Render_Message(3,0,NULL);
 for(i=0;i< shadow_buffer_size*shadow_buffer_size; i++){
   *(Lights[id].shadow_buffer+i) = 1.0e30;
   *(depth2+i) = 1.e30;
 }
 R_roty(t1,CamTheta);
 R_rotx(t2,CamAlpha);
 R_m4by4(t2,t1,t3);
 R_rotz(t1,CamPhi);
 R_m4by4(t1,t3,t2);
 R_tram(t1,ViewPoint[0],ViewPoint[1],ViewPoint[2]);
 R_m4by4(t1,t2,itr);   /* ITR is inverse viewing transformation */
 R_m4by1(itr,Lights[id].p[0],Lights[id].p[1],Lights[id].p[2],&x,&y,&z);
 R_tram(t1,-x,-y,-z);
 R_rotz(t2,-Lights[id].dPhi*PI/180.0);
 R_m4by4(t2,t1,t3);
 R_rotx(t1,-Lights[id].dAlpha*PI/180.0);
 R_m4by4(t1,t3,t2);
 R_roty(t1,0.0);   /* no light bank dTheta*/
 R_m4by4(t1,t2,tr);  /* transform for light as viewing point */
 R_m4by4(tr,itr,Lights[id].shadow_transform); /* to move point to sh. buffer*/
 shadow_buffer_size2 = shadow_buffer_size/2.0;
 /* cone2 is 1/2 the cone angle + edge andle - radians */
 Lights[id].scale=d=(double)shadow_buffer_size2/tan(Lights[id].cone2);
 if(R_Nground > 0){ /* get ground depth */
   double c,b,mu;
   vector gp0,gp1,gp2,gn;
   R_m4by1(Lights[id].shadow_transform,
               Ground.p[0],Ground.p[1],Ground.p[2],&x,&y,&z);
   gp0[0]=x; gp0[1]=y; gp0[2]=z;
   R_m4by1(Lights[id].shadow_transform,
               Ground.x[0],Ground.x[1],Ground.x[2],&x,&y,&z);
   gp1[0]=x; gp1[1]=y; gp1[2]=z;
   R_m4by1(Lights[id].shadow_transform,
               Ground.y[0],Ground.y[1],Ground.y[2],&x,&y,&z);
   gp2[0]=x; gp2[1]=y; gp2[2]=z;
   vecsub(gp1,gp0,gp1);
   vecsub(gp2,gp0,gp2);
   cross(gp1,gp2,gn);
   normalize(gn);
   c=DOT(gp0,gn);
   for(i=0,l=0;i<shadow_buffer_size;i++){
     d2=(double)(shadow_buffer_size2-i);
     for(j=0;j<shadow_buffer_size;j++){ /* light points (0,1,0) */
//       dr[0]=((double)j + 0.5); dr[1]=d; dr[2]=d2; // For direct distance
//       normalize(dr); b=DOT(gn,dr);                // ""     ""
       b=gn[0]*((double)j + 0.5)+gn[1]*d+gn[2]*d2;
       if(fabs(b) > 1.e-4){
         mu=c/b;
         if(mu > 0.0){
           *(depth2+l) = mu*d;
//           *(depth2+l) = mu;    // direct distance
         }
         l++;
       }
     }
   }
 }
 for(i=0;i<ObjectCount;i++){
   if(!Object[i].in_use)continue;
   v=Object[i].Vbase;
   f=Object[i].Fbase;
   if(!Object[i].cast_shadow)continue;
   for(j=0;j<Object[i].NoFaces;j++){
     vi[0]=(f+j)->V[0]; vi[1]=(f+j)->V[1]; vi[2]=(f+j)->V[2];
     for(l=0;l<3;l++){
       R_m4by1(Lights[id].shadow_transform,
               (v+vi[l])->p[0],(v+vi[l])->p[1],(v+vi[l])->p[2],
               &x,&y,&z);
       if(y > 0.0){  /* vertex in front of light */
         xx[l]=(double)shadow_buffer_size2 + d*x/y;
         yy[l]=(double)shadow_buffer_size2 - d*z/y;
         pp[l][0]=x; pp[l][1]=y; pp[l][2]=z;
       }
       else goto FACEBREAK;
     }
     xxmin=min(xx[1],xx[2]);     xxmin=min(xx[0],xxmin);
     xxmax=max(xx[1],xx[2]);     xxmax=max(xx[0],xxmax);
     yymin=min(yy[1],yy[2]);     yymin=min(yy[0],yymin);
     yymax=max(yy[1],yy[2]);     yymax=max(yy[0],yymax);
     if((yymax > yymin) && (xxmax > xxmin)){
       VECSUB(pp[2],pp[0],pp[2])
       VECSUB(pp[1],pp[0],pp[1])
       cross(pp[1],pp[2],en);
       if(normalize(en) == 0){
         ec=DOT(pp[0],en);
         for(l=max(0,(long)yymin);l<=min((long)yymax,shadow_buffer_size-1);l++){
           if(ShadowMirrorLine(l,xx,yy,&xxmin,&xxmax,0.0,
                               (double)shadow_buffer_size)){
             d2=(double)(shadow_buffer_size2-l);
             UpdateShadowZbuffer(Lights[id].shadow_buffer,depth2,
                                 l,xxmin,xxmax,d,d2,en,ec,
                                 (double)shadow_buffer_size2);
           }
         }
       }
     }
     FACEBREAK:;
   }
 }
 for(i=0;i< shadow_buffer_size*shadow_buffer_size; i++){
   d2f = *(Lights[id].shadow_buffer+i);
   if(*(depth2+i) < 1.e25 && d2f < 1.e25){
     *(Lights[id].shadow_buffer+i) = ( d2f + *(depth2+i) )*0.5;
   }
   else{
     *(Lights[id].shadow_buffer+i) += 1.e30;
   }
   *(depth2+i) = d2f;
 }
// X__Free(depth2);
 Lights[id].atmospheric_shadow_buffer=depth2;
 return;
}

#define SBS shadow_buffer_size
//#define TOLL1 1.0
//#define TOLL2 200.0       // OK for atmospherics for atmospherics

#define TOLL1 10.0
#define TOLL2 0.0

static void UpdateShadowZbuffer(float *id, float *iq,
                                long l, double xl, double xr,
                                double d, double d2, vector en, double ec,
                                double SBC){
 long i,ixl,ixr;
 double depth,bL,bR,aL,aR,dL,dR;
 xl += 0.01;
 xr -= 0.01;
 ixl = (long)xl;
 ixr = (long)xr;
 if(ixl == ixr){
   dL = xl - SBC;
   dR = xr - SBC;
   bL = en[0]*dL + en[2]*d2 + en[1]*d ;
   if(fabs(bL) < 1.e-30)aL=BIGP;
   else aL = ec/bL;
   bR = en[0]*dR + en[2]*d2 + en[1]*d ;
   if(fabs(bR) < 1.e-30)aR=BIGP;
   else aR = ec/bR;
   depth=max(aL,aR)*d;
   if(fabs(depth - *(id+l*SBS+ixl)) < TOLL1){;}
   else if(depth > 0 && depth <= *(id+l*SBS+ixl)){
     if(*(id+l*SBS+ixl) <= *(iq+l*SBS+ixl)) *(iq+l*SBS+ixl) = *(id+l*SBS+ixl);
     *(id+l*SBS+ixl) = depth;
   }
   else if(depth > 0 && depth < *(iq+l*SBS+ixl)) *(iq+l*SBS+ixl) = depth;
 }
 else{
   dL = xl - SBC;
   bL = en[0]*dL + en[2]*d2 + en[1]*d ;
   if(fabs(bL) < 1.e-30)aL=BIGP;
   else aL = ec/bL; /* mu for vector {dL,d,d2} */
   bR = en[0]*(double)(ixl+1-SBC) + en[2]*d2 + en[1]*d ;
   if(fabs(bR) < 1.e-30)aR=BIGP;
   else aR = ec/bR; /* mu for vector {dL,d,d2} */
   depth=max(aL,aR)*d;
   if(fabs(depth - *(id+l*SBS+ixl)) < TOLL1){;}
   else if(depth > 0 && depth <= *(id+l*SBS+ixl)){
     if(*(id+l*SBS+ixl) <= *(iq+l*SBS+ixl)) *(iq+l*SBS+ixl) = *(id+l*SBS+ixl);
     *(id+l*SBS+ixl) = depth;
   }
   else if(depth > 0 && depth < *(iq+l*SBS+ixl)) *(iq+l*SBS+ixl) = depth;
   bL = en[0]*(double)(ixr-SBC) + en[2]*d2 + en[1]*d ;
   if(fabs(bL) < 1.e-30)aL=BIGP;
   else aL = ec/bL; /* mu for vector {dL,d,d2} */
   dR = xr - SBC;
   bR = en[0]*dR + en[2]*d2 + en[1]*d ;
   if(fabs(bR) < 1.e-30)aR=BIGP;
   else aR = ec/bR; /* mu for vector {dL,d,d2} */
   depth=max(aR,aL)*d;
   if(fabs(depth - *(id+l*SBS+ixr)) < TOLL1){;}
   else if(depth > 0 && depth <= *(id+l*SBS+ixr)){
     if(*(id+l*SBS+ixr) <= *(iq+l*SBS+ixr)) *(iq+l*SBS+ixr) = *(id+l*SBS+ixr);
     *(id+l*SBS+ixr) = depth;
   }
   else if(depth > 0 && depth < *(iq+l*SBS+ixr)) *(iq+l*SBS+ixr) = depth;
   if(ixr-ixl > 1)for(i=ixl+1;i<ixr;i++){
     bL = en[0]*(double)(i-SBC) + en[2]*d2 + en[1]*d ;
     if(fabs(bL) < 1.e-30)aL=BIGP;
     else aL = ec/bL;  /* mu for vector {dL,d,d2} */
     bR = en[0]*(double)(i+1-SBC) + en[2]*d2 + en[1]*d ;
     if(fabs(bR) < 1.e-30)aR=BIGP;
     else aR = ec/bR;  /* mu for vector {dL,d,d2} */
     depth=max(aL,aR)*d;
     if(fabs(depth - *(id+l*SBS+i)) < TOLL1){;}
     else if(depth > 0 && depth <= *(id+l*SBS+i)){
       if(*(id+l*SBS+i) <= *(iq+l*SBS+i)) *(iq+l*SBS+i) = *(id+l*SBS+i);
       *(id+l*SBS+i) = depth;
     }
     else if(depth > 0 && depth < *(iq+l*SBS+i)) *(iq+l*SBS+i) = depth;
   }
 }
}

#ifdef _SUNSTYLE
void PointInShadow(p,id,intensity)
vector p; long id; double *intensity; {
#else
void PointInShadow(vector p, long id, double *intensity){
#endif
  double x,y,z,zz,dx,dy,i0,ix,iy,ii;
  long xx,yy,i,j,
//  edge=0,
  edge=2,
  edge3,sum; /* double xx,yy; */
  R_m4by1(Lights[id].shadow_transform,p[0],p[1],p[2],&x,&y,&z);
  if(y > 0.0){  /* point in front of light */
    i0 = ix = iy = *intensity;
    dx=((double)shadow_buffer_size2 + Lights[id].scale*x/y);
    dy=((double)shadow_buffer_size2 - Lights[id].scale*z/y);
    xx=(long)dx; yy=(long)dy; sum=0; edge3=0;
    for(i=xx-edge;i<=xx+edge;i++)for(j=yy-edge;j<=yy+edge;j++){
      if(i >= 0 && i < SBS && j >= 0 && j < SBS){
        zz = *(Lights[id].shadow_buffer+j*SBS+i);
        if(zz < y)sum++; edge3++;
      }
    }
    if(edge3 > 0)i0 = i0 *
      (1.0 - (double)sum/(double)edge3*(1.0-shadow_density));
    xx++; sum=0; edge3=0;
    for(i=xx-edge;i<=xx+edge;i++)for(j=yy-edge;j<=yy+edge;j++){
      if(i >= 0 && i < SBS && j >= 0 && j < SBS){
        zz = *(Lights[id].shadow_buffer+j*SBS+i);
        if(zz < y)sum++; edge3++;
      }
    }
    if(edge3 > 0)ix = ix *
      (1.0 - (double)sum/(double)edge3*(1.0-shadow_density));
    yy++; xx--; sum=0; edge3=0;
    for(i=xx-edge;i<=xx+edge;i++)for(j=yy-edge;j<=yy+edge;j++){
      if(i >= 0 && i < SBS && j >= 0 && j < SBS){
        zz = *(Lights[id].shadow_buffer+j*SBS+i);
        if(zz < y)sum++; edge3++;
      }
    }
    if(edge3 > 0)iy = iy *
      (1.0 - (double)sum/(double)edge3*(1.0-shadow_density));
    yy--;
    dx -= (double)xx;
    dy -= (double)yy;
    ii = i0 + (ix - i0)*dx + (iy - i0)*dy;
    ii = min(ii,*intensity);
    *intensity = max(0.0,ii);
  }
  return;
}

#ifdef TRACER

/* this function is nessary only if full trace of shadows etc. */

#define TOL    1.e-6
#define TOLE   1.e-4

int R_intersect_face_old(vector p, vector d,
                     vector pf0, vector pf1, vector pf2, vector n,
                     double *depth, double tollerance){
/* tollerance is to allow for no hit on origin face */
/* p is origin of ray         d is (normalise ray direction */
/* other params are for face being tested */
 vector v1,pi;
 double mu,mu1,det1,det2,det3,ve1,ve2,ve3,vp1,vp2,vp3,p1,p2,p3,a,b,dm,detmax;
 int k;
 mu1=DOT(d,n);
 if(fabs(mu1) < TOL)return 0;   /* parallel to face */
 VECSUB(pf0,p,v1)
 mu=DOT(n,v1)/mu1;
 if(mu < tollerance)return 0;          /* behind the ray   */
 VECSCALE(mu,d,v1)
 VECSUM(p,v1,pi)
 ve1=pf1[0]-pf0[0]; ve2=pf1[1]-pf0[1]; ve3=pf1[2]-pf0[2];
 vp1=pf2[0]-pf0[0]; vp2=pf2[1]-pf0[1]; vp3=pf2[2]-pf0[2];
 p1=pi[0]-pf0[0];
 p2=pi[1]-pf0[1];
 p3=pi[2]-pf0[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 == 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;
 }
 if( a < -TOLE || a > 1.0+TOLE || b < -TOLE || b > 1.0+TOLE || (a+b) > 1.0+TOLE)return 0;
 *depth=mu;
 return 1;
}

#define TOLF 0.0001
#define TOLP 0.025   // 0.01 = 89.4 degrees.
#define TOLQ 100.0     // much smaller than face
int R_intersect_face(vector p, vector d,
                     vector pf0, vector pf1, vector pf2, vector n,
                     double *depth, double tollerance){
/* tollerance is to allow for no hit on origin face */
/* p is origin of ray         d is (normalise ray direction */
/* other params are for face being tested */
 vector v1,pi,w,u,v;
 double mu,a,b,de,uu,vv,wu,wv,uv; 
 mu=DOT(d,n);
 if(fabs(mu) < TOLP)return 0;   /* parallel to face */
 VECSUB(pf0,p,v1)
 mu=DOT(n,v1)/mu;
 if(mu < tollerance)return 0;   /* behind the ray origin  */
 VECSCALE(mu,d,v1)
 VECSUM(p,v1,pi)
 VECSUB(pi,pf0,w)
 VECSUB(pf1,pf0,u)
 VECSUB(pf2,pf0,v)
 uv=DOT(u,v);
 uu=DOT(u,u);
 vv=DOT(v,v);
 de=uu*vv-uv*uv;
 if(fabs(de) < TOLQ)return 0; 
 wu=DOT(w,u);
 wv=DOT(w,v);
 a=(wu*vv-wv*uv)/de; 
 b=(wv*uu-wu*uv)/de;
 if( a < -TOLF || a > 1.0+TOLF || b < -TOLF || b > 1.0+TOLF || (a+b) > 1.0+TOLF)return 0;
 *depth=mu;
 return 1;
}




#endif