/*
  *     program : ray8.c
  *
  *   OLD INFO:
  *  [compile with "cc -32 ray8.c -o r8 -L/usr/lib -lm -lc -limage -lgl"
  *
  * description : displays and rotates a wireframe model to choose viewpoint
  *               for rendering by ray-tracing, implementing diffuse lighting
  *               by many lights, correct shadows, multiple reflection,
  *               transparency and refraction (including total internal
  *               reflection), texturing (+ mapping of multiple .rgb files),
  *		  bump mapping. Resolution of image increases by factors of two,
  *               then averages 4,16,64,etc. rays per pixel.
  *               Algorithm is recursive.
  *		  Allows easy inclusion of new object types, ie it's generalised.
  *
  *      author : A.H.Barnett  NDTAC, Nuclear Electric plc, Wythenshawe]
  *
  *        date : 28/7/1993
  *        recompiled : 17/6/1996 at Harvard, Heller Group.
  *        6/9/2008: extracted just the computational core, writing raw floats
  *          compile under fedora linux: gcc r8_core.c -o r8 -lm
  */

 /* header files for input/output, maths and Graphics Library commands... */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <png.h>              /* for image output */ 

#define Boolean int
#define FALSE 0
#define TRUE 1

 /* maximum array dimensions: */

 /* objects... */
#define MAXO 256
 /* verteces... */
#define MAXV 128
 /* facets... */
#define MAXF 128
 /* verteces per facet... */
#define MAXN  16
 /* spheres and cylinders... */
#define MAXS 128
 /* lights... */
#define MAXL 16
 /* random sequence length... */
#define NRND 1000
 /* image buffer side... */
#define NIMG 400
 /* number of images... */
#define MAXI 3
 /* 2^31 - 1 */
#define MAXINT 2147483647.0

 /* object type descriptors (...what a flashy word !): */
#define SPHERE    0
#define FACET     1
 /* not implemented yet... */
#define CYLINDER  2

 /* diffuse colour texture type descriptors: */
#define UNTEXTURED   0
#define WOOD         1
#define XWOOD        2
#define WOBBLE       3
#define ZCHEQUERS    4
#define GRANITE      5
#define ZTILES       6
#define XBRICKS      7
#define YBRICKS      8
#define MOTTLED      9
#define MARBLE       10
 /* at 50 and above, units give image number, tens give mapping type: */
#define ZIMAGE       50
#define ZSPHEREWRAP  60
#define ZIMGGOURAUD  70
#define ZSPHGOURAUD  80

 /* bump mapping texture descriptors: */
#define FROSTED1     1
#define FROSTED2     2
#define FROSTED3     3
#define BUMPY1       10
#define BUMPY2       11
#define ZPOOL1       20
#define ZPOOL2       21


 /* setting up arrays: */

 /* total number of objects... */
short nobj;
 /* object pointer arrays: object i has type object[i], and is the
  * pointer[i]'th in that type's coordinate arrays.
  * texture[i][0] gives diffuse colour texturing for object i,
  * texture[i][1] gives bump mapping texture, or untextured if zero... */
short object[MAXO];
short pointer[MAXO];
short texture[MAXO][2];

 /* title of scene... */
char title[100];


 /* OBJECT PROPERTIES ARRAYS */

 /* diffuse colour intensities of objects... */
float objrgb[MAXS][3];
 /* reflectivity... */
float objref[MAXS];
 /* transparency of surfaces... */
float objtra[MAXS];
 /* refractive index, must be >= 1.0 with current refraction algorithm... */
float objrin[MAXS];


 /* COORDINATE ARRAYS */

 /* number of spheres... */
short ns;
 /* XYZ coords of sphere centres... */
float cen[MAXS][3];
 /* radii... */
float rad[MAXS];

 /* number of cylinders... */
short ncyl;
 /* XYZ coords of centre of each end... */
float cyl[MAXS][6];
 /* radii... */
float cylrad[MAXS];
 /* normals, end plane equations... */
float cnor[MAXS][3],cd[MAXS][2];

 /* number of vertices... */
short nver;
 /* XYZ coords of verteces... */
float ver[MAXV][3];

 /* number of facets... */
short nfac;
 /* number of verteces joined by each facet... */
short nvert[MAXF];
 /* the vertex list for each facet (verteces numbered from 0)... */
short fa[MAXF][MAXN];
 /* facet equations... NOTE: vectors from 0-1 and 1-2 in vertex list mustn't be parallel! */
float fnor[MAXF][3], fd[MAXF];
 /* facet boundary equations... */
float fbnor[MAXF][MAXN][3], fbd[MAXF][MAXN];

 /* number of lights... */
short nlig;
 /* light-vectors (point towards light source), do not need to be unit length... */
float lv[MAXL][3];
 /* light colour intensities RGB...  */
float lrgb[MAXL][3];

 /* ambient, sky colour RGB... */
float ambrgb[3],skyrgb[3];


 /* initial viewing position... */
float ex = 0.0, ey = -15.0, ez = 0.0;
float dist = 15.0;
float th = 0.4, ph = -0.7, rh = 0, xof = 0, yof = 0;

 /* 2 * tangent of half-width of view in Y-direction... */
float iris = 0.6;

 /* maximum number of reflections for recursive ray-tracing algorithm... */
short maxdep = 5;   /* 2 */

 /* reflections menu list (must match text in menu string)... */
short deplist[9] = {0,1,2,3,4,5,10,20};

 /* on/off flags for controlling: */
 /* chrome-mapped sky... */
Boolean sky = TRUE;
 /* enhanced grazing reflectivity... */
Boolean graze = TRUE;
 /* allows diffuse texturing, bump mapping, transparency, depth-cueing... */
Boolean txt = TRUE, bmp = TRUE, trn = TRUE, dcue = FALSE;
 /* displaying facet's vertex normals... */
Boolean vertnorms = TRUE;
 /* displaying numbers... */
Boolean vertnums = TRUE;
 /* axes... */
Boolean axes = TRUE;

 /* colours for certain display items... */
#define TITLECOL 0x007FFFFF
#define  TEXTCOL 0x00FF7F00
#define   CENCOL 0x0000FFFF
#define   VERCOL 0x0000FF00
#define  WIRECOL 0x00007FFF
#define   FACCOL 0x00BFFF00
#define  NORMCOL 0x0080BF00
#define  AXESCOL 0x003F00FF

 /* XYZ coords for axes... */
float orig[3] = { 0.0,0.0,0.0 };
float axcoord[3][3] = { {1.0,0.0,0.0}, {0.0,1.0,0.0}, {0.0,0.0,1.0} };

 /* wire-frame sphere approximation: */
 /* number of points per circle... */
#define NSPH 24

 /* array for 3 circles... */
float sphmodel[3][NSPH][3];

#define PI 3.1415926536

 /* image buffers, info... */
short imgr[MAXI][NIMG][NIMG],imgg[MAXI][NIMG][NIMG],imgb[MAXI][NIMG][NIMG];
short imgx[MAXI],imgy[MAXI];
float imgsize[MAXI],imgasp[MAXI];
char iname[40];
short inamecount;

 /* help window information... */
#define HELPLEN 22
char *helptext[HELPLEN] = {

 "    Ray-Tracing Package",
 "    Version 8",
 "",
 "    A.H.Barnett 28/7/1993",
 "",
 " Diffuse multiple lighting,",
 "  ambient light and sky colour,",
 "  chrome environment mapping,",
 " Correct shadow generation, with",
 "  semi-opaque shadows,",
 " Multiple specular reflections",
 "  (with angle-dependence option)",
 "  including off inside surfaces,",
 " Transparency, with refraction",
 "  and total internal reflection,",
 " Texturing of diffuse colours,",
 "  including multiple .rgb image",
 "  mapping onto planes, spheres,",
 " Bump mapping of normals,",
 " Model read and write to file,",
 " General object representation,",
 "  currently spheres and facets."  };

 /* other global variables: */
 /* window size and old, new mouse positions... */
long xmax, ymax, omx, mx, omy, my;

 /* menu identifiers... */
long menu,toggmen,reflmen,filemen,imagmen;

 /* window identifiers... */
long win,helpwin;

 /* array for filenames... */
char *fname;

 /* some keys for input function: */
#define KEYR 13
#define KEYD 127
#define KEYBACK 8

short inp();

 /* rnd function seed, arrays for texture... */
unsigned long seed, seq[NRND+1][3];

unsigned long rnd();
long nois();
float mottle();
 /* add these three to get to compile 17/6/96 */
void texturing();
void gouraudimage();
void bumpmap();
void render(float *img, int samples);

 /* flag for main program loop exiting... */
Boolean progexit = FALSE;



 /* main program... MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM */
main(argc, argv)
int argc;
char *argv[];
{
  int pix, nimg, samples;
  float *img;
  char *oname;
  FILE *fp;

  fname = malloc(40);          /* input name */
  oname = malloc(40);           /* output name */ 

  if (argc!=12) {
    fprintf(stderr, "usage: r8 nx ny nsamples infile outfile ex ey ez th ph iris\n\neg: r8 1000 1000 2 bob.r bob.img 2 -10 5 0.2 -.4 .6\nr8 1000 1000 2 pool.r pool.img 2 -12 5 0.3 -.5 .8\n");
    return;
  }    
  sscanf(argv[1], "%ld", &xmax);
  sscanf(argv[2], "%ld", &ymax);
  sscanf(argv[3], "%d", &samples);
  sscanf(argv[4], "%s", fname);
  sscanf(argv[5], "%s", oname);
  sscanf(argv[6], "%f", &ex);
  sscanf(argv[7], "%f", &ey);
  sscanf(argv[8], "%f", &ez);
  sscanf(argv[9], "%f", &th);
  sscanf(argv[10], "%f", &ph);
  sscanf(argv[11], "%f", &iris);

  printf("%ld,%ld (%d): (%f,%f,%f) (%f,%f) %f\n", xmax, ymax, samples,
	 ex, ey, ez, th, ph, iris);
  
  nimg = 3*xmax*ymax;
  img = malloc(sizeof(float)*nimg);   /* RGB, RGB, etc */

 /* set up graphics, window, menus, etc... */
  initialize();
  
  /* choose viewing position ? ... */

  render(img, samples);        /* compute */

  printf("writing %s ...\n", oname);
  fp = fopen(oname, "w");
  fwrite(&xmax,4,1,fp);
  fwrite(&ymax,4,1,fp);
  fwrite(img,sizeof(float), nimg, fp);
  close(fp);

  free(img);
} /* MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM */



void calcsetup()

 /* calculate facet equations, facet boundary equations, normalise light-vectors. */
{
    short i,j,nextj,k;
    float x,y,z,l,vo[3],vi[3];

 /* i loops round all facets... */
    for (i=0; i<nfac; ++i) {
     /* use first two facet boundary vectors, find X product and normalise, to give facet normals... */
        for (k=0; k<3; ++k) {
            vo[k] = ver[fa[i][1]][k] - ver[fa[i][2]][k];
            vi[k] = ver[fa[i][1]][k] - ver[fa[i][0]][k];
        }
        x = vo[1]*vi[2] - vo[2]*vi[1];
        y = vo[2]*vi[0] - vo[0]*vi[2];
        z = vo[0]*vi[1] - vo[1]*vi[0];
        l = sqrt(x*x + y*y + z*z);
        fnor[i][0] = (x /= l);
        fnor[i][1] = (y /= l);
        fnor[i][2] = (z /= l);
     /* facet d = normal.(position of vertex 0) ... */
        fd[i] = ( x*ver[k = fa[i][0]][0] + y*ver[k][1] + z*ver[k][2] );

     /* j loops round facet's verteces... */
	for (j=0; j<nvert[i]; ++j) {
	 /* wrap nextj around, use modulus nvert[i]... */
	    nextj = (j+1) % nvert[i];
	 /* same as above, but product is between one boundary vector and the normal... */
            for (k=0; k<3; ++k) {
		vo[k] = ver[fa[i][j]][k] - ver[fa[i][nextj]][k];
		vi[k] = fnor[i][k];
	    }
	    x = vo[1]*vi[2] - vo[2]*vi[1];
	    y = vo[2]*vi[0] - vo[0]*vi[2];
	    z = vo[0]*vi[1] - vo[1]*vi[0];
	    l = sqrt(x*x + y*y + z*z);
	    fbnor[i][j][0] = (x /= l);
	    fbnor[i][j][1] = (y /= l);
	    fbnor[i][j][2] = (z /= l);
	 /* facet d = normal.(position of vertex fa[i][j]) ... */
	    fbd[i][j] = ( x*ver[k = fa[i][j]][0] + y*ver[k][1] + z*ver[k][2] );
	}
    }

 /* i loops around all cylinders... */
    for (i=0; i<ncyl; ++i) {
	x = cyl[i][3]-cyl[i][0];
	y = cyl[i][4]-cyl[i][1];
	z = cyl[i][5]-cyl[i][2];
        l = sqrt(x*x + y*y + z*z);
        cnor[i][0] = (x /= l);
        cnor[i][1] = (y /= l);
        cnor[i][2] = (z /= l);
     /* end planes: d = normal.(position of end) ... */
	cd[i][0] = x*cyl[i][0] + y*cyl[i][1] + z*cyl[i][2];
	cd[i][1] = x*cyl[i][3] + y*cyl[i][4] + z*cyl[i][5];
    }

    for (i=0; i<nlig; ++i) {
        l = sqrt(lv[i][0]*lv[i][0] + lv[i][1]*lv[i][1] + lv[i][2]*lv[i][2]);
        lv[i][0] /= l;
        lv[i][1] /= l;
        lv[i][2] /= l;
    }

}



initialize()

 /* sets up graphics, window, menus, devices to be queued, sphere approximation
  * coords, default filename, facet normals, and normalises light-vectors, loads default file */
{

    short i,ch;
    float s,c;
    FILE *fptr;

 /* open, load, close file... */
    if ((fptr = fopen(fname,"r")) == NULL) printf("file: %s, open error...\n",fname);
    load(fptr);
    fclose(fptr);

 /* set up sphere approximation coords... */
    for (i=0; i<NSPH; ++i) {
        s = sin(i*2.0*PI/NSPH);
        c = cos(i*2.0*PI/NSPH);
        sphmodel[0][i][0] = 0;
        sphmodel[0][i][1] = c;
        sphmodel[0][i][2] = s;
        sphmodel[1][i][0] = s;
        sphmodel[1][i][1] = 0;
        sphmodel[1][i][2] = c;
        sphmodel[2][i][0] = c;
        sphmodel[2][i][1] = s;
        sphmodel[2][i][2] = 0;
    }

 /* start-up rnd function... */
    seed = 174924501;
    for (i=0; i<10; ++i)
        rnd(seed);

 /* make random sequence arrays, with non-zero values, and wrap-around... */
    for (i=0; i<NRND; ++i) {
	while ((seq[i][0] = rnd(seed)) == 0);
	while ((seq[i][1] = rnd(seed)) == 0);
	while ((seq[i][2] = rnd(seed)) == 0);
    }
    seq[NRND][0] = seq[0][0];
    seq[NRND][1] = seq[0][1];
    seq[NRND][2] = seq[0][2];

 /* initialise image parameters... */
    for (i=0; i<MAXI; ++i) {
	imgx[i] = imgy[i] = 100;
	imgasp[i] = imgsize[i] = 1.0;
    }

    calcsetup();

}

load(ptr)
FILE *ptr;
 /* loads a model from already-opened file (pointer ptr), into program's arrays. */
{
    short i,j,c;
    float x,y,z,r,s,t;
    char *deftitle = "Untitled";

 /* defaults set before loading... */
    ns = nver = nfac = ncyl = nobj = nlig = 0;
 /* copy default title into title array... */
    for (i=0; (c = deftitle[i]); title[i++] = c);

 /* loop; interpret type character and read file accordingly... */
    while ((c = getc(ptr)) != 'E') {
        switch (c) {
            case '#':
             /* enclose comments between #s: they are ignored... */
                while (getc(ptr) != '#') ;
                break;
            case 'T':
             /* read title... */
                getc(ptr);
                for (i=0; (c = getc(ptr)) != '\n'; title[i++] = c) ;
                title[i] = 0;
                break;
	    case 'K':
	     /* read sky colour... */
                fscanf(ptr," %f %f %f",&x,&y,&z);
		skyrgb[0] = x;
		skyrgb[1] = y;
		skyrgb[2] = z;
		break;
	    case 'A':
	     /* read ambient colour... */
                fscanf(ptr," %f %f %f",&x,&y,&z);
		ambrgb[0] = x;
		ambrgb[1] = y;
		ambrgb[2] = z;
		break;
            case 'S':
		object[nobj] = SPHERE;
		pointer[nobj] = ns;
             /* read sphere data... */
                fscanf(ptr," %f %f %f %f",&x,&y,&z,&r);
                cen[ns][0] = x;
                cen[ns][1] = y;
                cen[ns][2] = z;
                rad[ns++] = r;
	     /* properties... */
                fscanf(ptr," %f %f %f %hd %hd %f %f %f",&x,&y,&z,&i,&j,&r,&s,&t);
                objrgb[nobj][0] = x;
                objrgb[nobj][1] = y;
                objrgb[nobj][2] = z;
		texture[nobj][0] = i;
		texture[nobj][1] = j;
                objref[nobj] = r;
                objtra[nobj] = s;
                objrin[nobj++] = t;
                break;
            case 'F':
		object[nobj] = FACET;
		pointer[nobj] = nfac;
             /* read facet vertex data until -1 found... */
		i = 0;
		fscanf(ptr," %hd",&c);
		while (c != -1 && i < MAXN) {
		    fa[nfac][i++] = c;
		    fscanf(ptr," %hd",&c);
		}
		nvert[nfac++] = i;
	     /* properties... */
                fscanf(ptr," %f %f %f %hd %hd %f %f %f",&x,&y,&z,&i,&j,&r,&s,&t);
                objrgb[nobj][0] = x;
                objrgb[nobj][1] = y;
                objrgb[nobj][2] = z;
		texture[nobj][0] = i;
		texture[nobj][1] = j;
                objref[nobj] = r;
                objtra[nobj] = s;
                objrin[nobj++] = t;
                break;
	    case 'V':
	     /* read vertex coord data... */
		fscanf(ptr," %hd",&i);
		if (i >= nver) nver = i+1;
		fscanf(ptr," %f %f %f",&x,&y,&z);
		ver[i][0] = x;
		ver[i][1] = y;
		ver[i][2] = z;
		break;
	    case 'C':
	     /* cylinder coord data... */
		object[nobj] = CYLINDER;
		pointer[nobj] = ncyl;
             /* read sphere data... */
                fscanf(ptr," %f %f %f %f %f %f",&x,&y,&z,&r,&s,&t);
                cyl[ncyl][0] = x;
                cyl[ncyl][1] = y;
                cyl[ncyl][2] = z;
		cyl[ncyl][3] = r;
		cyl[ncyl][4] = s;
		cyl[ncyl][5] = t;
                fscanf(ptr," %lf",cylrad+ncyl);
		++ncyl;
	     /* properties... */
                fscanf(ptr," %f %f %f %hd %hd %f %f %f",&x,&y,&z,&i,&j,&r,&s,&t);
                objrgb[nobj][0] = x;
                objrgb[nobj][1] = y;
                objrgb[nobj][2] = z;
		texture[nobj][0] = i;
		texture[nobj][1] = j;
                objref[nobj] = r;
                objtra[nobj] = s;
                objrin[nobj++] = t;
                break;
            case 'L':
             /* read lighting data... */
                fscanf(ptr," %f %f %f %f %f %f",&x,&y,&z,&r,&s,&t);
                lv[nlig][0] = x;
                lv[nlig][1] = y;
                lv[nlig][2] = z;
                lrgb[nlig][0] = r;
                lrgb[nlig][1] = s;
                lrgb[nlig++][2] = t;
                break;
            default:
                break;
        }
    }

    calcsetup();
}

void render(float *img, int samples)
/* samples is number of sub-pixel samples in each direction (eg 1,2,4 etc) 
 * gutted of all but caluclational content, barnett 6/9/08
*/
{
  short x,y,i,j,k,pix, val;
  long nxmax,nymax,dev;
  double qx,qy,qz,ox,oy,oz,st,ct,sp,cp,sr,cr;
  float prgb[3],rgb[3];

 /* set up sines, cosines for viewing transformation... */
    st = sin(th); ct = cos(th);
    sp = sin(ph); cp = cos(ph);
    sr = sin(rh); cr = cos(rh);

    pix = 1;
    
    for (x=0; x<xmax; x+=pix) {
      for (y=0; y<ymax; y+=pix) {
	prgb[0] = prgb[1] = prgb[2] = 0;
	for (i=0; i<samples; ++i)
	  for (j=0; j<samples; ++j) {
	    qx = ((x + i*pix/(float)samples)/xmax - 0.5)*iris*xmax/(float)ymax;
	    qy = 1.0;
	    qz = ((y + j*pix/(float)samples)/ymax - 0.5)*iris;
	    ox =  qx*cr + qz*sr;
	    oz = -qx*sr + qz*cr;
	    qz =  oz*cp + qy*sp;
	    oy = -oz*sp + qy*cp;
	    qy =  oy*ct + ox*st;
	    qx = -oy*st + ox*ct;
	    ox = sqrt(qx*qx + qy*qy + qz*qz);
	    qx /= ox;
	    qy /= ox;
	    qz /= ox;
	    /* printf("ray:\n"); */
	    raytrace(ex,ey,ez,qx,qy,qz,1,rgb);
	    for (k=0; k<3; ++k) prgb[k] += rgb[k];
	  }
	/* write into img array... */
	for (k=0; k<3; ++k) img[y+x*ymax+k*(xmax*ymax)] = \
			      (float) prgb[k] / (samples*samples);
			      /* (.5*k*x)/xmax;           was a test */
      }
    }


}



short intersect(px,py,pz,qx,qy,qz,laminptr)
double px,py,pz,qx,qy,qz;
double *laminptr;
 /* Returns j if line r = p + la.q intersects object j at a distance (ila) less than
  * lamin, in which case it sets lamin to be ila. Returns -1 if no intersections.
  * Needs algorithms for each object type. */
{
    short i,j,k,obj = -1;
    double la,fla,ila,ql,cx,cy,cz,dx,dy,dz,dl,rr,tiny = 0.0001;
    Boolean outoffacet;

  for (j=0; j<nobj; ++j) {
    i = pointer[j];
    switch (object[j]) {

    case SPHERE:
	la = qx*(cx=(cen[i][0]-px)) + qy*(cy=(cen[i][1]-py)) + qz*(cz=(cen[i][2]-pz));
	if ((dx = -cx + la*qx)<rad[i] && (dy = -cy + la*qy)<rad[i] && (dz = -cz + la*qz)<rad[i])
	    if ((dl = dx*dx + dy*dy + dz*dz) < (rr = rad[i]*rad[i])) {
		fla = sqrt(rr - dl);
		if ((ila = la-fla) > tiny && ila < *laminptr) {
		    *laminptr = ila;
		    obj = j;
		} else if ((ila = la+fla) > tiny && ila < *laminptr) {
		    *laminptr = ila;
		    obj = j;
		}
	    }
	break;

    case CYLINDER:
	break;

    case FACET:
	la = (fd[i] - (fnor[i][0]*px + fnor[i][1]*py + fnor[i][2]*pz)) / (fnor[i][0]*qx + fnor[i][1]*qy + fnor[i][2]*qz);
	if (la > tiny && la < *laminptr) {
	 /* c is point of intersection... */
	    cx = px + la*qx;
	    cy = py + la*qy;
	    cz = pz + la*qz;
	    outoffacet = FALSE;
	 /* check against all boundary planes to see if inside facet... */
	    for (k=0; k<nvert[i]; ++k) {
		if ((fbnor[i][k][0]*cx + fbnor[i][k][1]*cy + fbnor[i][k][2]*cz) < (fbd[i][k]-tiny)) {
		    outoffacet = TRUE;
		 /* skip to end of loop... */
		    k = nvert[i];
		}
	    }
	    if (!outoffacet) {
		*laminptr = la;
		obj = j;
	    }
	}
	break;
    }
  }
  return obj;
}

float shadow(px,py,pz,qx,qy,qz,currobj)
short currobj;
float qx,qy,qz;
double px,py,pz;
 /* returns shadowing factor (between 0 and 1) at point p for a light-vector q.
  * currobj is the object on which shadow is needed...
  * Needs algorithms for each object type, efficient with opaque objects.
  * single-precision used for speed. */
{
    short i,j,k;
    float sh = 1.0;
    float la,fla,ql,cx,cy,cz,dx,dy,dz,dl,rr,tiny = 0.0001;
    Boolean outoffacet;

    for (j=0; j<nobj; ++j) {
	i = pointer[j];
	switch (object[j]) {

	case SPHERE:
	if (objtra[j] < tiny) {
	 /* efficient handling of opaque spheres: don't deal with current one... */
	    if (j!=currobj) {
		la = qx*(cx=(cen[i][0]-px)) + qy*(cy=(cen[i][1]-py)) + qz*(cz=(cen[i][2]-pz));
		if (la>0.0 && (dx = -cx + la*qx)<rad[i] && (dy = -cy + la*qy)<rad[i] && (dz = -cz + la*qz)<rad[i])
		 /* check for intersection... */
		    if ((dx*dx + dy*dy + dz*dz) < rad[i]*rad[i]) {
			sh = 0.0;
		     /* skip to end of objects... */
			j = nobj;
		    }
	    }
	} else {
	 /* partial shadows for transparent spheres... */
	    la = qx*(cx=(cen[i][0]-px)) + qy*(cy=(cen[i][1]-py)) + qz*(cz=(cen[i][2]-pz));
	 /* la>-rad[i] test ignores spheres which are completely behind light-vector... */
	    if (la>-rad[i] && (dx = -cx + la*qx)<rad[i] && (dy = -cy + la*qy)<rad[i] && (dz = -cz + la*qz)<rad[i])
		if ((dl=dx*dx + dy*dy + dz*dz) < (rr=rad[i]*rad[i])) {
                        fla = sqrt(rr - dl);
		     /* shadowed by front & back surfaces of sphere... */
                        if (la-fla > tiny) sh *= objtra[j]*objtra[j];
		     /* shadowed only by front surface of sphere... */
                        else if (la+fla > tiny) sh *= objtra[j];
                    }
	}
	break;

	case CYLINDER:
	break;

	case FACET:
	la = (fd[i] - (fnor[i][0]*px + fnor[i][1]*py + fnor[i][2]*pz)) / (fnor[i][0]*qx + fnor[i][1]*qy + fnor[i][2]*qz);
	if (la > tiny) {
	 /* c is point of intersection... */
	    cx = px + la*qx;
	    cy = py + la*qy;
	    cz = pz + la*qz;
	    outoffacet = FALSE;
	 /* check against all boundary planes to see if inside facet... */
	    for (k=0; k<nvert[i]; ++k) {
		if ((fbnor[i][k][0]*cx + fbnor[i][k][1]*cy + fbnor[i][k][2]*cz) < (fbd[i][k]-tiny)) {
		    outoffacet = TRUE;
		 /* skip to end of loop... */
		    k = nvert[i];
		}
	    }
	    if (!outoffacet)
		if (objtra[j] < tiny) {
		    sh = 0.0;
		 /* skip to end of objects... */
		    j = nobj;
		} else sh *= objtra[j];
	}
	break;

	}
    }
    return sh;
}

raytrace(px,py,pz,qx,qy,qz,dep,cp)
double px,py,pz,qx,qy,qz;
float *cp;
short dep;
 /* Recursive procedure which traces along line r = p + la.q, until
  * a recursive depth maxdep is reached, and returns the RGB colour
  * in cp (colour pointer to 3-element array). la is lambda, the
  * distance along the line from p.
  * Contributions to RGB come from:
  *			1) diffuse intensity at the point, shadowed correctly,
  *	 if reflective, 2) raytrace(reflected ray),
  * and if transparent, 3) raytrace(refracted ray).
  * Apart from surface normals, code is essentially object type independent. */
{
    short i,j,k;
    double la,fla,ila,lamin,ql,cx,cy,cz,dx,dy,dz,dl,rr,ncp,cr,ssr,ni;
    double nx,ny,nz,cang,depth,dummy,tiny = 0.0001;
    float rgb[3],textrgb[3],irgb[3];

/*  printf("dep: %d\n",dep); */
/*printf("qx,y,z: %f %f %f\n",qx,qy,qz); */

 /* initialise minimum lambda to 'a long way away'... */
    lamin = 1.0E+10;

 /* check intersections with all objects: j is closest object, or
  * -1 if none intersected... */
    j = intersect(px,py,pz,qx,qy,qz,&lamin);

 /* zero the colour carrier... */
    irgb[0] = irgb[1] = irgb[2] = 0.0;

    if (j==-1) {
     /* no objects: sky is visible... */
        for (i=0; i<nlig; ++i) {
            rr = qx*lv[i][0] + qy*lv[i][1] + qz*lv[i][2];
            if (rr > 0.9) {
                for (k=0; k<7; ++k) rr *= rr;
                for (k=0; k<3; ++k) irgb[k] += rr * lrgb[i][k];
            }
        }
        if (!sky) {
     /* use (1+cos(angle from vertical))/2 as sky intensity... */
	/* Idea: check if dep=1, then can have non-reflective background sky! */
            rr = (1.0 + qz)/2;
            for (k=0; k<3; ++k) irgb[k] += rr * skyrgb[k];
        } else if (dep>1) {
	 /* primitive chrome mapping... */
	    irgb[0] += mottle(qx*10,qy*10,qz*10,0);
	    irgb[1] += mottle(qx*9,qy*9,qz*9,1);
	    irgb[2] += mottle(qy*7,qz*11,qx*4,0);
	}
     /* no depth-cueing for light sources... */
	lamin =0.0;
    } else {
     /* object j visible... 
      * c is new point for starting ray, ie new p... */
        cx = px + lamin*qx;
        cy = py + lamin*qy;
        cz = pz + lamin*qz;
/* printf("dep = %d,j = %d,c = %f %f %f",dep,j,cx,cy,cz); */
     /* find object normal at c...
      * Needs algorithms for each object type. */
	i = pointer[j];
	switch (object[j]) {
	case SPHERE:
            nx = (cx - cen[i][0])/rad[i];
            ny = (cy - cen[i][1])/rad[i];
            nz = (cz - cen[i][2])/rad[i];
	    break;
	case CYLINDER:
	    break;
	case FACET:
	    nx = fnor[i][0];
	    ny = fnor[i][1];
	    nz = fnor[i][2];
	    break;
	}
     /* modulate normal if bump mapping activated... */
	if (bmp && texture[j][1] != 0)
	    bumpmap(cx,cy,cz,j,&nx,&ny,&nz);

     /* texturing colours into array textrgb... */
	if (txt && texture[j][0] != 0)
	    texturing(cx,cy,cz,j,textrgb);
	else for (k=0; k<3; ++k) textrgb[k] = objrgb[j][k];

     /* diffuse light, shadowing calculation: */
	for (i=0; i<nlig; ++i) {
	    ncp = nx*lv[i][0] + ny*lv[i][1] + nz*lv[i][2];
	 /* illuminate diffusely front & back surfaces of facets, front only for others... */
	    if (object[j] == FACET) ncp = ((ncp>0.0) ? ncp : -ncp);
	    else if (ncp<0.0) ncp = 0.0;
	    ncp = shadow(cx,cy,cz,lv[i][0],lv[i][1],lv[i][2],j) * ncp;
	    for (k=0; k<3; ++k) irgb[k] += ncp * textrgb[k] * lrgb[i][k];
	}
     /* ambient lighting... */
	for (k=0; k<3; ++k) irgb[k] += textrgb[k] * ambrgb[k];

     /* braching reflected & refracted ray calculation: */
        ncp = 0.0;
     /* only reflect if not at bottom of recursive depth... */
        if (dep <= maxdep)
            if (objref[j] > tiny) {
	     /* calculate normal component n.q ... */
                ncp = nx*qx + ny*qy + nz*qz;
                raytrace(cx,cy,cz,qx-2*ncp*nx,qy-2*ncp*ny,qz-2*ncp*nz,dep+1,rgb);
	     /* if graze is TRUE, make reflectivity tend to 1 as ncp tends to 0... */
		rr = graze ? (objref[j] + (1.0 - objref[j])*(1.0 - sqrt(sqrt((ncp>0.0)?ncp:-ncp)))) : objref[j];
                for (k=0; k<3; ++k) irgb[k] += rr * rgb[k];
            }
     /* always refract regardless of recursive depth... */
        if (trn && objtra[j] > tiny) {
            if (objrin[j] > (1.0+tiny)) {
	     /* we're dealing with refraction... make sure ncp is known */
                if (ncp == 0.0) ncp = nx*qx + ny*qy + nz*qz;
	     /* sign of ncp determines in/out-ness of ray: */
                if (ncp < 0.0) {
                 /* ray is entering object... */
                    cr = sqrt(1.0 - (1.0 - ncp*ncp) /objrin[j]/objrin[j]);
                    ni = ncp + objrin[j]*cr;
                 /* although simply dividing vector q - n*ni by objrin[j] gives the
                  * correct refracted ray, the length is not 1.0 with sufficient
                  * accuracy for high objrin[j], so dividing by dl is used: */
                    dl = sqrt((dx=qx-nx*ni)*dx + (dy=qy-ny*ni)*dy + (dz=qz-nz*ni)*dz);
		 /* don't increase recursion depth for transmission... */
                    raytrace(cx,cy,cz,dx/dl,dy/dl,dz/dl,dep,rgb);
                } else {
                 /* ray is exiting object... don't sqrt(1-ssr) yet- it might be -ve! */
                    ssr = (1.0 - ncp*ncp) * objrin[j]*objrin[j];
                    if (ssr < 1.0) {
                     /* refraction...  */
                        ni = ncp - sqrt(1.0 - ssr)/objrin[j];
                        dl = sqrt((dx=qx-nx*ni)*dx + (dy=qy-ny*ni)*dy + (dz=qz-nz*ni)*dz);
                        raytrace(cx,cy,cz,dx/dl,dy/dl,dz/dl,dep,rgb);
                     /* ... or, total internal reflection... only when exiting, though */
                    } else {
                        if (dep <= maxdep) raytrace(cx,cy,cz,qx-2*ncp*nx,qy-2*ncp*ny,qz-2*ncp*nz,dep+1,rgb);
                     /* when internal reflections bottom-out, substitute black... */
                        else rgb[0] = rgb[1] = rgb[2] = 0.0;
                    }
                }
            } else
	     /* no bending of ray: continue from where left off... */
		raytrace(cx,cy,cz,qx,qy,qz,dep,rgb);
	 /* add in contribution from transmitted ray... */
            for (k=0; k<3; ++k) irgb[k] += objtra[j] * rgb[k];
        }

    }

 /* depth-cueing... */
    depth = dcue ? exp(-lamin/50.0) : 1.0;
 /* pass back colour carrier... */
    *cp     = irgb[0] * depth;
    *(cp+1) = irgb[1] * depth;
    *(cp+2) = irgb[2] * depth;
}

void texturing(x,y,z,j,trgb)
short j;
double x,y,z;
float *trgb;
 /* Feeds rgb value into pointer trgb, calculated from xyz coords using
  * algorithm chosen by texture[j][0], the diffuse texture type... */
{
    short t,i,imageno;
    long ix,iy,iz;
    Boolean flag;
    float a,r,theta,modx,mody,modz,size,border,tiny=0.0001;
 /* wood grain wave-vector... */
    float kx=20.0, ky=12.0, kz=8.0;
 /* standard grey colour... */
    *trgb     = 0.6;
    *(trgb+1) = 0.6;
    *(trgb+2) = 0.6;

 /* decode which image is used in mapping, set t as texture type... */
    t = texture[j][0];
    if (t >= 50) {
	imageno = t % 10;
	t = t-imageno;
    }

 /* a is factor to modulate rgb... -ve a gives grey colour (mortar?) unless
  * *trgb are chosen by routine... */
    switch (t) {

    case WOOD:
	a = (2. + cos(kx*x + ky*y + kz*z))/3;
    break;
    case XWOOD:
	a = (3. + cos(-kx*y + ky*x - kz*z))/4;
    break;
    case WOBBLE:
	a = (2. + cos(kx*x   -ky*y + kz*(z+ cos(8.0*(z+y))) ))/3;
    break;
    case ZCHEQUERS:
	size = 1.0;
	modx = (x<0) ? (1.-x)/size : x/size;
	mody = (y<0) ? (1.-y)/size : y/size;
	a = 0.2 + 0.8*(((int)modx % 2 + (int)mody % 2) % 2);
    break;
    case GRANITE:
	size = 0.02;
     /* integer parts... highest int less than float, ok for -ve values... */
	ix = (x>-tiny) ? (long)(x/size +tiny) : (long)(x/size +tiny) - 1;
	iy = (y>-tiny) ? (long)(y/size +tiny) : (long)(y/size +tiny) - 1;
	iz = (z>-tiny) ? (long)(z/size +tiny) : (long)(z/size +tiny) - 1;
	a = 1.0 - 0.8*nois(ix,iy,iz)/MAXINT;
    break;
    case MOTTLED:
	size = 0.1;
	a = 1.0 - 0.8*mottle(x/size,y/size,z/size,0);
    break;
    case MARBLE:
	a = 0.0;
     /* generate fractal texturing using sum of mottles on many scales... */
	for (size=2.0; size>0.02; size/=2) a += mottle(x/size,y/size,z/size,0) - 0.5;
	a = 0.5 + 0.4*a;
     /* clip */
	if (a<0.0) a=0.0;
	else if (a>1.0) a=1.0;
    break;
    case ZTILES:
	size = 1.5;
	border = 0.1;
	modx = x/size;
	mody = y/size/2.;
	modx = (x>0) ? modx-(int)modx : 1.+modx-(int)modx;
	mody = (y>0) ? mody-(int)mody : 1.+mody-(int)mody;
	if (mody<0.5) flag = mody<border/2. || modx<border;
	else flag = mody<(1.+border)/2. || (modx>0.5 && modx<(0.5+border));
     /* slight random gradation on flagstones... */
	a = flag ? 0.2 : 1.0 - 0.5*mottle(x,y,z,0);
    break;
    case XBRICKS:
	size = 0.5;
	border = 0.12;
	mody = y/size;
	modz = z/size;
	mody = (y>0) ? mody-(int)mody : 1.+mody-(int)mody;
	modz = (z>0) ? modz-(int)modz : 1.+modz-(int)modz;
	if (modz<0.5) flag = modz<border || mody<border;
	else flag = modz<(0.5+border) || (mody>0.5 && mody<(0.5+border));
	a = flag ? -1.0 : 1.0;
    break;
    case YBRICKS:
	size = 0.5;
	border = 0.12;
	modx = x/size;
	modz = z/size;
	modx = (x>0) ? modx-(int)modx : 1.+modx-(int)modx;
	modz = (z>0) ? modz-(int)modz : 1.+modz-(int)modz;
	if (modz<0.5) flag = modx<border || modz<border;
	else flag = modz<(0.5+border) || (modx>0.5 && modx<(0.5+border));
	a = flag ? -1.0 : 1.0;
    break;
    case ZIMAGE:
	modx = x/imgsize[imageno];
	mody = (y*imgx[imageno])/(imgsize[imageno]*imgy[imageno]*imgasp[imageno]);
	modx = (modx>0) ? modx-(int)modx : 1.+modx-(int)modx;
	mody = (mody>0) ? mody-(int)mody : 1.+mody-(int)mody;
     /* pointers to image buffer array... */
	ix = (long)(modx*imgx[imageno]);
	iy = (long)(mody*imgy[imageno]);
	*trgb = imgr[imageno][iy][ix]/255.;
	*(trgb+1) = imgg[imageno][iy][ix]/255.;
	*(trgb+2) = imgb[imageno][iy][ix]/255.;
	a = -1.0;
    break;
    case ZSPHEREWRAP:
     /* find k, vector from center of sphere[pointer[j]]... */
	i = pointer[j];
	kx = x - cen[i][0];
	ky = y - cen[i][1];
	kz = z - cen[i][2];
     /* find spherical coords... */
	r = kx*kx + ky*ky + tiny;
	theta = (kx>0) ? asin(ky/sqrt(r)) : PI - asin(ky/sqrt(r));
	modx = theta/(imgsize[imageno]*PI) + 0.5;
	size = sqrt(r + kz*kz);
	mody = imgx[imageno]*asin(kz/size)/(imgsize[imageno]*imgy[imageno]*imgasp[imageno]*PI) + 0.5;
	modx = (modx>0) ? modx-(int)modx : 1.+modx-(int)modx;
	mody = (mody>0) ? mody-(int)mody : 1.+mody-(int)mody;
     /* pointers to image buffer array... */
	ix = (long)(modx*imgx[imageno]);
	iy = (long)(mody*imgy[imageno]);
	*trgb = imgr[imageno][iy][ix]/255.;
	*(trgb+1) = imgg[imageno][iy][ix]/255.;
	*(trgb+2) = imgb[imageno][iy][ix]/255.;
	a = -1.0;
    break;
    case ZIMGGOURAUD:
	modx = x/imgsize[imageno];
	mody = (y*imgx[imageno])/(imgsize[imageno]*imgy[imageno]*imgasp[imageno]);
	modx = (modx>0) ? modx-(int)modx : 1.+modx-(int)modx;
	mody = (mody>0) ? mody-(int)mody : 1.+mody-(int)mody;
	gouraudimage(modx,mody,trgb,imageno);
	a = -1.0;
    break;
    case ZSPHGOURAUD:
     /* find k, vector from center of sphere[pointer[j]]... */
	i = pointer[j];
	kx = x - cen[i][0];
	ky = y - cen[i][1];
	kz = z - cen[i][2];
     /* find spherical coords... */
	r = kx*kx + ky*ky + tiny;
	theta = (kx>0) ? asin(ky/sqrt(r)) : PI - asin(ky/sqrt(r));
	modx = theta/(imgsize[imageno]*PI) + 0.5;
	size = sqrt(r + kz*kz);
	mody = imgx[imageno]*asin(kz/size)/(imgsize[imageno]*imgy[imageno]*imgasp[imageno]*PI) + 0.5;
	modx = (modx>0) ? modx-(int)modx : 1.+modx-(int)modx;
	mody = (mody>0) ? mody-(int)mody : 1.+mody-(int)mody;
	gouraudimage(modx,mody,trgb,imageno);
	a = -1.0;
    break;

    }
    if (a>-tiny) {
	*trgb     = a*objrgb[j][0];
	*(trgb+1) = a*objrgb[j][1];
	*(trgb+2) = a*objrgb[j][2];
    }
}

void gouraudimage(mx,my,trgb,i)
short i;
float mx,my,*trgb;
 /* mx,my give coords on image in range 0 to 1.0, trgb points to
  * returned rgb values, i tells it which image buffer to use... */
{
    short ox,oy,ix,iy;
    float x,y,fx,fy,w00,w01,w10,w11;

    ox = (short)(x = mx*imgx[i]);
    oy = (short)(y = my*imgy[i]);
    ix = (ox+1) % imgx[i];
    iy = (oy+1) % imgy[i];
    fx = x - (float)ox;
    fy = y - (float)oy;
    w00 = (1.-fx)*(1.-fy);
    w01 = fx*(1.-fy);
    w10 = (1.-fx)*fy;
    w11 = fx*fy;
    *trgb     = (w00*imgr[i][oy][ox] + w01*imgr[i][oy][ix] +
	w10*imgr[i][iy][ox] + w11*imgr[i][iy][ix])/255;
    *(trgb+1) = (w00*imgg[i][oy][ox] + w01*imgg[i][oy][ix] +
	w10*imgg[i][iy][ox] + w11*imgg[i][iy][ix])/255;
    *(trgb+2) = (w00*imgb[i][oy][ox] + w01*imgb[i][oy][ix] +
	w10*imgb[i][iy][ox] + w11*imgb[i][iy][ix])/255;
}

void bumpmap(x,y,z,j,nx,ny,nz)
double x,y,z;
double *nx,*ny,*nz;
short j;
 /* modulates normal using texture[j][1] type, and normalises it.
  * Ideas taken from Advanced Animation and Rendering Techniques by Watt &
  * Watt, p200... */
{
    double l,f;
    float tweak,size;

    switch (texture[j][1]) {

    case FROSTED1:
	tweak = 0.04;
	*nx += tweak*random()/MAXINT;
	*ny += tweak*random()/MAXINT;
	*nz += tweak*random()/MAXINT;
    break;
    case FROSTED2:
	tweak = 0.25;
	*nx += tweak*random()/MAXINT;
	*ny += tweak*random()/MAXINT;
	*nz += tweak*random()/MAXINT;
    break;
    case FROSTED3:
	tweak = 1.0;
	*nx += tweak*random()/MAXINT;
	*ny += tweak*random()/MAXINT;
	*nz += tweak*random()/MAXINT;
    break;
    case ZPOOL1:
	tweak = 0.5;
	size = 1.0;
     /* maybe use shear xform to make ripples directional... */
	x = (x)/size;
	y = (y)/size;
	*nx += tweak * (mottle(x,y,z/size,0) - 0.5);
	*ny += tweak * (mottle(x,y,z/size,1) - 0.5);
    break;
    case ZPOOL2:
	tweak = 0.4;
	l = sqrt(x*x + y*y);
	f = exp(-(l-4)*(l-4)/2.0)*cos(12*l);
	*nx += tweak * f * x/l;
	*ny += tweak * f * y/l;
    break;
    case BUMPY1:
	size = 0.08;
	tweak = 0.5;
	*nx += tweak*(mottle(x/size,y/size,z/size,0)-0.5);
	*ny += tweak*(mottle(x/size,y/size,z/size,1)-0.5);
	*nz += tweak*(mottle(y/size,z/size,x/size,0)-0.5);
    break;
    case BUMPY2:
	size = 0.07;
	tweak = 0.8;
	*nx += tweak*2*(l = mottle(x/size,y/size,z/size,0)-0.5)*l;
	*ny += tweak*2*(l = mottle(x/size,y/size,z/size,1)-0.5)*l;
	*nz += tweak*2*(l = mottle(y/size,z/size,x/size,0)-0.5)*l;
    break;

    }
    l = sqrt((*nx)*(*nx) + (*ny)*(*ny) + (*nz)*(*nz));
    *nx /= l;
    *ny /= l;
    *nz /= l;
}

unsigned long rnd(s)
unsigned long s;
 /* a pseudo-random number generator, 0 < rnd(seed) < 2^32 . */
{
 /* quick rnd algorithm for unsigned long integers... Numerical Recipes in Fortran, p275 */

    s = s * 1664525 + 1013904223;
    return (seed = s);
}

float mottle(x,y,z,v)
float x,y,z;
short v;
 /* returns smooth, deterministic, marbly random texture ( between 0 and 1),
  * x,y,z have been scaled so mottling is of unit size...
  * v (vary) is passed to nois. */
{
    short i,j,k;
    long ix,iy,iz;
    long p0,p1,p00,p01,p10,p11,p000,p001,p010,p011;
    float fx,fy,fz,tiny = 0.0001;

 /* integer parts... highest int less than float, ok for -ve values... */
    ix = (x>-tiny) ? (long)(x+tiny) : (long)(x+tiny) - 1;
    iy = (y>-tiny) ? (long)(y+tiny) : (long)(y+tiny) - 1;
    iz = (z>-tiny) ? (long)(z+tiny) : (long)(z+tiny) - 1;
 /* fractional parts... */
    fx = x - (float)ix;
    fy = y - (float)iy;
    fz = z - (float)iz;
 /* linear interpolation... */
    p000 = nois(ix,iy,iz,v);
    p001 = nois(ix,iy,iz+1,v);
    p010 = nois(ix,iy+1,iz,v);
    p011 = nois(ix,iy+1,iz+1,v);
    p00 = (long)((nois(ix+1,iy,iz,v) - p000)*fx) + p000;
    p01 = (long)((nois(ix+1,iy,iz+1,v) - p001)*fx) + p001;
    p10 = (long)((nois(ix+1,iy+1,iz,v) - p010)*fx) + p010;
    p11 = (long)((nois(ix+1,iy+1,iz+1,v) - p011)*fx) + p011;
    p0 = (long)((p10-p00)*fy) + p00;
    p1 = (long)((p11-p01)*fy) + p01;
    return ((p1-p0)*fz + p0) / MAXINT;
}

long nois(ix,iy,iz,vary)
long ix,iy,iz;
short vary;
 /* deterministic random value using look-up arrays, given discretised coords.
  * vary allows generation of many independently random values... (two for now) */
{
    short i,j,k;
    i = (ix>=0) ? ix % NRND : (ix % NRND) + NRND;
    j = (iy>=0) ? iy % NRND : (iy % NRND) + NRND;
    k = (iz>=0) ? iz % NRND : (iz % NRND) + NRND;
 /* use pre-stored sequences to generate seed, then iterated a few times... */
    if (vary==0) rnd(seq[i][0] * seq[j][1] * seq[k][2] + 137);
    else rnd(seq[i][2] * seq[j][0] * seq[k][1] + 257);
    rnd(seed);
    rnd(seed);
 /* lose lowest bit to fit in signed range... */
    return (long)(rnd(seed) /2);
}