/*************************************************************************
**									**
**	TRACE1.C	A wee little ray-tracing program for IBM-XT,	**
**			Cubicomp CS-5, and Microsoft-C.			**
**									**
**	K.A.Bjorke	10 July 1984					**
**	(C) 1984	5th Generation Digital Research			**
**									**
*************************************************************************/

/*
.fo 			/* TRACE1.C Page # */
(dot commands in some comments for listing with wordstar)
*/

#include <stdio.h>
#include <conio.h>	/* All I/O current thru monitors */

#define VERSION "1.00"

#define LIGHTSOURCE 1	/* enumerated object types */
#define GLASS 2
#define PLANE 3

#define MAXLITES 4	/* maximum Number of each enumerated type above */
#define MAXBALLS 6
#define MAXPOLYS 8

#define NULL 0
#define FALSE 0
#define TRUE 1
#define ZERO 0.0

#define AIR 0.985	/* per-unit transmission in open air (close to 1) */
#define INFINITY 65565	/* maximum scene-space value */
#define MAXLEVEL 10	/* maximum levels of recursive tracing */
#define BRILLIANCE transmit	/* to aid in "cheating" */
#define BACKGROUND ZERO	/* intensity if there's nothing there */

#define XREG 0x0300	/* Cubicomp Control-Register Ports */
#define YREG 0x0302
#define DATA 0x0304
#define RED 0x0306
#define GREEN 0x0308
#define BLUE 0x030A
#define CONTROL 0x030C
#define COMMAND 0x030E

#define MAXSCAN 512	/* miscellaneous Machine Equates */
#define CENTER (MAXSCAN/2)
#define MAXGREY 256	/* Number of gray scales for monochrome picture */

/*
.pa 
*/

/*********************************
** 	External Functions 	**
*********************************/

extern sqrt();	/* pointer to float as argument!! */

/*****************************************
** 	Global structures & types 	**
*****************************************/

/* typedef int (*PFI)(); */	   /* used in plane structure definition */
typedef float VECTOR[3];	   /* used a lot */

struct ball {
	int center[3];
	int radius, r2;		   /* r2 = radius * radius */
	float k, spec;		   /* refract, reflect */
	float transmit;		   /* cheat and use transmit as intensity */
} sphere[MAXBALLS], light[MAXLITES];

struct beam {
	VECTOR origin;
	VECTOR direction;
	float ac, bc, cc, dc;	   /* cc & dc are PARTIAL quadratic terms */
};

struct plane {
	float aco, bco, cco, dco;  /* as in ax+by+cz+d=0 */
	float diff;		   /* assume only diffuse reflections */
	VECTOR normal;		   /* normal-to-plane */
     /*	PFI style; */		   /* function which would define bounds */
} poly[MAXPOLYS];

struct strike {
	float t;
	int type;
	int index;
};

/************************************
** 	Other Global Variables 	   **
************************************/

VECTOR viewpoint;			/* camera position */
int lightsnow, ballsnow, polysnow;	/* how many for this pic? */ 
float intamb;				/* ambient light */
float greys[MAXSCAN+1];			/* used for anti-aliasing */
/* FILE *fpt; */			/* general-purpose file pointer */
/* int _fmode = 0x8000;	*/		/* binary file mode */
char title[32];				/* title for picture */
int px, py;				/* current screen coordinates */
/*
.pa 
*/

/*********************************
**				**
**	Code Begins Here	**
**				**
*********************************/

main(argc,argv)
	int argc;
	char *argv[];
{
	advertise();		/* (C) 1984 */
	if (argc != 1)
		tutor();	/* arguments? what? */
	/* else... */
	read_script();		/* set up scene */
	grey_map();		/* create lookup table */
	init_screen();		/* clear screen and set for continuous write */
	draw_pic();		/* trace like crazy */
}

/***********************
** Vector dot product **
***********************/
float dotprod(a,b)
	float *a, *b;
{
	int i;
	float q;
	q = ZERO;
	for (i = 0; i < 3; ++i)
		q += (*a++ * (*b++));
	return (q);
}

/*
.pa 
*/

/*********************************
** Recursive Ray-Trace Function **
*********************************/
float trace_ray(zap, medium, level)
	struct beam *zap; 
	float medium;
	int level;
{
	int i, j, m, ct;
	float med2, result, falloff;
	VECTOR util;		/* utility vector */
	struct strike intersect[MAXLITES+MAXBALLS+MAXPOLYS+1];
	struct beam newray;	/* created by GLASS objects */

	result = ZERO;
	if ((--level == 0) || (medium == ZERO))
		return(result);		/* to prevent endless reflections */
					/* & skip opaque objects */
	ct = 0;
	prepare_beam(zap);
	blank_intersects(&intersect, (MAXLITES+MAXBALLS+MAXPOLYS));
		/* now, search for intersections with all objects... */
	for (i = 0; i < lightsnow; ++i) {
		if (cross_ball(&light[i], zap, &intersect[ct].t) != FALSE) {
			intersect[ct].index = i;
			intersect[ct++].type = LIGHTSOURCE;
		}
	}
	for (i = 0; i < ballsnow; ++i) {
		if (cross_ball(&sphere[i], zap, &intersect[ct].t) != FALSE) {
			intersect[ct].index = i;
			intersect[ct++].type = GLASS;
		}
	}
	for (i = 0; i < polysnow; ++i) {
		if (cross_poly(&poly[i], zap, &intersect[ct].t) != FALSE) {
			intersect[ct].index = i;
			intersect[ct++].type = PLANE;
		}
	}
	if (ct == 0)
		return((result = BACKGROUND));		/* nuthin */
	/* else... */
	i = sort_by_t(&intersect);		/* now pick nearest */
	j = intersect[i].index;			/* for convenience */
	falloff = medium / intersect[i].t;	/* this can get used a lot */
/* [Continued]
.pa 
*/

	if (intersect[i].type == LIGHTSOURCE) {   /* simplest case */
		result = light[j].BRILLIANCE * falloff;
	} else if (intersect[i].type == GLASS) {
			/* going in or out of sphere? */
		med2 = (sphere[j].transmit == medium) ?
			AIR : sphere[j].transmit;
			/* calculate REFRACTION */
		for (m = 0; m < 3; ++m) 
			util[m] = sphere[j].center[m] -
				    (newray.origin[m] = (zap->origin[m] + 
				    intersect[i].t * (zap->direction[m])));
		normalize(util);	/* used as normal vector */
			/* note normal goes INTO the sphere */
			/* this is because zap does, too */
		for (m = 0; m < 3; ++m)
			newray.direction[m] = zap->direction[m] + 
					      sphere[j].k * util[m];
				/* deflect newray in direction of normal */
		result = trace_ray(&newray, med2, level) * falloff;

			/* now add REFLECTION */
		med2 = 2 * dotprod(util, &zap->direction[0]); 
				/* N.L = cos(angle_of_incidence) */
		for (m = 0; m < 3; ++m)
			newray.direction[m] = zap->direction[m] - 
					      util[m] * med2;
				/* for regular vectors, incid + reflect = */
				/* normal * 2 cos(angle_of_incidence) */
				/* note reversal of signs! */
		result += trace_ray(&newray, medium, level) * 
			   sphere[j].spec * falloff;
	} else {		/* must be a plane, then */
		result = intamb * poly[j].diff * falloff;
		for (m = 0; m < 3; ++m)
			util[m] = zap->origin[m] + 
				    intersect[i].t * zap->direction[m];
		for (m = 0; m < lightsnow; ++m) {
			visible(m, i, util, medium, falloff, &result);
		}
	}
	return(result);
}

/*
.pa
*/

/************************
** Actual drawing loop **
************************/
draw_pic()
{
	int i;
	struct beam laser;		/* to stress point sampling */
	float intensity, prev, avg;

	cprintf("\n\n\n\nNow Drawing %s.\n",title);
	for (i = 0; i < 3; ++i)
		laser.origin[i] = viewpoint[i];		/* consistent */
	for (py = 0; py <= MAXSCAN; ++py) {		/* scan rows */
		for (px = 0; px <= MAXSCAN; ++px) {	/* scan columns */
			prev = intensity;
			laser.direction[0] = px - CENTER;
			laser.direction[1] = py - CENTER;
			laser.direction[2] = -viewpoint[2];
				/* now do all the real work: */
			intensity = trace_ray(&laser, AIR, MAXLEVEL);
			if ((px != 0) && (py != 0)) {
				avg = (intensity + prev + greys[px] +
				       greys[px - 1]) / 4;
				if (avg > (MAXGREY-1))
					avg = (MAXGREY-1);
				write_pixel((px-1), (py-1), (int)avg);
					/* avg is anti-aliased */
			}
			greys[((px == 0) ? MAXSCAN : (px-1))] = prev;
		}
	}
}

/*
.pa 
*/

/***********************************************************************
** find out if a lightsource is visible from some point and add it in **
***********************************************************************/
int visible(i,f,point, medium, falloff, value)
	int i,f;			/* light & poly indices */
	float point[];			/* point of incidence */
	float medium;			/* local transmission value */
	float falloff;			/* falloff to origin of first ray */
	float *value;			/* cumulative intensity */
{
	int m, j, ct;
	float dist;
	struct strike testint[MAXBALLS+MAXPOLYS+1];
	struct beam testray;
		/* we only have to check intersects of light-stoppers */

	ct = 0;
	blank_intersects(testint,MAXLITES);
		/* point testray at light[i] */
	for (m = 0; m < 3; ++m) {
		testray.origin[m] = point[m];
		testray.direction[m] = light[i].center[m] - point[m];
	}
	prepare_beam(&testray);
	cross_ball(&light[i], &testray, &dist);	/* get distance to light */
	for (j = 0; j < ballsnow; ++j) {
		if (cross_ball(&sphere[j], &testray, &testint[ct].t) != NULL) {
			testint[ct].index = j;
			testint[ct++].type = GLASS;
		}
	}
	for (j = 0; j < polysnow; ++j) {
		if (cross_poly(&poly[j], &testray, &testint[ct].t) != NULL) {
			testint[ct].index = j;
			testint[ct++].type = PLANE;
		}
	}
		/* if (ct == 0) then trivial accept -- all clear */
	if (ct != 0) {
		j = sort_by_t(testint);
		if (testint[j].t < dist)
			return(FALSE);	/* something between us and light */
	}
	/* light is visible, so... */
	/* Lambert: ir = ia*ka+ip*kd(L.N) */
	*value += (light[i].BRILLIANCE * poly[f].diff * 
		  dotprod(&poly[f].normal[0], &testray.direction[0]) *
		  medium / dist) * falloff;
		  	 /* note no assumption of AIR made */
	return(TRUE);
}

/*
.pa
*/

/**********************************************************
** Normalize direction vector & Pre-calculate Quadratics **
**********************************************************/
prepare_beam(illum)
	struct beam *illum;
{
	int i;
	illum->ac = illum->bc = illum->cc = 0;
	normalize(&illum->direction[0]);
	for (i = 0; i < 3; ++i) {
		illum->ac += illum->direction[i] * illum->direction[i];
		illum->bc += illum->direction[i] * illum->origin[i];
		illum->cc += illum->origin[i] * illum->origin[i];
	}
	illum->bc *= 2;	/* quadratics will be the same for every ball struc */
	illum->ac *= 2;	/* only change in equation will be center[] & r2... */
	illum->dc = illum->bc * illum->bc;
}

/*
.pa
*/

/*********************************************************************
** calculate ray-ball intersections - return TRUE if any valid ones **
*********************************************************************/
int cross_ball(globe, ray, t1)
	struct ball *globe;
	struct beam *ray;
	float *t1;
{
	float c, d, t2;
	c = ray->cc - globe->r2;
	d = ray->dc - 2 * ray->ac * c;	/* quadratic determinant */
	if (d < 0)
		return (FALSE);		/* no real solutions */
	/* else... */
	if (d == 0) {			/* one real solution */
		*t1 = -(ray->bc / ray->ac);
		/* return ( (*t1 <= 0) ? FALSE : TRUE); */
		if (*t1 <= 0)
			return(FALSE);	/* but in the wrong direction */
		/* else... */
		return(TRUE);		/* if it's okay */
	}
	/* else two real solutions */
	sqrt(&d);			/* now in the form: */
					/* - b +- sqrt(d)/(2*a) */
	*t1 = -(ray->bc + d) / ray->ac;
	t2  = -(ray->bc - d) / ray->ac;
	if ((t2 <= 0) && (*t1 <= 0))
		return(FALSE);		/* both behind current position */
	if (t2 <= 0)
		return(TRUE);		   /* t1 must be okay */
	*t1 = (*t1 <= 0) ? t2 :		   /* t2 must be okay */
		           min(*t1, t2);   /* else just take the closer one */
	return(TRUE);
}

/*
.pa
*/

/***********************************************************************
** calculate ray-plane intersections - return TRUE if valid one found **
***********************************************************************/
int cross_poly(flat, ray, t)
	struct plane *flat;
	struct beam *ray;
	float *t;
{
	float dt;
	dt = flat->aco * ray->direction[0] +
	     flat->bco * ray->direction[1] +
	     flat->cco * ray->direction[2];
	if (dt == 0)
		return (FALSE);		/* ray and plane are parallel */
	/* else */
	/* simple linear math */
	*t = -( flat->aco * ray->origin[0] +
		flat->bco * ray->origin[1] +
		flat->cco * ray->origin[2] +
		flat->dco) / dt;
	return(TRUE);
}

/*
.pa 
*/

/**************************
** Normalize a 3d vector **
**************************/
normalize(vect)
	float vect[];
{
	int i;
	float magnitude;
	magnitude = 0;
	for (i = 0; i <3; ++i)
		magnitude += vect[i] * vect[i];
	sqrt(&magnitude);	/* length of [a b c] = sqrt(a^2 + b^2 + c^2) */
	for (i = 0; i < 3; ++i)
		vect[i] /= magnitude;
}

/*
.pa
*/

/********************************************
** find the nearest intersection in a list **
********************************************/
int sort_by_t(intersect)
	struct strike intersect[];
{
	int i, j;
	float t2;
	t2 = INFINITY;
	for (i = 0; intersect[i].type != NULL; ++i)
		j = (intersect[i].t < t2) ? i : j;
	return(j);
}

/*********************************************************
** blank out types in an intersection list (initialize) **
*********************************************************/
blank_intersects(intersect, top)
	struct strike intersect[];
	int top;
{
	int i;
	for (i = 0; i <= top; ++i)
		intersect[i].type = NULL;
}

/*******************************************
** set up the Cubicomp screen for drawing **
*******************************************/
init_screen()
{
	outw(CONTROL, 0x0000);		/* frame 0, no masks, normal */
	outw(DATA, 0);
	outp(CONTROL, 0x0010);		/* clear screen to black */
	while ((inp(CONTROL+1) && 0x0080) == NULL); 	/* wait for ready */
	outw(XREG, 0);
	outw(YREG, 0);		/* point to upper-left corner for counting */
}

/*******************************************************
** write a value to the current pixel & count forward **
*******************************************************/
write_pixel(x, y, value)
	int x, y, value;	/* this routine actually throws away x & y */
				/* since the Cubicomp auto-increments */
{
	outw(DATA, value);
	outp(CONTROL, 0x0021);	/* write and increment */
}

/*
.pa
*/

/*****************************************
** Make a grey-scale map, 0-(MAXGREY-1) **
*****************************************/
grey_map()
{
	int i;
	for (i = 0; i < MAXGREY; ++i) {
		outp(RED, i);			/* same value to all regs */
		outp(GREEN, i);
		outp(BLUE, i);
		outw(DATA, i);			/* and as table address */
		outp(CONTROL, 4);		/* Map-write */
		while((inp(CONTROL+1) && 0x0080) == NULL);	/* wait... */
	}
}

/********************************************************
** send a 16-bit value thru the IBM's 8-bit data lines **
********************************************************/
outw(port,wrd)
	int port, wrd;
{
	outp(port,(wrd & 0x00FF));
	outp(++port,(wrd >> 8));
}

/************************
** toot one's own horn **
************************/
advertise()
{
	int i;
	for (i = 0; i < 24; ++i)
		putch('\n',stderr);	/* crude clear-monitor */
	cputs("Trace1\tVersion ");
	cputs(VERSION);
	cputs("\tK. A. Bjorke\n");
	cputs("(C) 1984 5th Generation Digital\n\n");
}

/*
.pa
*/
	
/**********************************
** respond to a bad command-line **
**********************************/
tutor()
{
	cputs("Correct Command Line Format:\n");
	cputs("\tTRACE1 [No Args]\n");
	exit(1);
}

/*****************************************************
** Read script file (from keyboard in this version) **
*****************************************************/
read_script()
{
	int i, j;
	char xyz[3];
	xyz[0] = 'X'; xyz[1] = 'Y'; xyz[2] = 'Z';
	lightsnow = ballsnow = polysnow = 0;
	cputs("Enter Title of Image: ");
	cgets(title);
	cputs("\nEnter \'Camera\' Viewpoint as XYZ triplet:\n");
	for (j = 0; j < 3; ++j) {
		putch(xyz[j]);
		cputs(": ");
		cscanf("%d", &viewpoint[j]);
	}
	cputs("\nEnter the Intensity of Ambient Light: ");
	cscanf("%f",intamb);
	cputs("\nNumber of Light Sources: ");
	cscanf("%d", &lightsnow);
	for (i = 0; i < lightsnow; ++i) {
		cprintf("\n\nLIGHT SOURCE %1d:",(i + 1));
		cputs("\nEnter Center Point as XYZ triplet:\n");
		for (j = 0; j < 3; ++j) {
			putch(xyz[j]);
			cputs(": ");
			cscanf("%d", &light[i].center[j]);
		}
		cputs("\nRadius: ");
		cscanf("%d", &light[i].radius);
		cputs("\nIntensity: ");
		cscanf("%f", &light[i].BRILLIANCE);
		light[i].r2 = light[i].radius * light[i].radius;
	}
/*		[Continued]
.pa
*/
	cputs("\nNumber of \'Glass\' Spheres: ");
	cscanf("%d", &ballsnow);
	for (i = 0; i < lightsnow; ++i) {
		cprintf("\n\nSPHERE %1d:",(i + 1));
		cputs("\nEnter Center Point as XYZ triplet:\n");
		for (j = 0; j < 3; ++j) {
			putch(xyz[j]);
			cputs(": ");
			cscanf("%d", &sphere[i].center[j]);
		}
		cputs("\nRadius: ");
		cscanf("%d", &sphere[i].radius);
		cputs("\nIndex of Reflection: ");
		cscanf("%f", &sphere[i].spec);
		cputs("\nIndex of Refraction: ");
		cscanf("%f", &sphere[i].k);
		cputs("\nTransmission per unit of material: ");
		cscanf("%f", &sphere[i].transmit);
		sphere[i].r2 = sphere[i].radius * light[i].radius;
	}
	cputs("\n\nEnter Number of Arbitrary Planes: ");
	cscanf("%d",&polysnow);
	for (i = 0; i < polysnow; ++i) {
		cprintf("\n\nPLANE %1d\n",(i + 1));
		cputs("\tEnter coefficients of the form:\n");
		cputs("\t  Ax + By + Cz + D = 0"); 
		cputs("\nA: ");
		cscanf("%f", &poly[i].aco);
		cputs("\nB: ");
		cscanf("%f", &poly[i].bco);
		cputs("\nC: ");
		cscanf("%f", &poly[i].cco);
		cputs("\nD: ");
		cscanf("%f", &poly[i].dco);
		cputs("\n Plane Grey-Scale Value: ");
		cscanf("%f", &poly[i].diff);
		poly[i].normal[0] = poly[i].aco;
		poly[i].normal[1] = poly[i].bco;
		poly[i].normal[2] = poly[i].cco;
		normalize(&poly[i].normal[0]);	/* precalculate normal */
	}
}

/*
.pa
*/

/***************************
** If file not found, etc **
***************************/
nofile(s)
	char s[];
{
	cputs("Error! - Can't open ");
	cputs(s);
	cputs("!\n");
	exit(2);
}


/* eof */