#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#		"End of archive 17 (of 19)."
# Contents:  libray/libobj/blob.c rayview/glmethods.c
# Wrapped by kolb@woody on Wed Jul 17 17:57:00 1991
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'libray/libobj/blob.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'libray/libobj/blob.c'\"
else
echo shar: Extracting \"'libray/libobj/blob.c'\" \(18207 characters\)
sed "s/^X//" >'libray/libobj/blob.c' <<'END_OF_FILE'
X/*
X * blob.c
X *
X * Copyright (C) 1990, 1991, Mark Podlipec, Craig E. Kolb
X * All rights reserved.
X *
X * This software may be freely copied, modified, and redistributed
X * provided that this copyright notice is preserved on all copies.
X *
X * You may not distribute this software, in whole or in part, as part of
X * any commercial product without the express consent of the authors.
X *
X * There is no warranty or other guarantee of fitness of this software
X * for any purpose.  It is provided solely "as is".
X *
X * $Id: blob.c,v 4.0 91/07/17 14:36:07 kolb Exp Locker: kolb $
X *
X * $Log:	blob.c,v $
X * Revision 4.0  91/07/17  14:36:07  kolb
X * Initial version.
X * 
X */
X#include "geom.h"
X#include "blob.h"
X
Xstatic Methods *iBlobMethods = NULL;
Xstatic char blobName[] = "blob";
X
Xunsigned long BlobTests, BlobHits;
X
X/*
X * Blob/Metaball Description
X *
X * In this implementation a Blob is made up of a threshold T, and a 
X * group of 1 or more Metaballs.  Each metaball 'i'  is defined by 
X * its position (xi,yi,zi), its radius ri, and its strength si.
X * Around each Metaball is a density function di(Ri) defined by:
X *
X *         di(Ri) =  c4i * Ri^4  +  c2i * Ri^2  + c0i     0 <= Ri <= ri
X *         di(Ri) =  0                                    ri < Ri 
X *
X * where
X *       c4i  = si / (ri^4)
X *       c2i  = -(2 * si) / (ri^2)
X *       c0i  = si
X *       Ri^2 is the distance from a point (x,y,z) to the center of
X *            the Metaball - (xi,yi,zi).
X *
X * The density function looks sort of like a W (Float-u) with the 
X * center hump crossing the d-axis at  si  and the two humps on the
X * side being tangent to the R-axis at  -ri  and  +ri. Only the 
X * range  [0,ri] is being used. 
X * I chose this function so that derivative = 0  at the points  0 and +ri. 
X * This keeps the surface smooth when the densities are added. I also 
X * wanted no  Ri^3  and  Ri  terms, since the equations are easier to 
X * solve without them. This led to the symmetry about the d-axis and
X * the derivative being equal to zero at -ri as well.
X *
X * The surface of the Blob is defined by:
X *
X *
X *                  N
X *                 ---  
X *      F(x,y,z) = \    di(Ri)  = T
X *                 /
X *                 ---
X *                 i=0
X *
X * where
X *
X *     di(Ri)    is   given above
X *     Ri^2      =    (x-xi)^2  +  (y-yi)^2  + (z-zi)^2
X *      N        is   number of Metaballs in Blob
X *      T        is   the threshold
X *    (xi,yi,zi) is   the center of Metaball i
X *
X */
X
X/*****************************************************************************
X * Create & return reference to a metaball.
X */
XBlob *
XBlobCreate(T, mlist, npoints)
XFloat T;
XMetaList *mlist;
Xint npoints;
X{
X	Blob *blob;
X	int i;
X	MetaList *cur;
X
X/* 
X * There has to be at least one metaball in the blob.
X * Note: if there's only one metaball, the blob 
X * will be a sphere.
X */
X	if (npoints < 1)
X	{
X		RLerror(RL_WARN, "blob field not correct.\n");
X		return (Blob *)NULL;
X	}
X
X/*  
X * Initialize primitive and Geom structures
X */
X	blob = (Blob *)Malloc(sizeof(Blob));
X	blob->T = T;
X	blob->list=(MetaVector *)
X	    Malloc( (unsigned)(npoints*sizeof(MetaVector)) );
X	blob->num = npoints;
X
X/*
X * Initialize Metaball list
X */
X	for(i=0;i<npoints;i++)
X	{
X		cur = mlist;
X		if ( (cur->mvec.c0 < EPSILON) || (cur->mvec.rs < EPSILON) )
X		{
X			RLerror(RL_WARN,"Degenerate metaball in blob (sr).\n");
X			return (Blob *)NULL;
X		}
X		/* store radius squared */
X		blob->list[i].rs = cur->mvec.rs * cur->mvec.rs;
X		/* Calculate and store coefficients for each metaball */
X		blob->list[i].c0 = cur->mvec.c0;
X		blob->list[i].c2 = -(2.0 * cur->mvec.c0) / blob->list[i].rs;
X		blob->list[i].c4 = cur->mvec.c0 /
X				(blob->list[i].rs * blob->list[i].rs);
X		blob->list[i].x = cur->mvec.x;
X		blob->list[i].y = cur->mvec.y;
X		blob->list[i].z = cur->mvec.z;
X		mlist = mlist->next;
X		free((voidstar)cur);
X	}
X	/* 
X	 * Allocate room for Intersection Structures and
X	 * Allocate room for an array of pointers to point to
X	 * the Intersection Structures.
X	 */
X	blob->ilist = (MetaInt *)Malloc( 2 * blob->num * sizeof(MetaInt));
X	blob->iarr = (MetaInt **)Malloc( 2 * blob->num * sizeof(MetaInt *));
X	return blob;
X}
X
XMethods *
XBlobMethods()
X{
X	if (iBlobMethods == (Methods *)NULL) {
X		iBlobMethods = MethodsCreate();
X		iBlobMethods->create = (GeomCreateFunc *)BlobCreate;
X		iBlobMethods->methods = BlobMethods;
X		iBlobMethods->name = BlobName;
X		iBlobMethods->intersect = BlobIntersect;
X		iBlobMethods->normal = BlobNormal;
X		iBlobMethods->bounds = BlobBounds;
X		iBlobMethods->stats = BlobStats;
X		iBlobMethods->checkbounds = TRUE;
X		iBlobMethods->closed = TRUE;
X	}
X	return iBlobMethods;
X}
X
X/*****************************************************************************
X * Function used by qsort() when sorting the Ray/Blob intersection list
X */
Xint
XMetaCompare(A,B)
Xchar *A,*B;
X{
X	MetaInt **AA,**BB;
X
X	AA = (MetaInt **) A;
X	BB = (MetaInt **) B;
X	if (AA[0]->bound == BB[0]->bound) return(0);
X	if (AA[0]->bound <  BB[0]->bound) return(-1);
X	return(1);  /* AA[0]->bound is > BB[0]->bound */
X}
X
X/*****************************************************************************
X * Ray/metaball intersection test.
X */
Xint
XBlobIntersect(blob, ray, mindist, maxdist)
XBlob *blob;
XRay *ray;
XFloat mindist, *maxdist;
X{
X	double c[5], s[4];
X	Float dist;
X	MetaInt *ilist,**iarr;
X	register int i,j,inum;
X	extern void qsort();
X
X	BlobTests++;
X
X	ilist = blob->ilist;
X	iarr = blob->iarr;
X
X/*
X * The first step in calculating the Ray/Blob intersection is to 
X * divide the Ray into intervals such that only a fixed set of 
X * Metaballs contribute to the density function along that interval. 
X * This is done by finding the set of intersections between the Ray 
X * and each Metaball's Sphere/Region of influence, which has a 
X * radius  ri  and is centered at (xi,yi,zi).
X * Intersection information is kept track of in the MetaInt 
X * structure and consists of:
X *
X *   type    indicates whether this intersection is the start(R_START)
X *           of a Region or the end(R_END) of one. 
X *   pnt     the Metaball of this intersection
X *   bound   the distance from Ray origin to this intersection
X *
X * This list is then sorted by  bound  and used later to find the Ray/Blob 
X * intersection.
X */
X
X	inum = 0;
X	for(i=0; i < blob->num; i++)
X	{
X		register Float xadj, yadj, zadj;
X		register Float b, t, rs;
X		register Float dmin,dmax;
X
X		rs   = blob->list[i].rs;
X		xadj = blob->list[i].x - ray->pos.x;
X		yadj = blob->list[i].y - ray->pos.y;
X		zadj = blob->list[i].z - ray->pos.z;
X
X	/*
X	 * Ray/Sphere of Influence intersection
X	 */
X		b = xadj * ray->dir.x + yadj * ray->dir.y + zadj * ray->dir.z;
X		t = b * b - xadj * xadj - yadj * yadj - zadj * zadj + rs;
X
X	/* 
X	 * don't except imaginary or single roots. A single root is a ray
X	 * tangent to the Metaball's Sphere/Region. The Metaball's
X	 * contribution to the overall density function at this point is
X	 * zero anyway.
X	 */
X		if (t > 0.0) /* we have two roots */
X		{
X			t = sqrt(t);
X			dmin = b - t;
X	/* 
X	 * only interested in stuff in front of ray origin 
X	 */
X			if (dmin < mindist) dmin = mindist;
X			dmax = b + t;
X			if (dmax > dmin) /* we have a valid Region */
X			{
X		/*
X	         * Initialize min/start and max/end Intersections Structures
X	         * for this Metaball
X	         */
X				ilist[inum].type = R_START;
X				ilist[inum].pnt = i;
X				ilist[inum].bound = dmin;
X				for(j=0;j<5;j++) ilist[inum].c[j] = 0.0;
X				iarr[inum] = &(ilist[inum]);
X				inum++;
X
X				ilist[inum].type = R_END;
X				ilist[inum].pnt = i;
X				ilist[inum].bound = dmax;
X				for(j=0;j<5;j++) ilist[inum].c[j] = 0.0;
X				iarr[inum] = &(ilist[inum]);
X				inum++;
X			} /* End of valid Region */
X		} /* End of two roots */
X	} /* End of loop through metaballs */
X
X	/*
X	 * If there are no Ray/Metaball intersections there will 
X	 * not be a Ray/Blob intersection. Exit now.
X	 */
X	if (inum == 0)
X	{
X		return FALSE;
X	}
X
X	/* 
X	 * Sort Intersection list. No sense using qsort if there's only
X	 * two intersections.
X	 *
X	 * Note: we actually aren't sorting the Intersection structures, but
X	 * an array of pointers to the Intersection structures. 
X	 * This is faster than sorting the Intersection structures themselves.
X	 */
X	if (inum==2)
X	{
X		MetaInt *t;
X		if (ilist[0].bound > ilist[1].bound)
X		{
X			t = iarr[0];
X			iarr[0] = iarr[1];
X			iarr[1] = t;
X		}
X	}
X	else qsort((voidstar)(iarr), (unsigned)inum, sizeof(MetaInt *),
X			MetaCompare);
X
X/*
X* Finding the Ray/Blob Intersection
X* 
X* The non-zero part of the density function for each Metaball is
X* 
X*   di(Ri) = c4i * Ri^4  +  c2i * Ri^2  +  c0i 
X*
X* To find find the Ray/Blob intersection for one metaball
X* substitute for distance   
X*
X*     Ri^2 = (x-xi)^2 + (y-yi)^2 + (z-zi)^2
X* 
X* and then substitute the Ray equations:
X*   
X*     x  = x0 + x1 * t
X*     y  = y0 + y1 * t
X*     z  = z0 + z1 * t
X*
X* to get a big mess :^). Actually, it's a Quartic in t and it's fully 
X* listed farther down. Here's a short version:
X*
X*   c[4] * t^4  +  c[3] * t^3  +  c[2] * t^2  +  c[1] * t  +  c[0]  =  T
X*
X* Don't forget that the Blob is defined by the density being equal to 
X* the threshold T.
X* To calculate the intersection of a Ray and two or more Metaballs, 
X* the coefficients are calculated for each Metaball and then added 
X* together. We can do this since we're working with polynomials.
X* The points of intersection are the roots of the resultant equation.
X*
X* The algorithm loops through the intersection list, calculating
X* the coefficients if an intersection is the start of a Region and 
X* adding them to all intersections in that region.
X* When it detects a valid interval, it calculates the 
X* roots from the starting intersection's coefficients and if any of 
X* the roots are in the interval, the smallest one is returned.
X*
X*/
X
X	{
X		register Float *tmpc;
X		MetaInt *strt,*tmp;
X		register int istrt,itmp;
X		register int num,exitflag,inside;
X
X	/*
X	 * Start the algorithm outside the first interval
X	 */
X		inside = 0;
X		istrt = 0; 
X		strt = iarr[istrt];
X		if (strt->type != R_START)
X			RLerror(RL_WARN,"MetaInt sanity check FAILED!\n");
X
X		/*
X		 * Loop through intersection. If a root is found the code
X		 * will return at that point.
X		 */
X		while( istrt < inum )
X		{
X			/* 
X			 * Check for multiple intersections  at the same point.
X			 * This is also where the coefficients are calculated
X			 * and spread throughout that Metaball's sphere of
X			 * influence.
X			 */
X			do
X			{
X				/* find out which metaball */
X				i = strt->pnt;
X				/* only at starting boundaries do this */
X				if (strt->type == R_START)
X				{
X					register MetaVector *ml;
X					register Float a1,a0;
X					register Float xd,yd,zd;
X					register Float c4,c2,c0;
X
X					/* we're inside */
X					inside++;
X
X	/*  As promised, the full equations
X	 *
X	 *   c[4] = c4*a2*a2; 
X	 *   c[3] = 4.0*c4*a1*a2;
X	 *   c[2] = 4.0*c4*a1*a1 + 2.0*c4*a2*a0 + c2*a2;
X	 *   c[1] = 4.0*c4*a1*a0 + 2.0*c2*a1;
X	 *   c[0] = c4*a0*a0 + c2*a0 + c0;
X	 *
X	 * where
X	 *        a2 = (x1*x1 + y1*y1 + z1*z1) = 1.0 because the ray 
X	 *                                           is normalized
X	 *        a1 = (xd*x1 + yd*y1 + zd*z1)
X	 *        a0 = (xd*xd + yd*yd + zd*zd)
X	 *        xd = (x0 - xi)
X	 *        yd = (y0 - yi)
X	 *        zd = (z0 - zi)
X	 *        (xi,yi,zi) is center of Metaball
X	 *        (x0,y0,z0) is Ray origin
X	 *        (x1,y1,z1) is normalized Ray direction
X	 *        c4,c2,c0   are the coefficients for the
X	 *                       Metaball's density function
X	 *
X	 */
X
X					/*
X					 * Some compilers would recalculate
X					 * this each time its used.
X					 */
X					ml = &(blob->list[i]);
X
X					xd = ray->pos.x - ml->x;
X					yd = ray->pos.y - ml->y;
X					zd = ray->pos.z - ml->z;
X					a1 = (xd * ray->dir.x + yd * ray->dir.y
X						+ zd * ray->dir.z);
X					a0 = (xd * xd + yd * yd + zd * zd);
X
X					c0 = ml->c0;
X					c2 = ml->c2;
X					c4 = ml->c4;
X
X					c[4] = c4;
X					c[3] = 4.0*c4*a1;
X					c[2] = 2.0*c4*(2.0*a1*a1 + a0) + c2;
X					c[1] = 2.0*a1*(2.0*c4*a0 + c2);
X					c[0] = c4*a0*a0 + c2*a0 + c0;
X
X		/* 
X		 * Since we've gone through the trouble of calculating the 
X		 * coefficients, put them where they'll be used.
X		 * Starting at the current intersection(It's also the start of
X		 * a region), add the coefficients to each intersection until
X		 * we reach the end of this region.
X		 */
X					itmp = istrt; 
X					tmp = strt;
X					while( (tmp->pnt != i) ||
X					       (tmp->type != R_END) )
X					{
X						tmpc = tmp->c;
X						for(j=0;j<5;j++)
X							*tmpc++ += c[j];
X						itmp++; 
X						tmp = iarr[itmp];
X					}
X
X				} /* End of start of a Region */ 
X				/* 
X				 * Since the intersection wasn't the start
X				 * of a region, it must the end of one.
X				 */
X				else inside--;
X
X		/*
X		 * If we are inside a region(or regions) and the next
X		 * intersection is not at the same place, then we have
X		 * the start of an interval. Set the exitflag.
X		 */
X				if ((inside > 0)
X				    && (strt->bound != iarr[istrt+1]->bound) )
X					exitflag = 1;
X				else 
X				/* move to next intersection */
X				{
X					istrt++; 
X					strt = iarr[istrt];
X					exitflag = 0;
X				}
X		/* 
X		 * Check to see if we're at the end. If so then exit with
X		 * no intersection found.
X		 */
X				if (istrt >= inum)
X				{
X					return FALSE;
X				}
X			} while(!exitflag);
X				/* End of Search for valid interval */
X
X		/*
X		 * Find Roots along this interval
X		 */
X
X			/* get coefficients from Intersection structure */
X			tmpc = strt->c;
X			for(j=0;j<5;j++) c[j] = *tmpc++;
X
X			/* Don't forget to put in threshold */
X			c[0] -= blob->T;
X
X			/* Use Graphic Gem's Quartic Root Finder */
X			num = SolveQuartic(c,s);
X
X		/*
X		 * If there are roots, search for roots within the interval.
X		 */
X			dist = 0.0;
X			if (num>0)
X			{
X				for(j=0;j<num;j++)
X				{
X		/* 
X		 * Not sure if EPSILON is truly needed, but it might cause
X		 * small cracks between intervals in some cases. In any case 
X		 * I don't believe it hurts.
X		 */
X					if ((s[j] >= (strt->bound - EPSILON))
X					    && (s[j] <= (iarr[istrt+1]->bound +
X					        EPSILON) ) )
X					{
X						if (dist == 0.0)
X						/* first valid root */
X							dist = s[j];
X						else if (s[j] < dist)
X	 					/* else only if smaller */
X							dist = s[j];
X					}
X				}
X				/*
X				 * Found a valid root 
X				 */
X				if (dist > mindist && dist < *maxdist)
X				{
X					*maxdist = dist;
X					BlobHits++;
X					return TRUE;
X					/* Yeah! Return valid root */
X				}
X			}
X
X		/* 
X		 * Move to next intersection
X		 */
X			istrt++;
X			strt = iarr[istrt];
X
X		} /* End of loop through Intersection List */
X	} /* End of Solving for Ray/Blob Intersection */
X
X	/* 
X	 * return negative
X	 */
X	return FALSE;
X}
X
X
X/***********************************************
X * Find the Normal of a Blob at a given point
X *
X * The normal of a surface at a point (x0,y0,z0) is the gradient of that
X * surface at that point. The gradient is the vector 
X *
X *            Fx(x0,y0,z0) , Fy(x0,y0,z0) , Fz(x0,y0,z0)
X *
X * where Fx(),Fy(),Fz() are the partial dervitives of F() with respect
X * to x,y and z respectively. Since the surface of a Blob is made up
X * of the sum of one or more polynomials, the normal of a Blob at a point
X * is the sum of the gradients of the individual polynomials at that point.
X * The individual polynomials in this case are di(Ri) i = 0,...,num .
X *
X * To find the gradient of the contributing polynomials
X * take di(Ri) and substitute U = Ri^2;
X *
X *      di(U)    = c4i * U^2  +  c2i * U  + c0i
X *
X *      dix(U)   = 2 * c4i * Ux * U  +  c2i * Ux
X *
X *        U      = (x-xi)^2 + (y-yi)^2 + (z-zi)^2 
X *
X *        Ux     = 2 * (x-xi) 
X *
X * Rearranging we get
X *  
X *    dix(x0,y0,z0) = 4 * c4i * xdi * disti + 2 * c2i * xdi
X *    diy(x0,y0,z0) = 4 * c4i * ydi * disti + 2 * c2i * ydi
X *    diz(x0,y0,z0) = 4 * c4i * zdi * disti + 2 * c2i * zdi
X * 
X * where
X *         xdi       =   x0-xi
X *         ydi       =   y0-yi
X *         zdi       =   y0-yi
X *        disti      =   xdi*xdi + ydi*ydi + zdi*zdi   
X *       (xi,yi,zi)  is  the center of Metaball i
X *       c4i,c2i,c0i are the coefficients of Metaball i's density function
X */
Xint
XBlobNormal(blob, pos, nrm, gnrm)
XBlob *blob;
XVector *pos, *nrm, *gnrm;
X{
X	register int i;
X
X	/*  
X	 * Initialize normals to zero 
X	 */
X	nrm->x = nrm->y = nrm->z = 0.0;
X	/*
X	 * Loop through Metaballs. If the point is within a Metaball's
X	 * Sphere of influence, calculate the gradient and add it to the
X	 * normals
X	 */
X	for(i=0;i < blob->num; i++)
X	{
X		register MetaVector *sl;
X		register Float dist,xd,yd,zd;
X
X		sl = &(blob->list[i]);
X		xd = pos->x - sl->x;
X		yd = pos->y - sl->y;
X		zd = pos->z - sl->z;
X
X		dist  = xd*xd + yd*yd + zd*zd;
X		if (dist <= sl->rs )
X		{
X			register Float temp;
X
X			/* temp is negative so normal points out of blob */
X			temp = -2.0 * (2.0 * sl->c4 * dist  +  sl->c2);
X			nrm->x += xd * temp;
X			nrm->y += yd * temp;
X			nrm->z += zd * temp;
X
X		}
X	}
X	(void)VecNormalize(nrm);
X	*gnrm = *nrm;
X	return FALSE;
X}
X
X
X/*****************************************************************************
X * Calculate the extent of the Blob
X */
Xvoid
XBlobBounds(blob, bounds)
XBlob *blob;
XFloat bounds[2][3];
X{
X	Float r;
X	Float minx,miny,minz,maxx,maxy,maxz;
X	int i;
X
X	/*
X	 * Loop through Metaballs to find the minimun and maximum in each
X	 * direction.
X	 */
X	for( i=0;  i< blob->num;  i++)
X	{
X		register Float d;
X
X		r = sqrt(blob->list[i].rs);
X		if (i==0)
X		{
X			minx = blob->list[i].x - r;
X			miny = blob->list[i].y - r;
X			minz = blob->list[i].z - r;
X			maxx = blob->list[i].x + r;
X			maxy = blob->list[i].y + r;
X			maxz = blob->list[i].z + r;
X		}
X		else
X		{
X			d = blob->list[i].x - r; 
X			if (d<minx) minx = d;
X			d = blob->list[i].y - r; 
X			if (d<miny) miny = d;
X			d = blob->list[i].z - r; 
X			if (d<minz) minz = d;
X			d = blob->list[i].x + r; 
X			if (d>maxx) maxx = d;
X			d = blob->list[i].y + r; 
X			if (d>maxy) maxy = d;
X			d = blob->list[i].z + r; 
X			if (d>maxz) maxz = d;
X		}
X	}
X	bounds[LOW][X]  = minx;
X	bounds[HIGH][X] = maxx;
X	bounds[LOW][Y]  = miny;
X	bounds[HIGH][Y] = maxy;
X	bounds[LOW][Z]  = minz;
X	bounds[HIGH][Z] = maxz;
X}
X
Xchar *
XBlobName()
X{
X	return blobName;
X}
X
Xvoid
XBlobStats(tests, hits)
Xunsigned long *tests, *hits;
X{
X	*tests = BlobTests;
X	*hits = BlobHits;
X}
X
Xvoid
XBlobMethodRegister(meth)
XUserMethodType meth;
X{
X	if (iBlobMethods)
X		iBlobMethods->user = meth;
X}
END_OF_FILE
if test 18207 -ne `wc -c <'libray/libobj/blob.c'`; then
    echo shar: \"'libray/libobj/blob.c'\" unpacked with wrong size!
fi
# end of 'libray/libobj/blob.c'
fi
if test -f 'rayview/glmethods.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'rayview/glmethods.c'\"
else
echo shar: Extracting \"'rayview/glmethods.c'\" \(17863 characters\)
sed "s/^X//" >'rayview/glmethods.c' <<'END_OF_FILE'
X/*
X * glmethods.c
X *
X * Support routines for SGI/RS6000 machines.
X *
X * Copyright (C) 1989, 1991, Craig E. Kolb, Allan Snyder
X *
X * This software may be freely copied, modified, and redistributed
X * provided that this copyright notice is preserved on all copies.
X *
X * You may not distribute this software, in whole or in part, as part of
X * any commercial product without the express consent of the authors.
X * 
X * There is no warranty or other guarantee of fitness of this software
X * for any purpose.  It is provided solely "as is".
X *
X * $Id: glmethods.c,v 4.0 91/07/17 17:38:43 kolb Exp Locker: kolb $
X *
X * $Log:	glmethods.c,v $
X * Revision 4.0  91/07/17  17:38:43  kolb
X * Initial version
X * 
X */
X#include <gl.h>
X#include <device.h>
X#include "rayshade.h" 
X#include "options.h"
X#include "viewing.h"
X
X#include "libsurf/surface.h"
X
X#include "libobj/box.h"
X#include "libobj/cone.h"
X#include "libobj/csg.h"
X#include "libobj/cylinder.h"
X#include "libobj/disc.h"
X#include "libobj/grid.h"
X#include "libobj/instance.h"
X#include "libobj/list.h"
X#include "libobj/plane.h"
X#include "libobj/poly.h"
X#include "libobj/sphere.h"
X#include "libobj/triangle.h"
X
X#include "liblight/light.h"
X#include "liblight/extended.h"
X#include "liblight/infinite.h"
X#include "liblight/point.h"
X#include "liblight/spot.h"
X
X#define CIRCLE_SAMPLES	20	/* # of samples around circle (cone/cyl/etc) */
X#define DEFAULT_NEAR	.2	/* default value for near clipping plane */
X#define DEFAULT_FAR	350.	/* default far clipping plane */
X#define PLANE_RAD	450.	/* radius of disc used to represent planes */
X
X#ifndef SPHERELIB
X#define SPHERE_LEVEL	3
X#endif
X
X/*
X * Pass a normal stored in a Vector to the geometry pipeline.
X */
X#define GLNormal(w)	(glnrm[0] = (w)->x, glnrm[1] = (w)->y, \
X			 glnrm[2] = (w)->z, n3f(glnrm))
X
Xstatic void	GLBoxDraw(), GLConeDraw(), GLCsgDraw(), GLCylinderDraw(),
X		GLDiscDraw(),
X		GLGridDraw(), GLHfDraw(), GLInstanceDraw(),	
X		GLListDraw(), GLPlaneDraw(), GLPolygonDraw(),
X		GLSphereDraw(), GLTorusDraw(), GLTriangleDraw(),
X		GLBoundsDraw(),
X		GLInfiniteLight(), GLPointLight(), GLSpotLight(),
X		GLExtendedLight();
X
Xstatic short	cursurf = 1,
X		curlight = 1;
X
Xstatic Object	GLBoxObject,
X		GLCylObject,
X		GLSphereObject,
X		GLBoxObjectDefine(),
X		GLCylObjectDefine();
X
Xextern Object	GLSphereObjectDefine();
X
Xstatic int	Doublebuffered;
X
Xstatic float	Ident[4][4] =	{1., 0., 0., 0.,
X				 0., 1., 0., 0.,
X				 0., 0., 1., 0.,
X				 0., 0., 0., 1.},
X		glnrm[3];
X
Xstatic float **CirclePoints;
X
Xstatic void LightDraw(), GeomDraw(), GLMultMatrix(),
X		GLPushSurface(), GLPopSurface(), LightDrawInit(),
X		ObjectInit(), GLDrawFrame(), ScreenDrawInit(), DrawInit(),
X		GLUnitCircle(), GLDisc(), ComputeClippingPlanes();
X
XFloat CurrentTime;
X
Xstatic unsigned long BackPack;	/* Packed background color */
Xstatic long Zinit;		/* maximum Zbuffer value */
X
XMethodsRegister()
X{
X	BoxMethodRegister(GLBoxDraw);
X	ConeMethodRegister(GLConeDraw);
X	CsgMethodRegister(GLCsgDraw);
X	CylinderMethodRegister(GLCylinderDraw);
X	DiscMethodRegister(GLDiscDraw);
X	GridMethodRegister(GLGridDraw);
X	/*HfMethodRegister(GLHfDraw);*/
X	InstanceMethodRegister(GLInstanceDraw);
X	ListMethodRegister(GLListDraw);
X	PlaneMethodRegister(GLPlaneDraw);
X	PolygonMethodRegister(GLPolygonDraw);
X	SphereMethodRegister(GLSphereDraw);
X	/* Have to write sphere routine for RS6000... */
X	/*TorusMethodRegister(GLTorusDraw);*/
X	TriangleMethodRegister(GLTriangleDraw);
X
X	InfiniteMethodRegister(GLInfiniteLight);
X	PointMethodRegister(GLPointLight);
X	ExtendedMethodRegister(GLExtendedLight);
X	SpotMethodRegister(GLSpotLight);
X}
X
Xvoid
XRender()
X{
X	short val;
X	float tmp;
X
X	DrawInit();
X	qdevice(ESCKEY);
X	qdevice(SPACEKEY);
X	qdevice(REDRAW);
X
X	DrawFrames();
X
X	while (TRUE) {
X		switch (qread(&val)) {
X			case ESCKEY:
X				exit(0);
X				break;
X			case REDRAW:
X				reshapeviewport();
X				DrawFrames();
X				break;
X			case SPACEKEY:
X				if (!val)
X					DrawFrames();
X				break;
X		}
X	}
X}
X
XDrawFrames()
X{
X	int i;
X	for (i = 0; i < Options.totalframes; i++)
X		GLDrawFrame(i);
X}
X
Xstatic void
XDrawInit()
X{
X	extern Surface DefaultSurface;
X
X	ScreenDrawInit();
X	ObjectInit();
X	LightDrawInit();
X	/*
X	 * Push the default surface.
X	 */
X	GLPushSurface(&DefaultSurface);
X}
X
Xstatic void
XScreenDrawInit()
X{
X	/*
X	 * Open window, initialize graphics, etc.
X	 */
X#ifdef sgi
X	foreground();
X#endif
X	prefsize(Screen.xsize, Screen.ysize);
X	winopen("rayview");
X	RGBmode();
X	mmode(MVIEWING);
X
X	if (Options.totalframes > 1) {
X		Doublebuffered = TRUE;
X		doublebuffer();
X	}
X
X	gconfig();
X	/*
X	 * Initialize viewing matrix.
X	 */
X	GLViewingInit();
X
X	zbuffer(TRUE);
X
X	BackPack = (unsigned char)(255*Screen.background.r) |
X		((unsigned char)(255*Screen.background.g) << 8) |
X		((unsigned char)(255*Screen.background.b) << 16);
X	Zinit = getgdesc(GD_ZMAX);
X	lsetdepth(0, Zinit);
X}
X
X/*
X * Draw the World object.
X */
Xstatic void
XGLDrawFrame(frame)
Xint frame;
X{
X	extern Geom *World;
X
X	RSStartFrame(frame);
X	CurrentTime = Options.framestart;
X	TimeSet(CurrentTime);
X
X	czclear(BackPack, Zinit);	
X
X	/*
X	 * Draw the World object
X	 */
X	GeomDraw(World);
X	if (Doublebuffered)
X		swapbuffers();
X}
X
XGLViewingInit()
X{
X	float near, far, aspect, dist, T[4][4];
X	short ang;
X	extern Geom *World;
X
X	ang = Camera.vfov * 10 + 0.5;
X	aspect = Camera.hfov / Camera.vfov;
X
X	T[0][0]=Screen.scrni.x; T[0][1]=Screen.scrnj.x; T[0][2]= -Camera.dir.x;
X	T[1][0]=Screen.scrni.y; T[1][1]=Screen.scrnj.y; T[1][2]= -Camera.dir.y;
X	T[2][0]=Screen.scrni.z; T[2][1]=Screen.scrnj.z; T[2][2]= -Camera.dir.z;
X
X	T[0][3] = T[1][3] = T[2][3] = 0.;
X
X	T[3][0] = -dotp(&Camera.lookp, &Screen.scrni);
X	T[3][1] = -dotp(&Camera.lookp, &Screen.scrnj);
X	T[3][2] = dotp(&Camera.lookp, &Camera.dir);
X	T[3][3] = 1.;
X
X	ComputeClippingPlanes(&near, &far, World->bounds);
X
X	loadmatrix(Ident);
X	perspective(ang, aspect, near, far);
X	polarview((float)Camera.lookdist, 0, 0, 0);
X	multmatrix(T);
X}
X
X
Xstatic void
XObjectInit()
X{
X	CircleDefine();
X	GLBoxObject = GLBoxObjectDefine();
X
X#ifdef SPHERELIB
X	GLSphereObject = GLSphereObjectDefine();
X#else
X	GLSphereObject = GLSphereObjectDefine(SPHERE_LEVEL);
X#endif
X
X	GLCylObject = GLCylObjectDefine();
X}
X
Xstatic void
XGeomDraw(obj)
XGeom *obj;
X{
X	Trans *ct;
X	/*
X	 * If object has a surface associated with it,
X	 * start using it.
X	 */
X	if (obj->surf)
X		GLPushSurface(obj->surf);
X	if (obj->trans) {
X		/*
X		 * Take care of any animated transformations that
X		 * exist.
X		 */
X		if (obj->animtrans && !equal(obj->timenow, CurrentTime)) {
X			TransResolveAssoc(obj->trans);
X			obj->timenow = CurrentTime;
X		}
X		pushmatrix();
X		/*
X	 	 * Apply in reverse order.
X		 */
X		for (ct = obj->transtail; ct; ct = ct->prev)
X			GLMultMatrix(&ct->trans);
X	}
X
X	if (obj->methods->user) {
X		/*
X	 	 * Call proper method
X	 	 */
X		(*obj->methods->user)(obj->obj);
X	} else {
X		/*
X		 * Draw bounding box.
X		 */
X		GLBoundsDraw(obj->bounds);
X	}
X
X	if (obj->surf)
X		GLPopSurface();
X	if (obj->trans)
X		popmatrix();
X}
X
Xstatic float surfprops[] =	{AMBIENT, 0, 0, 0,
X				 DIFFUSE, 0, 0, 0,
X				 SPECULAR, 0, 0, 0,
X				 SHININESS, 0,
X				 LMNULL};
Xstatic float	*amb = &surfprops[1],
X		*diff = &surfprops[5],
X		*spec = &surfprops[9],
X		*shine = &surfprops[13];
X
Xstatic void
XGLPushSurface(surf)
XSurface *surf;
X{
X	static Surface *lastsurf = NULL;
X
X	/*
X	 * Start using the given surface.
X	 * In the case of the use of an "applysurf",
X	 * the same surface will be applied to many
X	 * primitives individually.  By saving the
X	 * pointer to the last surface, we keep from
X	 * lmdef'ing the surface when we don't need to.
X	 */
X	if (surf != lastsurf) {
X		amb[0] = surf->amb.r; amb[1] = surf->amb.g;
X			amb[2] = surf->amb.b;
X		diff[0] = surf->diff.r; diff[1] = surf->diff.g;
X			diff[2] = surf->diff.b;
X		spec[0] = surf->spec.r; spec[1] = surf->spec.g;
X			spec[2] = surf->spec.b;
X		shine[0] = surf->srexp;
X		lmdef(DEFMATERIAL, cursurf, 15, surfprops);
X		lastsurf = surf;
X	}
X	lmbind(MATERIAL, cursurf);
X	cursurf++;
X}
X
Xstatic void
XGLMultMatrix(trans)
XRSMatrix *trans;
X{
X	float newmat[4][4];
X	/*
X	 * Multiply in the given transformation.
X	 */
X	newmat[0][0] = trans->matrix[0][0]; newmat[0][1] = trans->matrix[0][1];
X	newmat[0][2] = trans->matrix[0][2]; newmat[0][3] = 0.;
X	newmat[1][0] = trans->matrix[1][0]; newmat[1][1] = trans->matrix[1][1];
X	newmat[1][2] = trans->matrix[1][2]; newmat[1][3] = 0.;
X	newmat[2][0] = trans->matrix[2][0]; newmat[2][1] = trans->matrix[2][1];
X	newmat[2][2] = trans->matrix[2][2]; newmat[2][3] = 0.;
X	newmat[3][0] = trans->translate.x; newmat[3][1] = trans->translate.y;
X	newmat[3][2] = trans->translate.z ; newmat[3][3] = 1.;
X	multmatrix(newmat);
X}
X
Xstatic void
XGLPopSurface()
X{
X	cursurf--;
X	lmbind(MATERIAL, cursurf-1);
X}
X
Xstatic void
XGLBoundsDraw(bounds)
XFloat bounds[2][3];
X{
X	float sx, sy, sz;
X
X	pushmatrix();
X
X	translate(bounds[LOW][X], bounds[LOW][Y], bounds[LOW][Z]);
X	scale(	bounds[HIGH][X] - bounds[LOW][X],
X		bounds[HIGH][Y] - bounds[LOW][Y],
X		bounds[HIGH][Z] - bounds[LOW][Z]);
X	pushmatrix();
X	callobj(GLBoxObject);
X	popmatrix();
X	popmatrix();
X}
X
Xstatic void
XGLBoxDraw(box)
XBox *box;
X{
X	GLBoundsDraw(box->bounds);
X}
X
X
Xstatic void
XGLConeDraw(cone)
XCone *cone;
X{
X	int i, j;
X	float norm[3], normconst, ZeroVect[3], tmpv[3];
X
X	ZeroVect[0] = ZeroVect[1] = ZeroVect[2] = 0.;
X	/*
X	 * Sides.
X	 * For the normal, assume we are finding the normal at
X	 * (C_P[i].x, C_P[i].y, 1.) at this point, the unnormalized
X	 * normal = (C_P[i].x, C_P[i].y, -tantheta^2).  When we normalize,
X	 * then length of the vector is simply equal to:
X	 * sqrt(CP[i].x^2 + CP[i].y^2 + tantheta^4)  ==
X	 * sqrt(1. + tantheta^4)
X	 * In our case, tantheta = 1., so...
X	 */
X
X	normconst = 1. / sqrt(2.);
X	norm[2] = -normconst;
X
X	pushmatrix();
X	GLMultMatrix(&cone->trans.trans);
X
X	for (i = 0; i < CIRCLE_SAMPLES; i++) {
X		j = (i + 1) % CIRCLE_SAMPLES;
X		norm[0] = CirclePoints[i][0] * normconst;
X		norm[1] = CirclePoints[i][1] * normconst;
X		n3f(norm);
X		bgnpolygon();
X		if (cone->start_dist > EPSILON) {
X			tmpv[2] = cone->start_dist;
X			tmpv[0] = CirclePoints[i][0] * cone->start_dist;
X			tmpv[1] = CirclePoints[i][1] * cone->start_dist;
X			v3f(tmpv);
X			tmpv[0] = CirclePoints[j][0] * cone->start_dist;
X			tmpv[1] = CirclePoints[j][1] * cone->start_dist;
X			v3f(tmpv);
X		} else
X			v3f(ZeroVect);
X		tmpv[2] = 1.;
X		tmpv[0] = CirclePoints[j][0];
X		tmpv[1] = CirclePoints[j][1];
X		v3f(tmpv);
X		tmpv[0] = CirclePoints[i][0];
X		tmpv[1] = CirclePoints[i][1];
X		v3f(tmpv);
X		endpolygon();
X	}
X	popmatrix();
X}
X
Xstatic void
XGLCsgDraw(csg)
XCsg *csg;
X{
X	/*
X	 * Punt by drawing both objects.
X	 */
X	GeomDraw(csg->obj1);
X	GeomDraw(csg->obj2);
X}
X
Xstatic void
XGLCylinderDraw(cyl)
XCylinder *cyl;
X{
X	pushmatrix();
X	GLMultMatrix(&cyl->trans.trans);
X	callobj(GLCylObject);
X	popmatrix();
X}
X
Xstatic void
XGLDiscDraw(disc)
XDisc *disc;
X{
X	GLDisc((Float)sqrt(disc->radius), &disc->pos, &disc->norm);
X}
X
Xstatic void
XGLDisc(rad, pos, norm)
XFloat rad;
XVector *pos, *norm;
X{
X	RSMatrix m, tmp;
X	Vector atmp;
X
X	/*
X	 * This is kinda disgusting...
X	 */
X	ScaleMatrix(rad, rad, 1., &m);
X	if (fabs(norm->z) == 1.) {
X		atmp.x = 1.;
X		atmp.y = atmp.z = 0.;
X	} else {
X		atmp.x = norm->y;
X		atmp.y = -norm->x;
X		atmp.z= 0.;
X	}
X
X	RotationMatrix(atmp.x, atmp.y, atmp.z, -acos(norm->z), &tmp);
X	MatrixMult(&m, &tmp, &m);
X	TranslationMatrix(pos->x, pos->y, pos->z, &tmp);
X	MatrixMult(&m, &tmp, &m);
X	pushmatrix();
X	GLMultMatrix(&m);
X	/*
X	 * Draw unit circle.
X	 */
X	GLUnitCircle(0., TRUE);
X	popmatrix();
X}
X
Xstatic void
XGLGridDraw(grid)
XGrid *grid;
X{
X	Geom *ltmp;
X
X	for (ltmp = grid->unbounded; ltmp; ltmp = ltmp->next)
X		GeomDraw(ltmp);
X	for (ltmp = grid->objects; ltmp; ltmp = ltmp->next)
X		GeomDraw(ltmp);
X}
X
Xstatic void
XGLHfDraw(){}
X
Xstatic void
XGLInstanceDraw(inst)
XInstance *inst;
X{
X	GeomDraw(inst->obj);
X}
X
Xstatic void
XGLListDraw(list)
XList *list;
X{
X	Geom *ltmp;
X
X	for (ltmp = list->unbounded; ltmp; ltmp = ltmp->next)
X		GeomDraw(ltmp);
X	for (ltmp = list->list; ltmp; ltmp = ltmp->next)
X		GeomDraw(ltmp);
X}
X
Xstatic void
XGLPlaneDraw(plane)
XPlane *plane;
X{
X	/*
X	 * Draw a really big disc.
X	 */
X	GLDisc((Float)PLANE_RAD, &plane->pos, &plane->norm);
X}
X
Xstatic void
XGLPolygonDraw(poly)
Xregister Polygon *poly;
X{
X	register int i;
X
X	GLNormal(&poly->norm);
X
X	bgnpolygon();
X	for (i = 0; i < poly->npoints; i++)
X		v3d(&poly->points[i]);
X	endpolygon();
X}
X
Xstatic void
XGLSphereDraw(sph)
XSphere *sph;
X{
X	pushmatrix();
X	translate(sph->x, sph->y, sph->z);
X	scale(sph->r, sph->r, sph->r);
X	callobj(GLSphereObject);
X	popmatrix();
X}
X
Xstatic void
XGLTorusDraw(){}
X
Xstatic void
XGLTriangleDraw(tri)
Xregister Triangle *tri;
X{
X	/*
X	 * If Float is a double, use v3d,
X	 * otherwise use v3f.
X	 */
X	if (tri->type != PHONGTRI) {
X		GLNormal(&tri->nrm);	
X		bgnpolygon();
X		v3d(&tri->p[0]); v3d(&tri->p[1]); v3d(&tri->p[2]);
X		endpolygon();
X	} else {
X		bgnpolygon();
X		GLNormal(&tri->vnorm[0]);
X		v3d(&tri->p[0]);
X		GLNormal(&tri->vnorm[1]);
X		v3d(&tri->p[1]);
X		GLNormal(&tri->vnorm[2]);
X		v3d(&tri->p[2]);
X		endpolygon();
X	}
X}
X
Xfloat lightprops[] =	{POSITION, 0., 0., 0., 0.,
X			 LCOLOR, 0, 0, 0,
X			 LMNULL};
X
Xfloat 	*lpos = &lightprops[1],
X	*lcolor = &lightprops[6];
X
Xfloat lmodel[] =	{AMBIENT, 1., 1., 1.,
X			 ATTENUATION, 1., 0.,
X			 LOCALVIEWER, 0.,
X#ifdef sgi
X			 ATTENUATION2, 0.,
X			 TWOSIDE, 1.,
X#endif
X			 LMNULL};
X
Xstatic void
XLightDrawInit()
X{
X	Light *light;
X	extern Light *Lights;
X
X	for (light = Lights; light; light = light->next)
X		LightDraw(light);
X
X	switch (curlight-1) {
X		case 8:
X			lmbind(LIGHT7, 8);
X		case 7:
X			lmbind(LIGHT6, 7);
X		case 6:
X			lmbind(LIGHT5, 6);
X		case 5:
X			lmbind(LIGHT4, 5);
X		case 4:
X			lmbind(LIGHT3, 4);
X		case 3:
X			lmbind(LIGHT2, 3);
X		case 2:
X			lmbind(LIGHT1, 2);
X		case 1:
X			lmbind(LIGHT0, 1);
X	}
X
X	lmodel[1] = Options.ambient.r;
X	lmodel[2] = Options.ambient.g;
X	lmodel[3] = Options.ambient.b;
X
X#ifdef sgi
X	lmdef(DEFLMODEL, 1, 14, lmodel);
X#else
X	lmdef(DEFLMODEL, 1, 10, lmodel);
X#endif
X	lmbind(LMODEL, 1);
X}
X
Xstatic void
XLightDraw(light)
XLight *light;
X{
X	if (!light || !light->methods->user)
X		return;
X	lcolor[0] = light->color.r;
X	lcolor[1] = light->color.g;
X	lcolor[2] = light->color.b;
X	(*light->methods->user)(light->light);
X}
X
Xstatic void
XGLExtendedLight(ext)
XExtended *ext;
X{
X	lpos[0] = ext->pos.x;
X	lpos[1] = ext->pos.y;
X	lpos[2] = ext->pos.z;
X	lpos[3] = 1.;
X	lmdef(DEFLIGHT, curlight++, 10, lightprops);
X}
X
X
Xstatic void
XGLInfiniteLight(inf)
XInfinite *inf;
X{
X	lpos[0] = inf->dir.x;
X	lpos[1] = inf->dir.y;
X	lpos[2] = inf->dir.z;
X	lpos[3] = 0.;
X	lmdef(DEFLIGHT, curlight++, 10, lightprops);
X}
X
Xstatic void
XGLPointLight(pt)
XPointlight *pt;
X{
X	lpos[0] = pt->pos.x;
X	lpos[1] = pt->pos.y;
X	lpos[2] = pt->pos.z;
X	lpos[3] = 1.;
X	lmdef(DEFLIGHT, curlight++, 10, lightprops);
X}
X
Xstatic void
XGLSpotLight(spot)
XSpotlight *spot;
X{
X}
X
Xstatic float boxfaces[6][4][3] = {
X{ 	{1., 1., 1.},
X	{0., 1., 1.},
X	{0., 0., 1.},
X	{1., 0., 1.}	},
X{	{1., 0., 1.},
X	{0., 0., 1.},
X	{0., 0., 0.},
X	{1., 0., 0.}	},
X{	{1., 0., 0.},
X	{0., 0., 0.},
X	{0., 1., 0.},
X	{1., 1., 0.}	},
X{	{0., 1., 0.},
X	{0., 1., 1.},
X	{1., 1., 1.},
X	{1., 1., 0.}	},
X{	{1., 1., 1.},
X	{1., 0., 1.},
X	{1., 0., 0.},
X	{1., 1., 0.}	},
X{	{0., 0., 1.},
X	{0., 1., 1.},
X	{0., 1., 0.},
X	{0., 0., 0.}}};
X
Xfloat boxnorms[6][3] = {
X	{0, 0, 1},
X	{0, -1, 0},
X	{0, 0, -1},
X	{0, 1, 0},
X	{1, 0, 0,},
X	{-1, 0, 0}};
X
X/*
X * Define a unit cube centered at (0.5, 0.5, 0.5)
X */
Xstatic Object
XGLBoxObjectDefine()
X{
X	int i, box;
X
X	makeobj(box = genobj());
X	for (i = 0; i < 6; i++)	{
X		n3f(boxnorms[i]);
X		polf(4, boxfaces[i]);
X	}
X	closeobj();
X	return box;
X}
X
Xstatic void
XGLUnitCircle(z, up)
Xfloat z;
Xint up;
X{
X	int i;
X	float norm[3];
X
X	norm[0] = norm[1] = 0.;
X	bgnpolygon();
X
X	if (up) {
X		norm[2] = 1.;
X		n3f(norm);	
X		for (i = 0; i < CIRCLE_SAMPLES; i++) {
X			CirclePoints[i][2] = z;
X			v3f(CirclePoints[i]);
X		}
X	} else {
X		norm[2] = -1.;
X		n3f(norm);	
X		for (i = CIRCLE_SAMPLES -1; i; i--) {
X			CirclePoints[i][2] = z;
X			v3f(CirclePoints[i]);
X		}
X	}
X
X	endpolygon();
X}
X
Xstatic Object
XGLCylObjectDefine()
X{
X	int i, j, cyl;
X	float norm[3];
X
X	makeobj(cyl = genobj());
X	norm[2] = 0;
X	for (i = 0; i < CIRCLE_SAMPLES; i++) {
X		j = (i+1)%CIRCLE_SAMPLES;
X		norm[0] = CirclePoints[i][0];
X		norm[1] = CirclePoints[i][1];
X#ifdef sgi
X		n3f(norm);
X		bgnpolygon();
X		CirclePoints[i][2] = 0;
X		v3f(CirclePoints[i]);
X		CirclePoints[j][2] = 0;
X		v3f(CirclePoints[j]);
X		CirclePoints[j][2] = 1;
X		v3f(CirclePoints[j]);
X		CirclePoints[i][2] = 1;
X		v3f(CirclePoints[i]);
X		endpolygon();
X#else
X		n3f(norm);
X		pmv(CirclePoints[i][0], CirclePoints[i][1], 0.);
X		pdr(CirclePoints[j][0], CirclePoints[j][1], 0.);
X		pdr(CirclePoints[j][0], CirclePoints[j][1], 1.);
X		pdr(CirclePoints[i][0], CirclePoints[i][1], 1.);
X		pclos();
X#endif
X	}
X	closeobj();
X	return cyl;
X}
X
X/*
X * Fill CirclePoints[t] with X and Y values cooresponding to points
X * along the unit circle.  The size of CirclePoints is equal to
X * CIRCLE_SAMPLES.
X */
XCircleDefine()
X{
X	int i;
X	double theta, dtheta;
X
X	dtheta = 2.*M_PI / (double)CIRCLE_SAMPLES;
X	CirclePoints = (float **)malloc(CIRCLE_SAMPLES * sizeof(float *));
X	for (i = 0, theta = 0; i < CIRCLE_SAMPLES; i++, theta += dtheta) {
X		CirclePoints[i] = (float *)malloc(3 * sizeof(float));
X		CirclePoints[i][0] = cos(theta);
X		CirclePoints[i][1] = sin(theta);
X		/* Z is left unset. */	
X	}
X}
X
X/*
X * Given world bounds, determine near and far
X * clipping planes.  Only problem occurs when there are
X * unbbounded objects (planes)...
X */
Xstatic void
XComputeClippingPlanes(near, far, bounds)
Xfloat *near, *far;
XFloat bounds[2][3];
X{
X	Vector e, c;
X	double rad, d;
X
X
X	/*
X	 * Compute 'radius' of scene.
X	 */
X	rad = max(max((bounds[HIGH][X] - bounds[LOW][X])*0.5,
X		      (bounds[HIGH][Y] - bounds[LOW][Y])*0.5),
X		      (bounds[HIGH][Z] - bounds[LOW][Z])*0.5);
X
X	c.x = (bounds[LOW][X] + bounds[HIGH][X])*0.5;
X	c.y = (bounds[LOW][Y] + bounds[HIGH][Y])*0.5;
X	c.z = (bounds[LOW][Z] + bounds[HIGH][Z])*0.5;
X
X	/*
X	 * Compute position of eye relative to the center of the sphere.
X	 */
X	VecSub(Camera.pos, c, &e);
X	d = VecNormalize(&e);
X	if (d <= rad)
X		/*
X		 * Eye is inside sphere.
X		 */
X		*near = DEFAULT_NEAR;
X	else
X		*near = (d - rad)*0.2;
X	*far = (d + rad)*2.;
X}
END_OF_FILE
if test 17863 -ne `wc -c <'rayview/glmethods.c'`; then
    echo shar: \"'rayview/glmethods.c'\" unpacked with wrong size!
fi
# end of 'rayview/glmethods.c'
fi
echo shar: End of archive 17 \(of 19\).
cp /dev/null ark17isdone
MISSING=""
for I in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ; do
    if test ! -f ark${I}isdone ; then
	MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 19 archives.
    rm -f ark[1-9]isdone ark[1-9][0-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0