/*9%&LaTeX
%Copyright 1998 by Jeffrey Lassahn.  See below.
%
%create LaTeX source from this file with
% docstrip -fraytrace.tex raytrace.c

\documentclass{report}

\input{doccode.tex}

\begin{document}
\begin{titlepage}

\vspace*{3in}

{\Huge\bf Ray Tracer}

\vspace{4pt}

Version 1.1

\vspace*{3in}
{\small This document and the program it describes are 
	copyright \copyright\ 1998 by Jeffrey Lassahn.  I give permission to
	distribute them to anyone as long as they are not modified.
	Modified versions may also be distributed as long as they bear this notice
	and clearly state the name of the person who made the modification, and the
	date of the modification.  I make no warranties of any sort regarding the
	reliability or usefulness of either this document or this computer program.}
	
\end{titlepage}

\chapter{The Ray Tracer}
This program is a fairly simple ray tracer, which reads a file describing the
scene to be rendered and outputs a BMP format image file.  To use the program,
type:

\begin{center}
\textc{raytrace <input file> <output file>}
\end{center}
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define XINFINITY (1.0e30)
#define EPSILON (1.0e-9)
#define MINCOLOR 0.01

/*9
Objects in the scene are built up by union and intersection of volumes
defined by quadratic inequalities.  Each inequality has information attached
about it's surface properties.  These properties are diffuse reflection,
specular reflection, and transparency.  Each surface property has three
numbers attached to it, specifying how much light of each of the three colors
red, green, and blue is reflected or transmitted.  This gives a total of
nine numbers to describe the surface color.

Each volume is specified by an inequality:
$$ \sum_{i,j} v_i \cdot T_{ij} \cdot v_j \leq 0 $$
where $T$ is a $4\times 4$ symmetric matrix representing the volume and
$v$ is a vector whose first $3$ components are the $x$, $y$, and $z$
components of a point in space and whose 4th component is always $1$.

All calculations in the program are done using double precision floating
point.

The individual volume descriptions are collected into an array which holds
all volumes.  Each volume is also part of a circular linked list which
specifies which other volumes it should be intersected with.  The entire
scene is, therefore, this set of points:
$$ 
\bigcup_{\hbox{\it lists}} \left( 
\bigcap_{\hbox{\it volumes}}
\cal V
\right)
$$

Any combination of unions and intersections can be expressed in this form,
although it may require the repetition of some terms.
*/

struct part {
	struct part *andlist;
	double col[3];            /* color: R,G,B */
	double trans[3];          /* transparencyto R,G,B */
	double ref[3];            /* reflection of R,G,B */
	double T[4][4];
};

struct part *fulllist;
int fulllength;

/*9
The scene is illuminated by ambient light which has brightness specified by
three numbers giving the amount of red, blue, and green light, and by $0$
or more lamps which have brightness specified by three numbers, and also
a position specified by three coordinates.  The lamps are point light sources
and are not directly visible.
*/

struct lamp {
	double intens[3];
	double r[3];
};
struct lamp *light;
int lightlength;

int res;

/*9
The camera viewpoint is specified by three coordinates for the view position,
and a set of three vectors giving the $x$ axis and $y$ axis of the surface
onto which the scene is projected, and the direction from the viewpoint to the
center of the projection plane.
*/
 
double origin[3],vx[3],vy[3],vz[3];   /* viewpoint... */

/* intersect computes the smallest positive t for which v*T*v=0
    where v=dv*t+vo.  returns 0 if there is no such point. */
/*9
\section{Intersection Computations}
The basic operation performed by the ray tracer is finding the intersection
of a ray with the surface of a volume.  This operation is performed by the
\textc{intersect} and \textc{intersectboth} functions.  The difference
between the two functions is that \textc{intersect} assumes that only the
closer of the two possible intersection points is needed, while
\textc{intersectboth} provides both solutions.

A ray is represented by
$$ R_i = V0_i + t \cdot dV_i $$
where $V0$ is the starting point of the ray, $dV$ is a vector pointing
in the direction of the ray, and $R$ traces out the ray as the scalar $t$
varies from $0$ to $\infty$.

The intersection points are computed by substituting $R$ into the equation
for the volume's surface and solving for $t$:
\begin{eqnarray*}
0 &=& \sum_{i,j} R_i \cdot T_{ij} \cdot R_j \\
  &=& \sum_{i,j} (V0_i + t \cdot dV_i) \cdot T_{ij}
                \cdot (V0_j + t \cdot dV_j) \\
  &=& \sum_{i,j} V0_i \cdot T_{ij} \cdot V0_j +
                2t(V0_i \cdot T_{ij} \cdot dV_j) +
                t^2(dV_i \cdot T_{ij} \cdot dV_j)
\end{eqnarray*}
The final step above requires that $T$ be symmetric.  At this point,
$t$ can be found by applying the quadratic formula.  The actual program
must explicitly check for special cases where $t^2$ has a coefficient of
$0$ or where the equation has no real solutions (which indicates that
the ray misses the surface entirely).

The intersection points must then be tested to be sure that $t$ is positive.
*/

double intersect(double T[4][4],double dv[3],double vo[3])
{
	double a,b,c,s;
	double x[3];

	x[0]=dv[0]*T[0][0] + dv[1]*T[0][1] + dv[2]*T[0][2];
	x[1]=dv[0]*T[1][0] + dv[1]*T[1][1] + dv[2]*T[1][2];
	x[2]=dv[0]*T[2][0] + dv[1]*T[2][1] + dv[2]*T[2][2];

	a=dv[0]*x[0] + dv[1]*x[1] + dv[2]*x[2];

	b=2*(vo[0]*x[0] + vo[1]*x[1] + vo[2]*x[2]
	+dv[0]*T[3][0] + dv[1]*T[3][1] + dv[2]*T[3][2]);


	c=vo[0]*(vo[0]*T[0][0] + 2*(vo[1]*T[0][1] + vo[2]*T[0][2] + T[0][3]))
	+ vo[1]*                (vo[1]*T[1][1] + 2*(vo[2]*T[1][2] + T[1][3]))
	+ vo[2]*                                 (vo[2]*T[2][2] + 2*T[2][3])
																														 +T[3][3];


	/* FIXME do we need these EPSILONs to avoid overflow problems?
	if ((-EPSILON<a)&&(a<EPSILON)) {
		if ((-EPSILON<b)&&(b<EPSILON)) {
	*/
	if (a==0) {
		if (b==0) {
			return 0.0;
		}
		return (-c/b);
	}
	if ((b*b-4.0*a*c)<0) {
		return 0.0;
	}
	s=sqrt(b*b-4.0*a*c);
	if (((-b-s)/(2*a))>EPSILON)
		return ((-b-s)/(2*a));
	if (((-b+s)/(2*a))>EPSILON)
		return ((-b+s)/(2*a));
	return 0.0;
}

double intersectboth(double T[4][4],double dv[3],double vo[3], double *t2)
{
	double a,b,c,s;
	double x[3];

	(*t2)=0.0;

	x[0]=dv[0]*T[0][0] + dv[1]*T[0][1] + dv[2]*T[0][2];
	x[1]=dv[0]*T[1][0] + dv[1]*T[1][1] + dv[2]*T[1][2];
	x[2]=dv[0]*T[2][0] + dv[1]*T[2][1] + dv[2]*T[2][2];

	a=dv[0]*x[0] + dv[1]*x[1] + dv[2]*x[2];

	b=2*(vo[0]*x[0] + vo[1]*x[1] + vo[2]*x[2]
	+dv[0]*T[3][0] + dv[1]*T[3][1] + dv[2]*T[3][2]);

	c=vo[0]*(vo[0]*T[0][0] + 2*(vo[1]*T[0][1] + vo[2]*T[0][2] + T[0][3]))
	+ vo[1]*                  (vo[1]*T[1][1] + 2*(vo[2]*T[1][2] + T[1][3]))
	+ vo[2]*                                    (vo[2]*T[2][2] + 2*T[2][3])
																																+T[3][3];
	if ((-0.0001<a)&&(a<0.0001)) {
		if ((-0.0001<b)&&(b<0.0001)) {
			return (0.0);
		}
		return (-c/b);
	}
	if ((b*b-4.0*a*c)<0) {
		return (0.0);
	}
	s=sqrt(b*b-4.0*a*c);
	(*t2)=(-b-s)/(2*a);
	return ((-b+s)/(2*a));
}

/*9
\section{Shadow Computation}
Shadows are computed by tracing the ray from each point that is hit on each
surface to each lamp.  The ray is found by setting $V0$ to the
surface point and $dV$ to the vector from the surface to the lamp.
The amount of light that each lamp contributes to the point is the lamp's
brightness multiplied by the transparency of each surface which intersects
the ray with $0 < t \leq 1$ (it is important not to allow $t=0$ because this
will allow the object's surface to cast a shadow on itself).  
*/ 

void shadow(double dv[3],double vo[3],double mul[3])
{
	struct part *andpt;
	double t,x,y,z,t2;
	short i,f;

	for (i=0;i<fulllength;i++) {
		/* check for intersection with part i */
		t=intersectboth (fulllist[i].T,dv,vo,&t2);
		if ((t>=EPSILON) && (t<=1.0)) {
			x=vo[0]+dv[0]*t;
			y=vo[1]+dv[1]*t;
			z=vo[2]+dv[2]*t;

			f=-1;
			andpt=fulllist[i].andlist;
			while((andpt!=&fulllist[i]) && (f)) {
				if (0.0 < x*(x*(andpt->T[0][0]) 
										+2*y*(andpt->T[0][1]) 
										+2*z*(andpt->T[0][2]) 
										+2*(andpt->T[0][3]))
							 	 +y*(y*(andpt->T[1][1]) 
										+2*z*(andpt->T[1][2])
										+2*(andpt->T[1][3]))
								 +z*(z*(andpt->T[2][2])
										+2*(andpt->T[2][3]))
								 +(andpt->T[3][3]))

					f=0;
				andpt=andpt->andlist;
			}
			if (f) {
				mul[0]*=fulllist[i].trans[0];
				mul[1]*=fulllist[i].trans[1];
				mul[2]*=fulllist[i].trans[2];
			}
		}

		if ((t2>=EPSILON) && (t2<=1.0)) {
			x=vo[0]+dv[0]*t2;
			y=vo[1]+dv[1]*t2;
			z=vo[2]+dv[2]*t2;

			f=-1;
			andpt=fulllist[i].andlist;
			while((andpt!=&fulllist[i]) && (f)) {
				if (0.0 < x*(x*(andpt->T[0][0]) 
									+2*y*(andpt->T[0][1])
									+2*z*(andpt->T[0][2])
									+2*(andpt->T[0][3]))
								+y*(y*(andpt->T[1][1]) 
									+2*z*(andpt->T[1][2])
									+2*(andpt->T[1][3]))
								+z*(z*(andpt->T[2][2])
									+2*(andpt->T[2][3]))
								+(andpt->T[3][3]))
					f=0;
				andpt=andpt->andlist;
			}
			if (f) {
				mul[0]*=fulllist[i].trans[0];
				mul[1]*=fulllist[i].trans[1];
				mul[2]*=fulllist[i].trans[2];
			}
		}
	}   /* end of (for all parts) loop */
	return;
}

/* finds the intersection of the ray
    vo+t*dv (t>0)
with the world in general.
c is a pointer to the color of the result.  It must be
initialized to 0 by the calling routine, or the result will be
added to what is already there.
mul is a pointer to the multiplier for colors. It should start
initialized to 1. don't expect any inputs to be unmodified after the call.

*/
/*9
\section{Tracing a Ray}
The color of a ray is determined recursively.  The ray maintains three
numbers representing the current color (which start out as $0$) and
three numbers which are multipliers that represent how much objects
intersected by this ray can change the color of the final pixel.  The
multipliers start out as $1$.

For each pixel in the image, a ray is sent out through the pixel.  The closest
point where the ray intersects a surface is found by checking each surface
for intersection, and then comparing the intersection point with each volume
in the intersection list to see that it is part of an actual object.
This process is done by brute force; no clever algorithms are currently
implemented to speed up the process of searching all surfaces and volumes.

Once an intersection point is found, the diffuse reflection contribution to
the color is calculated.  This is done by finding the color contributions
of each lamp.  The color contribution from each lamp is:
$$ \frac{s \cdot \cos(\theta)}{r^2} $$
where $s$ is the brightness of the lamp adjusted for shadows as described
above, $r$ is the distance between the point and the lamp, and $\theta$ is
the angle between the surface's normal vector and the lamp's position
vector.  The contributions from all lamps are added together, and the sum is
multiplied by the diffuse color components of the surface and the components
of the multiplier.  The resulting color values are added to the ray color.

The contributions from specular reflection and transparency to the color of
the point are then found by creating two new rays; one at the reflection
angle and one in the same direction as the incoming ray.  The multipliers
for these rays are adjusted by multiplying with the transparency and
specular reflection terms of the surface color, and the process is repeated
by tracing the new rays.  These new rays share the color accumulators with
the original ray, so that they automatically add their color contributions to
the final pixel color.

The tree of rays this recursion produces is terminated whenever a ray's
multiplier drops below a preset limit.  If perfectly reflective surfaces
are arranged incorrectly, the multiplier for some rays might never be
reduced, resulting in an infinite recursion.  This should be avoided by
not using perfectly reflective surfaces.

Both reflection calculations use the normal vector of the surface.  This
can be found by taking the gradient of the surface equation as follows:
\begin{eqnarray*}
n &=& \nabla \sum_{i,j} v_i \cdot T_{ij} \cdot v_j \\
  &=& \frac{\partial}{\partial v_k} \sum_{i,j} v_i \cdot T_{ij} \cdot v_j \\
  &=& \sum_{i,j}\frac{\partial}{\partial v_k} v_i \cdot T_{ij} \cdot v_j \\
  &=& \sum_{i,j}\left(\frac{\partial}{\partial v_k} v_i \right) 
						\cdot T_{ij} \cdot v_j  +
  		\sum_{i,j} v_i \cdot 
  		\left(\frac{\partial}{\partial v_k} T_{ij} \cdot v_j\right)\\
  &=& \sum_i\left(\frac{\partial}{\partial v_k} v_k \right) 
						\cdot T_{kj} \cdot v_j  +
  		\sum_i v_i \cdot 
  		\left(\frac{\partial}{\partial v_k} T_{ik} \cdot v_k\right)\\
	&=& 2 \sum_i T_{kj} \cdot v_j
\end{eqnarray*}
*/

void trace( double dv[3],double vo[3],double c[3],double mul[3])
{
	struct part *andpt;
	double mint;
	double t,x,y,z,nn;
	double n[3],dvt[3],mult[3];
	short i,bestpt,f;

	/* FIXME should this be mul[0] + mul[1] + mul[2] ? */
	if ((mul[0]+mul[1]+mul[1])<MINCOLOR) {
		c[0]+=mul[0]*light[0].intens[0];
		c[1]+=mul[1]*light[0].intens[1];
		c[2]+=mul[2]*light[0].intens[2];
		return;
	}
	bestpt=fulllength;
	mint=XINFINITY;
	for (i=0;i<fulllength;i++) {
		/* check for intersection with part i */
		t=intersect (fulllist[i].T,dv,vo);
		if ((t>=EPSILON) && (t<mint)) {
			x=vo[0]+dv[0]*t;
			y=vo[1]+dv[1]*t;
			z=vo[2]+dv[2]*t;

			f=-1;
			andpt=fulllist[i].andlist;
			while((andpt!=&fulllist[i]) && (f)) {
				if (0.0 < x*(  x*(andpt->T[0][0])
											+2*y*(andpt->T[0][1])
											+2*z*(andpt->T[0][2])
											+2*(andpt->T[0][3]))
								 +y*(  y*(andpt->T[1][1]) 
											+2*z*(andpt->T[1][2])
											+2*(andpt->T[1][3]))
								 +z*(  z*(andpt->T[2][2]) 
											+2*(andpt->T[2][3]))
								 +(andpt->T[3][3]))

					f=0;
				andpt=andpt->andlist;
			}
			if (f) {
				mint=t;
				bestpt=i;
			}
		}
	}   /* end of (for all parts) loop */
	/* do color things for fullist[bestpt] and return */
	if (bestpt!=fulllength) {
		/* for convenience redefine andpt... */
		andpt=&fulllist[bestpt];
		/* coords. */
		x=mint*dv[0]+vo[0];
		y=mint*dv[1]+vo[1];
		z=mint*dv[2]+vo[2];

		/*normal vector */
		n[0]=andpt->T[0][0]*x+andpt->T[1][0]*y+andpt->T[2][0]*z+andpt->T[3][0];
		n[1]=andpt->T[0][1]*x+andpt->T[1][1]*y+andpt->T[2][1]*z+andpt->T[3][1];
		n[2]=andpt->T[0][2]*x+andpt->T[1][2]*y+andpt->T[2][2]*z+andpt->T[3][2];
		nn=n[0]*n[0]+n[1]*n[1]+n[2]*n[2];

		/* transparency */
		vo[0]=x;
		vo[1]=y;
		vo[2]=z;
		/* copy dv and mul so they don't get trashed... */
		dvt[0]=dv[0];
		dvt[1]=dv[1];
		dvt[2]=dv[2];
		mult[0]=mul[0]*andpt->trans[0];
		mult[1]=mul[1]*andpt->trans[1];
		mult[2]=mul[2]*andpt->trans[2];
		trace(dvt,vo,c,mult);

		/* diffuse color */
		vo[0]=x;
		vo[1]=y;
		vo[2]=z;

		for (i=1;i<lightlength;i++) {
			dvt[0]=(light[i].r[0]-x);
			dvt[1]=(light[i].r[1]-y);
			dvt[2]=(light[i].r[2]-z);
			mult[0]=mul[0]*light[i].intens[0]*andpt->col[0];
			mult[1]=mul[1]*light[i].intens[1]*andpt->col[1];
			mult[2]=mul[2]*light[i].intens[2]*andpt->col[2];
			shadow(dvt,vo,mult);
			if ((mult[0]+mult[1]+mult[2])>= MINCOLOR) {
				/* c+=col*intensity*shadow* |<n,l>| /(|l|*|n|*<l,l>) */
				t=n[0]*dvt[0]+n[1]*dvt[1]+n[2]*dvt[2];
				t=sqrt(t*t/(nn*(dvt[0]*dvt[0]+dvt[1]*dvt[1]+dvt[2]*dvt[2])));

	/* the next line is intensity=brightness/r^2.  It's right, but may
		not be a good idea in practice... */
				t/=(dvt[0]*dvt[0]+dvt[1]*dvt[1]+dvt[2]*dvt[2]);

				c[0]+=mult[0]*t;
				c[1]+=mult[1]*t;
				c[2]+=mult[2]*t;
			}
		}
		c[0]+=andpt->col[0]*mul[0]*light[0].intens[0];
		c[1]+=andpt->col[1]*mul[1]*light[0].intens[1];
		c[2]+=andpt->col[2]*mul[2]*light[0].intens[2];

		/* reflection */
		vo[0]=x;
		vo[1]=y;
		vo[2]=z;
		/* dv <= dv-2*n(<n,dv>/<n,n>)  */
		t=(n[0]*dv[0]+n[1]*dv[1]+n[2]*dv[2])/nn;
		dv[0]-=2*n[0]*t;
		dv[1]-=2*n[1]*t;
		dv[2]-=2*n[2]*t;
		mul[0]*=andpt->ref[0];
		mul[1]*=andpt->ref[1];
		mul[2]*=andpt->ref[2];
		trace(dv,vo,c,mul);
	}
	else {
		c[0]+=mul[0]*light[0].r[0];
		c[1]+=mul[1]*light[0].r[1];
		c[2]+=mul[2]*light[0].r[2];
	}
	return;
}


/*9
\section{File Format for Scene Models}
The format for files which describe a scene to be rendered is quite simple,
in the sense that it is easy for the machine to process.  Unfortunately,
this means that these files are not very easy for the user to read.
Each scene is specified by a single disk file, which consists of nothing
but white space and numbers in decimal format, except for the resolution
specifier which is in hexadecimal.  The white space is ignored, which
means that the numbers may be indented and arranged into lines in whatever
way is most convenient.

The file layout is as follows:

\begincode
<number of lamps + 1>
<number of surfaces>

Background color: <red> <blue> <green>
Ambient light: <red> <blue> <green>

Repeat once for each lamp:
	Position: <x> <y> <z>
	Color: <red> <blue> <green>

Repeat once for each group of intersecting volumes:
	<number of volumes to intersect>
	Repeat for each volume:
		Diffuse reflection color: <red> <green> <blue>
		transparency color: <red> <green> <blue>
		Specular reflection color: <red> <green> <blue>

		Surface Tensor:
			<Txx> <Txy> <Txz> <Tx>
			<Txy> <Tyy> <Tyz> <Ty>
			<Txz> <Tyz> <Tzz> <Tz>
			<Tx>  <Ty>  <Tz>  <Tc>

ViewPoint: <x> <y> <z>
Image Plane:
    <Vxx> <Vxy> <Vxz>
    <Vyx> <Vyy> <Vyz>
    <Vzx> <Vzy> <Vzz>

<Resolution>
\endcode

The image produced by the ray tracer is always square.  The width of the
image in pixels is $2^{16}/r$ where $r$ is the resolution parameter at the
end of the scene file.
*/

void freeit(void)
{
	if (fulllist)
		free(fulllist);
	if (light)
		free(light);
}


short loadit(char s[])
{
int i,j,hold,x,y;
FILE *f;

	if(!(f=fopen(s,"r")))
			return 0;

	fulllength=0;
	fulllist=NULL;
	if (!fscanf(f,"%d",&lightlength)) {
		fclose(f);
		return 0;
	}
	fscanf(f,"%d",&fulllength);
	if (!fulllength) {
		return 0;
		fclose(f);
	}
	fulllist=malloc(((long)fulllength)*sizeof(struct part));
	if (!fulllist) {
		fclose(f);
		return 0;
	}
	light=malloc(((long)lightlength)*sizeof(struct lamp));
	if (!light) {
		fclose(f);
		freeit();
		return 0;
	}

	for (i=0;i<lightlength;i++) {
		for (j=0;j<3;j++)
			fscanf(f,"%lf",&light[i].r[j]);
		for (j=0;j<3;j++)
			fscanf(f,"%lf",&light[i].intens[j]);
	}

	i=0;
	while (i<fulllength) {
		if (!fscanf(f,"%d",&j)) {
			printf("Malformed file!\n");
			fclose(f);
			return 0;
		}
		j+=i;
		hold=i;
		if (j>fulllength) {
			fclose(f);
			printf("Malformed file!\n");
			return 0;
		}
		while (i<j) {
			fulllist[i].andlist=&fulllist[i+1];
			/* read data */
			for (x=0;x<3;x++)
				fscanf(f,"%lf",&fulllist[i].col[x]);

			for (x=0;x<3;x++)
				fscanf(f,"%lf",&fulllist[i].trans[x]);

			for (x=0;x<3;x++)
				fscanf(f,"%lf",&fulllist[i].ref[x]);

			for (y=0;y<4;y++)
				for (x=0;x<4;x++)
					fscanf(f,"%lf",&fulllist[i].T[x][y]);

			i++;
		}
		fulllist[i-1].andlist=&fulllist[hold];
	}

	for (i=0;i<3;i++)
		fscanf(f,"%lf",&origin[i]);

	for (i=0;i<3;i++)
		fscanf(f,"%lf",&vx[i]);
	for (i=0;i<3;i++)
		fscanf(f,"%lf",&vy[i]);
	for (i=0;i<3;i++)
		fscanf(f,"%lf",&vz[i]);

	fscanf(f,"%x",&res);
	fclose (f);

	return -1;
}



void printpart(struct part *s)
{
short x,y;
	printf("@ %ld\n",s);
	printf("next @ %ld\n",s->andlist);
	printf("color: %lf %lf %lf\n",
	s->col[0],s->col[1],s->col[2]);
	for (y=0;y<4;y++) {
		for (x=0;x<4;x++)
			printf ("%lf ",s->T[x][y]);
		printf ("\n");
	}
}

#define HEADSIZE 54

unsigned char bmphead[]={
	'B','M',
	0,0,0,0,
	0,0,0,0,
	HEADSIZE,0,0,0,
	0x28,0,0,0,
	0,0,0,0,
	0,0,0,0,
	1,0,
	24,0,
	0,0,0,0,
	0,0,0,0,
	0,0,0,0,
	0,0,0,0,
	0,0,0,0,
	0,0,0,0	
};

int WriteBMPHead(FILE *fp,int dx,int dy,int bits)
{
long t,bytes;

	if (dx<1)
		return 0;
	if (dy<1)
		return 0;

	if (bits!=24)
		return 0;
	bytes=3;

	t=(long)dx*(long)dy*bytes;
	/* file size */
	bmphead[2]=(t+HEADSIZE)&0xFF;
	bmphead[3]=((t+HEADSIZE)>>8)&0xFF;
	bmphead[4]=((t+HEADSIZE)>>16)&0xFF;
	bmphead[5]=((t+HEADSIZE)>>24)&0xFF;

	/* width */
	bmphead[18]=dx&0xFF;
	bmphead[19]=(dx>>8)&0xFF;
	bmphead[20]=(dx>>16)&0xFF;
	bmphead[21]=(dx>>24)&0xFF;

	/* height */
	bmphead[22]=dy&0xFF;
	bmphead[23]=(dy>>8)&0xFF;
	bmphead[24]=(dy>>16)&0xFF;
	bmphead[25]=(dy>>24)&0xFF;

	/* image size */
	0,0,0,0,
	bmphead[34]=t&0xFF;
	bmphead[35]=(t>>8)&0xFF;
	bmphead[36]=(t>>16)&0xFF;
	bmphead[37]=(t>>24)&0xFF;

	if (fwrite(bmphead,1,sizeof(bmphead),fp)!=sizeof(bmphead)) {
		puts("bad write");
		return 0;
	}
	return 1;
}

void writeheader(FILE *f) {
	int width;
	width=(0xFFFDL+res)/res;
	WriteBMPHead(f,width,width,24);
}

int main(int argc,char *argv[])
{
double cols[3],muls[3];
long x,y;
double vo[3],dv[3];
FILE *f;

	if (argc<2) {
		printf("no command line args\n");
		exit(-1);
	}
	if (!loadit(argv[1])) {
		printf("we've got a problem!\n");
		exit(-1);
	}


	if (argc>2)
			f=fopen(argv[2],"wb");
	else {
		f=NULL;
		printf("no output file specified\n");
		exit(-1);
	}

	if (f)
		writeheader(f);
	else {
		printf("Couldn't open output file\n");
		exit(-1);
	}

	for (y=-0x7FFFL;y<=0x7FFFL;y+=res)
	for (x=-0x7FFFL;x<=0x7FFFL;x+=res) {
		vo[0]=origin[0];
		vo[1]=origin[1];
		vo[2]=origin[2];

		dv[0]=vx[0]*((double)x)/0x7FFFL+vy[0]*((double)y)/0x7FFFL+vz[0];
		dv[1]=vx[1]*((double)x)/0x7FFFL+vy[1]*((double)y)/0x7FFFL+vz[1];
		dv[2]=vx[2]*((double)x)/0x7FFFL+vy[2]*((double)y)/0x7FFFL+vz[2];

		cols[0]=0.0;
		cols[1]=0.0;
		cols[2]=0.0;
		muls[0]=1.0;
		muls[1]=1.0;
		muls[2]=1.0;

		trace(dv,vo,cols,muls);

		if (cols[0]>1.0) cols[0]=1.0;
		if (cols[1]>1.0) cols[1]=1.0;
		if (cols[2]>1.0) cols[2]=1.0;

		if (f) {
			putc((int)(cols[2]*255),f);
			putc((int)(cols[1]*255),f);
			putc((int)(cols[0]*255),f);
		}
	}
	fclose(f);
	freeit();
	return 0;
}

/*9
\end{document}
*/