
# This is a shell archive.  Remove anything before this line,
# then unpack it by saving it in a file and typing "sh file".
#
# Wrapped by 	Eric Haines <erich@juniper> on Mon Dec 17 16:05:03 1990
#
# This archive contains:
#	README		NFF		makefile	patchlevel.h	
#	def.h		lib.h		lib.c		balls.c		
#	gears.c		mount.c		rings.c		teapot.c	
#	tetra.c		tree.c		
#

LANG=""; export LANG
PATH=/bin:/usr/bin:$PATH; export PATH

echo x - README
cat >README <<'@EOF'
STANDARD PROCEDURAL DATABASES, by Eric Haines, 3D/Eye, Inc.

[Created while under contract to Hewlett-Packard FSD and HP Laboratories]
Version 3.0, as of 12/17/90
    address: 3D/Eye, Inc., 2359 N. Triphammer Rd, Ithaca, NY 14850
    email: erich@eye.com, wrath.cs.cornell.edu!eye!erich

This software package is not copyrighted and can be used freely.

History
-------
Versions 1.0 to 1.5 released February to June, 1987 for testing.
Version 2.0 released July, 1987.
Version 2.1 released August, 1987 - corrected info on speed of the HP 320,
    other minor changes to README.
Version 2.2 released November, 1987 - shortened file names to <=12 characters,
    procedure names to <= 32 characters, and ensured that all lines are <= 80
    characters (including return character).
Version 2.3 released March, 1988 - corrected gears.c to avoid interpenetration,
    corrected and added further instructions and global statistics for ray
    tracing to README.
Version 2.4 released May, 1988 - fixed hashing function for mountain.c.
Version 2.5 released October, 1988 - added NFF documentation.
Version 2.6 released August, 1989 - lib_output_cylcone fix (start_norm.w was
    not initialized).
Version 2.7 released July, 1990 - comment correction in lib.c, NFF file
    clarifications.
Version 3.0 released October, 1990 - added teapot.c database, mountain.c
    changed to mount.c, additions to def.h and lib.c, changes to README and
    NFF, added patchlevel.h file.

{These files use tab characters worth 8 spaces}


Introduction
------------

    This software is meant to act as a set of basic test images for ray
tracing algorithms.  The programs generate databases of objects which are
fairly familiar and "standard" to the graphics community, such as the teapot,
a fractal mountain, a tree, a recursively built tetrahedral structure, etc.  I
originally created them for my own testing of ray tracing efficiency schemes.
Since their first release other researchers have used them to test new
algorithms.  In this way, research on algorithmic improvements can be compared
in a more standardized way.  If one researcher ray-traces a car, another a
tree, the question arises, "How many cars to the tree?"  With these databases
we may be comparing oranges and apples ("how many hypercubes to a timeshared
VAX?"), but it's better than comparing oranges and orangutans.  In addition,
this document outlines some statistics that are more meaningful to researchers
than raw timing tests.  Using these statistics along with the same scenes
allows us to compare results in a more meaningful way.

    With the development and release of the Anderson benchmarks for graphics
hardware, the use of the SPD package for hardware testing is somewhat
redundant.  Therefore I have deleted references to testing hidden-surface
algorithms in this version.  However, another interesting use for the SPD has
been noted:  debugging.  By comparing the images and the statistics with the
output of your own ray tracer, you can detect program errors.  For example,
"mount" is useful for checking if refraction rays are generated correctly, and
"balls" can check for the correctness of eye and reflection rays.

    The images for these databases and other information about them can be
found in "A Proposal for Standard Graphics Environments," IEEE Computer
Graphics and Applications, vol. 7, no. 11, November 1987, pp.  3-5.  See
IEEE CG&A, vol. 8, no. 1, January 1988, p. 18 for the correct image of the
tree database (the only difference is that the sky is blue, not orange).  The
teapot database was added later, and consists of a shiny teapot on a shiny
checkerboard.

    The SPD package is available via anonymous FTP from:

	freedom.graphics.cornell.edu [128.84.247.85]
	weedeater.math.yale.edu [130.132.23.17]
	cs.uoregon.edu [128.223.4.13]

among others.  For those without FTP access, write to the netlib automatic
mailer:  research!netlib and netlib@ornl.gov are the sites.  Send a one line
message "send index" for more information, or "send haines from graphics" for
the latest version of the SPD package.


File Structure
--------------

    Eight different procedural database generators are included.  These were
designed to span a fair range of primitives, modeling structures, lighting and
surface parameters, background conditions, and other factors.  A complexity
factor is provided within each program to control the size of the database
generated.

This software package contains the following files:

    README - what you are now reading
    NFF - a description of the Neutral File Format used in the SPD
    patchlevel.h - keeps track of patch level
    makefile - used to make the programs (in HP-UX Unix)
    def.h - some useful "C" definitions
    lib.h - globals and library routine declarations
    lib.c - library of routines
    balls.c - fractal ball object (a.k.a. sphereflake) generator
    gears.c - 3D array of interlocking gears generator
    mount.c - fractal mountain and 4 glass ball generator
    rings.c - pyramid of dodecahedral rings generator
    teapot.c - the famous teapot on a checkerboard generator
    tetra.c - recursive tetrahedra generator
    tree.c - tree generator

    The compiled and linked programs will output a database in ASCII to stdout
containing viewing parameters, lighting conditions, material properties, and
geometric information.  The data format is called the 'neutral file format'
(or NFF) and is described in the `NFF' file.  This format is meant to be
minimal and easy to attach to a user-written filter program which will convert
the output into a file format of your choosing.

    Either of two sets of primitives can be selected for output.  If
OUTPUT_FORMAT is defined as OUTPUT_CURVES, the primitives are spheres, cones,
cylinders, and polygons.  If OUTPUT_FORMAT is set to OUTPUT_PATCHES, the
primitives are output as polygonal patches and polygons (i.e.  all other
primitives are polygonalized).  In this case OUTPUT_RESOLUTION is used to set
the amount of polygonalization of non-polygonal primitives.  In general,
OUTPUT_CURVES is used for ray-trace timing tests, and OUTPUT_PATCHES for
polygon-based algorithm timings.  Note that sphere primitives are not
polygonalized using the simple longitude/latitude tessellation, but rather are
chopped along the six faces of a cube projected onto the sphere.  This
tessellation avoids generation of the cusp points at the sphere's poles and so
eliminates discontinuities due to them.

    The size factor is used to control the overall size of the database.
Default values have been chosen such that the maximum number of primitives
less than 10,000 is output.  One purpose of the size factor is to avoid
limiting the uses of these databases.  Depending on the research being done
and the computing facilities available, a larger or smaller number of
primitives may be desired.  The size factor can also be used to show how an
algorithm's time changes as the complexity increases.

    To generate the databases, simply type the name of the database and direct
the output as desired, e.g. "balls > balls.nff" creates the default sized
database and sends the output to balls.nff.  A new feature in 3.0 is that you
can enter the size factor on the command line, e.g. "balls 2 > balls.nff"
gives a much smaller database of 91 spheres.  See the header of the database C
file (e.g. "balls.c") for how the size factor affects the output.

    Other parameters in the code itself (for example, branching angles for
"tree.c" and the fractal dimension in "mount.c") are included for your own
enjoyment, and so normally should not be changed if the database is used for
timing tests.  Because the hashing function in the original release of
"mount.c" is not very good (at low resolutions there is patterning), there is
a better hashing function provided which can be turned on by defining
NEW_HASH.  Use the old hashing function (i.e.  don't change anything) for
consistency with previously published results.

    Since the SPD package is designed to test efficiency, the actual shading
of the images is mostly irrelevant.  All that matters is that reflective
surfaces spawn reflection rays and transmitters spawn refraction rays.  For
this reason the effect of the other shading parameters is up to the
implementer.  Note that light intensities, ambient components, etc. are not
specified.  These may be set however you prefer.  Feel free to change any
colors you wish (especially the garish colors of `rings').  The thrust of
these databases is the testing of rendering speeds, and so the actual color
should not affect these calculations.  An ambient component should be used for
all objects, so that the time to compute it is included in the ray tracing
statistics.  A simple formula for a relative intensity (i.e. between 0 and
1) for each light and for the ambient component is the following:

	sqrt(# lights) / (# lights * 2).

    If you desire a model for a good scene description language, NFF is not
it.  Instead, you should look at Craig Kolb's ray-tracer language, which is
part of his excellent RayShade ray tracer, available via anonymous FTP from
weedeater.math.yale.edu [130.132.23.17].  Of the public domain ray tracers,
this is the one to beat in speed.  Also included in his distribution is an awk
filter which converts NFF files to his language.

    The programs all attempt to minimize program size (for ease in
distributing) and memory usage, and none allocate any dynamic memory.  As
such, many of the programs are relatively inefficient, regenerating various
data again and again instead of saving intermediate results.  This behavior
was felt to be unimportant, since generating the data is a one-time affair
which normally takes much less time than actually rendering the scene.


Database Analysis
-----------------

    The databases were designed with the idea of diversity in mind.  The
variables considered important were the amount of background visible, the
number of lights, the distribution of sizes of objects, the amount of
reflection and refraction, and the depth complexity (how many objects a ray
from the eye passes through).

balls:  This database is sometimes called "sphereflake", as it is generated
    like a snowflake curve.  This database consists mostly of various sized
    spheres.  It has no eye rays which hit the background, and the three light
    sources cause a large number of shadow rays to be generated.

gears:  This database consists of a set of meshed gears.  Some of the gears
    are transmitters, making this database lengthy to render.  The gear faces
    each have 144 vertices, and thus tests polygon inside/outside test
    efficiency.  Depth complexity is medium.

mount:  The fractal mountain generator is derived from Loren Carpenter's
    method, and the composition with the four glass spheres is inspired by
    Peter Watterberg's work.  Most objects are tiny (i.e. fractal facets),
    but rendering time is dominated by the rendering of the four large
    spheres.  Depth complexity is low, and there is much background area.

rings:  Objects with six pentagonal rings form a pyramid against a background
    polygon.  With a high amount of interreflection and shadowing, this scene
    is fairly lengthy to render.  Depth complexity is also high, with
    all of the objects partially or fully obscured.

teapot:  The famous teapot on a checkerboard.  There are a number of
    variations on the teapot database, i.e.  whether the bottom is included
    (the IEEE CG&A article added a bottom), variations in the control points
    on the lid (which creates thin, almost degenerate triangles), etc.  For
    this database generator, the bottom is created, and the degenerate
    polygonal patches at the center of the lid and bottom are not output.  The
    bottom of the teapot is not flat, interestingly enough.  The resolution of
    the checkerboard shows the resolution of the teapot meshing (with each
    teapot quadrilateral cut into two triangles), e.g. an 8x8 checkerboard
    means that 8x8x2 triangles are generated per patch.  The definitions for
    the 32 Bezier patches are a part of the program.  All objects are
    reflective and there are two light sources.  Depth complexity is low.

tetra:  A recursive tetrahedral pyramid, first visualized by Benoit Mandelbrot
    and Alan Norton.  This scene is dominated by background (around 80%).  With
    the objects not being reflective and there being only one light source,
    this database is particularly quick to render, with various forms of
    coherency being very useful.  Depth complexity is medium, though some
    rays must pass by many triangles for some of the background pixels.

tree:  A tree formed using Aono and Kunii's tree generation method.  With
    seven light sources, the emphasis is on shadow testing.  Shadow caching
    yields little improvement due to the narrow primitives, and many shadow
    rays pass through complex areas without hitting any objects.  There is a
    fair amount of background area.  Depth complexity is low.


                balls      gears       mount       rings
                -----      -----       -----       -----
primitives        SP         P           PS         YSP
total prim.      7382       9345        8196        8401
poly/patches    1417K       9345        8960        874K

lights            3          5           1           3
background        0%         7%         34%          0%
specular         yes        yes         yes         yes
transmitter       no        yes         yes          no

eye hit rays   263169     245086      173125      263169
reflect rays   175095     304643      354769      315236
refract rays        0     207564      354769           0
shadow rays    954368    2246955      412922     1085002


               teapot       tetra        tree
               ------       -----        ----
primitives       TP           P           OSP
total prim.      9264        4096        8191
poly/patches     9264        4096        852K

lights            2           1           7
background       39%         81%         35%
reflector        yes          no          no
transmitter       no          no          no

eye hit rays   161120       49788      169836
reflect rays   225248           0           0
refract rays        0           0           0
shadow rays    407656       46112     1097419


    "primitives" are S=sphere, P=polygon, T=polygonal patch, Y=cylinder,
O=cone, listed from most in database to least.

    "total prim." is the total number of ray-tracing primitives (polygons,
spheres, cylinders and cones) in the scene.  The number of polygons and vectors
generated is a function of the OUTPUT_RESOLUTION.  The default value for this
parameter is 4 for all databases.

    "poly/patches" is the total number of polygons and polygonal patches
generated when using OUTPUT_PATCHES.

    "lights" is simply the number of lights in a scene.  "background" is the
percentage of background color (empty space) seen directly by the eye for the
given view.  "reflector" tells if there are reflective objects in the scene,
and "transmitter" if there are transparent objects.

    "eye hit rays" is the number of rays from the eye which actually hit an
object (i.e. not the background).  513x513 eye rays are assumed to be shot
(i.e. one ray per pixel corner).  "reflect rays" is the total number of rays
generated by reflection off of reflecting and transmitting surfaces.  "refract
rays" is the number of rays generated by transmitting surfaces.  "shadow rays"
is the sum total of rays shot towards the lights.  Note that if a surface is
facing away from a light, or the background is hit, a light ray is not formed.
The numbers given can vary noticeably from a given ray tracer, but should all
be within about 10%.

    "K" means exactly 1000 (not 1024), with number rounded to the nearest K.
All of the above statistics should be approximately the same for all
classical ray tracers.


Testing Procedures
------------------

Below are listed the requirements for testing various algorithms.  These test
conditions should be realizable by most renderers, and are meant to represent
a common mode of operation for each algorithm.  Special features which the
software supports (or standard features which it lacks) should be noted in
your statistics.

    1)  The non-polygon (OUTPUT_CURVES) format should normally be used for
	ray-tracing tests.

    2)  All opaque (non-transmitting) primitives can be considered one-sided
	for rendering purposes.  Only the outside of primitives are visible in
	the scenes.  The only exception to this is the "teapot" database, in
	which the teapot itself should normally be double sided (this is
	necessary because the lid of the teapot does not fit tightly, allowing
	the viewer to see back faces).

    3)  Polygonal patches (which are always triangles in the SPD) should have
	smooth shading, if available.

    4)  Specular highlighting should be performed for surfaces with a
	reflective component.  The simple Phong distribution model is
	sufficient.

    5)  Light sources are positional.  If unavailable, assign the directional
	lights a vector given by the light position and the viewing "lookat"
	position.

    6)  Render at a resolution of 512 x 512 pixels, shooting rays at the
	corners (meaning that 513 x 513 eye rays will be created).  The four
	corner contributions are averaged to arrive at a pixel value.  If
	rendering is done differently, note this fact.  No pixel subdivision
	is performed.

    7)  The maximum tree depth is 5 (where the eye ray is of depth 1).  Beyond
	this depth rays do not have to be spawned.

    8)  All rays hitting only reflective and transmitting objects spawn
	reflection rays, unless the maximum ray depth was reached by the
	spawning ray.  No adaptive tree depth cutoff is allowed; that is, all
	rays must be spawned (adaptive tree depth is a proven time saver and
	is also dependent on the color model used - see Roy Hall's work for
	details).

    9)  All rays hitting transmitting objects spawn refraction rays, unless
	the maximum ray depth was reached or total internal reflection occurs.
	Transmitting rays should be refracted using Snell's Law (i.e. should
	not pass straight through an object).  If total internal reflection
	occurs, a reflection ray should still be generated at this node.  Note
	that all transmitters in the SPD are also reflective, but this is
	not a requirement of the file format itself.

    10) A shadow ray is not generated if the surface normal points away from
	the light.  This is true even on transmitting surfaces.  Note any
	changes from this condition.

    11) Assume no hierarchy is given with the database (for example, color
	change cannot be used to note a clustering).  The ray tracing program
	itself can create its own hierarchy, but this process should be
	automatic.  Note any exceptions to this (e.g. not including the
	background polygon in the efficiency structure, changing the
	background polygon into an infinite plane, etc).  Such changes can be
	critical for the efficiency of some schemes, so an explanation of why
	changes were made is important.

    12) Timing costs should be separated into at least two areas: preprocessing
	and ray-tracing.  Preprocessing includes all time spent initializing,
	reading the database, and creating data structures needed to ray-trace.
	Preprocessing should be all the constant cost operations--those that
	do not change with the resolution of the image.  Ray-tracing is the
	time actually spent tracing the rays (i.e. everything that is not
	preprocessing).

    13) Other timing costs which would be of interest are in a breakdown of
	times spent in preprocessing and during actual ray-tracing.  Examples
	include time spent reading in the data itself, creating a hierarchy,
	octree, or item buffer, and times spent on intersection the various
	primitives and on calculating the shade.

    14) Time-independent statistics on the performance of the algorithm should
	be gathered.  Some examples are the number of ray/object intersection
	tests, the number of ray/object tests which actually hit objects,
	number of octree or grid nodes accessed, and the number of successful
	shadow cache hits.


Timings
-------

Rendering time of the ray-traced test set on an HP-835, optimized (512 x 512):

       Input   Setup  Ray-Tracing  |  Polygon   Sphere    Cyl/Cone    Bounding
            (hr:min:sec)           |   Tests     Tests      Tests    Box Tests
-------------------------------------------------------------------------------
balls   0:14    0:19      24:56    |    822K     6197K        0        51726K
gears   0:56    0:18    1:02:31    |  13703K       0          0       107105K
mount   0:31    0:14      18:49    |   4076K     3978K        0        31106K
rings   0:41    0:33      53:41    |   1045K     5315K     16298K      91591K
teapot  1:02    0:18      25:38    |   7281K       0          0        57050K
tetra   0:15    0:06       3:48    |    965K       0          0         7637K
tree    0:40    0:31      11:36    |    479K      524K      1319K      22002K

Input:  time spent reading in the NFF file.

Setup:  time spent creating the database and any ray tracing efficiency
structures, and cleaning up after ray trace (does not include "Input" time).

Ray-Tracing:  time spent traversing and rendering the pixel grid.


A typical set of ray tracing intersection statistics for the tetra database is:

[these statistics should be the same for all users]
  image size: 512 x 512
  total number of pixels: 262144			[ 512 x 512 ]
  total number of trees generated: 263169		[ 513 x 513 ]
  total number of tree rays generated: 263169		[ no rays spawned ]
  number of eye rays which hit background: 213381	[ 81% - might vary ]
  average number of rays per tree: 1.000000
  average number of rays per pixel: 1.003910
  total number of shadow rays generated: 46111		[ might vary a bit ]

[these tests vary depending on the ray-tracing algorithm used]
Intersector performance
  Bounding box intersections:      7636497  -    26.71 usec/test
  Polygon intersections:            964567  -    36.77 usec/test

Ray generation
  Eye rays generated:              263169 ( 49788 hit - 18.92% )
  Reflection rays generated:            0
  Refraction rays generated:            0
  Shadow rays generated:            46111
      Coherency hits:                3407  -  7.39 % of total
      Fully tested:                 42704


    The ray-tracer which generated these statistics is based on hierarchical
bounding boxes generated using Goldsmith & Salmon's automatic hierarchy method
(see IEEE CG&A May 1987).  It no longer uses an item buffer, hence the higher
number of overall intersection tests from earlier SPD versions.

    One problem worth analyzing when using the SPD for efficiency tests is how
octree, SEADS, and other space dividing algorithms perform when the background
polygon dimensions are changed (thus changing the size of the outermost
enclosing box, which changes the encoding of the environment).  One analysis of
the effect of the background polygon on RayShade can be found in "Ray Tracing
News", volume 3, number 1.


Future Work
-----------

    These databases are not meant to be the ultimate in standards, but are
presented as an attempt at providing representative modeled environments.  A
number of extensions to the file format could be provided someday, along with
new database generators which use them.  The present databases do not contain
polygons with holes, spline patches, polygonal mesh, triangular strip, or
polyhedron data structures.  Modeling matrices are not output, and CSG
combined primitives are not included.  Light sources are particularly
simplistic.  For a richer, more user-friendly scene description language,
take a look at RayShade or RenderMan.

    As far as database geometry is concerned, most scenes have a preponderance
of small primitives.  If you find that these databases do not reflect the type
of environments you render, please write and explain why (or better yet, write
one or more programs that will generate your "typical" environments--maybe it
will get put in the next release).


Acknowledgements
----------------

    I originally heard of the idea of standard scenes from Don Greenberg back
in 1984.  Some time earlier he and Ed Catmull had talked over coming up with
some standard databases for testing algorithmic claims, and to them must go
the credit for the basic concept.  The idea of making small programs which
generate the data came to me after Tim Kay generously sent huge files of his
tree database via email - I felt there had to be a better way.  Adding the
teapot was inspired by the repeated demand on comp.graphics for this database,
certainly the most famous of all.

    Many thanks to the reviewers, listed alphabetically:  Kells Elmquist, Jeff
Goldsmith, Donald Greenberg, David Hook, Craig Kolb, Susan Spach, Rick Speer,
K.R. Subramanian, J. Eric Townsend, Mark VandeWettering, John Wallace, Greg
Ward, and Louise Watson.  Other people who have kindly offered their ideas and
opinions on this project include Brian Barsky, Andrew Glassner, Roy Hall, Chip
Hatfield, Tim Kay, John Recker, Paul Strauss, and Chan Verbeck.  These names
are mentioned mostly as a listing of people interested in this idea.  They do
not necessarily agree (and in some cases strongly disagree) with the validity
of the concept or the choice of databases.

    Your comments and suggestions on these databases are appreciated.  Please
do send me any timing results for software and hardware which you test, or any
publications which use the SPD package.
@EOF

chmod 644 README

echo x - NFF
cat >NFF <<'@EOF'
Neutral File Format,
    by Eric Haines, 3D/Eye Inc, 2359 North Triphammer Road, Ithaca, NY  14850
    email: erich@eye.com, wrath.cs.cornell.edu!eye!erich

Draft #1, 10/3/88
Version 2.7, as of 5/22/90 - added information on hither, light color
Version 3.0, 12/17/90 - minor information changes

The NFF (Neutral File Format) is designed as a minimal scene description
language.  The language was designed in order to test various rendering
algorithms and efficiency schemes.  It is meant to describe the geometry and
basic surface characteristics of objects, the placement of lights, and the
viewing frustum for the eye.  Some additional information is provided for
esthetic reasons (such as the color of the objects, which is not strictly
necessary for testing the efficiency of rendering algorithms).

Future enhancements could include:  circle and torus objects, spline surfaces
with trimming curves, directional lights, characteristics for positional
lights, CSG descriptions, and probably more by the time you read this.
Comments, suggestions, and criticisms are all welcome.

At present the NFF file format is used in conjunction with the SPD (Standard
Procedural Database) software, a package designed to create a variety of
databases for testing rendering schemes.  For more information about SPD see
"A Proposal for Standard Graphics Environments," IEEE Computer Graphics and
Applications, vol. 7, no. 11, November 1987, pp. 3-5.  See IEEE CG&A, vol. 8,
no. 1, January 1988, p. 18 for the correct image of the tree database (the
only difference is that the sky is blue, not orange).

The SPD package is available via anonymous FTP from:

	freedom.graphics.cornell.edu [128.84.247.85]
	weedeater.math.yale.edu [130.132.23.17]
	cs.uoregon.edu [128.223.4.13]

among others.  For those without FTP access, write to the netlib automatic
mailer:  research!netlib and netlib@ornl.gov are the sites.  Send a one line
message "send index" for more information, or "send haines from graphics" for
the latest version of the SPD package.  If you're looking for a good scene
description language, see Craig Kolb's RayShade software (found on weedeater).
NFF is meant for testing efficiency schemes and so has minimal support for
lighting and shading.

By providing a minimal interface, NFF is meant to act as a simple format to
allow the programmer to quickly write filters to move from NFF to the
local file format.  Presently the following entities are supported:
     A simple perspective frustum
     A background color description
     A positional (vs. directional) light source description
     A surface properties description
     Polygon, polygonal patch, cylinder/cone, and sphere descriptions

Files are output as lines of text.  For each entity, the first field defines
its type.  The rest of the line and possibly other lines contain further
information about the entity.  Entities include:

"v"  - viewing vectors and angles
"b"  - background color
"l"  - positional light location
"f"  - object material properties
"c"  - cone or cylinder primitive
"s"  - sphere primitive
"p"  - polygon primitive
"pp" - polygonal patch primitive


These are explained in depth below.

--------

Viewpoint location.  Description:
    "v"
    "from" Fx Fy Fz
    "at" Ax Ay Az
    "up" Ux Uy Uz
    "angle" angle
    "hither" hither
    "resolution" xres yres

Format:

    v
    from %g %g %g
    at %g %g %g
    up %g %g %g
    angle %g
    hither %g
    resolution %d %d

The parameters are:

    From:  the eye location in XYZ.
    At:    a position to be at the center of the image, in XYZ world
	   coordinates.  A.k.a. "lookat".
    Up:    a vector defining which direction is up, as an XYZ vector.
    Angle: in degrees, defined as from the center of top pixel row to
	   bottom pixel row and left column to right column.
    Hither: distance of the hither plane (if any) from the eye.  Mostly
	   needed for hidden surface algorithms.
    Resolution: in pixels, in x and in y.

  Note that no assumptions are made about normalizing the data (e.g. the
  from-at distance does not have to be 1).  Also, vectors are not
  required to be perpendicular to each other.

  For all databases some viewing parameters are always the same:
    Yon is "at infinity."
    Aspect ratio is 1.0.

  A view entity must be defined before any objects are defined (this
  requirement is so that NFF files can be displayed on the fly by hidden
  surface machines).

--------

Background color.  A color is simply RGB with values between 0 and 1:
    "b" R G B

Format:
    b %g %g %g

    If no background color is set, assume RGB = {0,0,0}.

--------

Positional light.  A light is defined by XYZ position.  Description:
    "l" X Y Z [R G B]

Format:
    l %g %g %g [%g %g %g]

    All light entities must be defined before any objects are defined (this
    requirement is so that NFF files can be used by hidden surface machines).
    Lights have a non-zero intensity of no particular value, if not specified
    (i.e. the program can determine a useful intensity as desired); the
    red/green/blue color of the light can optionally be specified.

--------

Fill color and shading parameters.  Description:
    "f" red green blue Kd Ks Shine T index_of_refraction

Format:
    f %g %g %g %g %g %g %g %g

    RGB is in terms of 0.0 to 1.0.

    Kd is the diffuse component, Ks the specular, Shine is the Phong cosine
    power for highlights, T is transmittance (fraction of contribution of the
    transmitting ray).  Usually, 0 <= Kd <= 1 and 0 <= Ks <= 1, though it is
    not required that Kd + Ks == 1.  Note that transmitting objects ( T > 0 )
    are considered to have two sides for algorithms that need these (normally
    objects have one side).

    The fill color is used to color the objects following it until a new color
    is assigned.

--------

Objects:  all objects are considered one-sided, unless the second side is
needed for transmittance calculations (e.g. you cannot throw out the second
intersection of a transparent sphere in ray tracing).

Cylinder or cone.  A cylinder is defined as having a radius and an axis
    defined by two points, which also define the top and bottom edge of the
    cylinder.  A cone is defined similarly, the difference being that the apex
    and base radii are different.  The apex radius is defined as being smaller
    than the base radius.  Note that the surface exists without endcaps.  The
    cone or cylinder description:

    "c"
    base.x base.y base.z base_radius
    apex.x apex.y apex.z apex_radius

Format:
    c
    %g %g %g %g
    %g %g %g %g

    A negative value for both radii means that only the inside of the object is
    visible (objects are normally considered one sided, with the outside
    visible).  Note that the base and apex cannot be coincident for a cylinder
    or cone.  Making them coincident could be used to define endcaps, but none
    of the SPD scenes currently make use of this definition.

--------

Sphere.  A sphere is defined by a radius and center position:
    "s" center.x center.y center.z radius

Format:
    s %g %g %g %g

    If the radius is negative, then only the sphere's inside is visible
    (objects are normally considered one sided, with the outside visible).
    Currently none of the SPD scenes make use of negative radii.

--------

Polygon.  A polygon is defined by a set of vertices.  With these databases,
    a polygon is defined to have all points coplanar.  A polygon has only
    one side, with the order of the vertices being counterclockwise as you
    face the polygon (right-handed coordinate system).  The first two edges
    must form a non-zero convex angle, so that the normal and side visibility
    can be determined by using just the first three vertices.  Description:

    "p" total_vertices
    vert1.x vert1.y vert1.z
    [etc. for total_vertices vertices]

Format:
    p %d
    [ %g %g %g ] <-- for total_vertices vertices

--------

Polygonal patch.  A patch is defined by a set of vertices and their normals.
    With these databases, a patch is defined to have all points coplanar.
    A patch has only one side, with the order of the vertices being
    counterclockwise as you face the patch (right-handed coordinate system).
    The first two edges must form a non-zero convex angle, so that the normal
    and side visibility can be determined.  Description:

    "pp" total_vertices
    vert1.x vert1.y vert1.z norm1.x norm1.y norm1.z
    [etc. for total_vertices vertices]

Format:
    pp %d
    [ %g %g %g %g %g %g ] <-- for total_vertices vertices

--------

Comment.  Description:
    "#" [ string ]

Format:
    # [ string ]

    As soon as a "#" character is detected, the rest of the line is considered
    a comment.
@EOF

chmod 644 NFF

echo x - makefile
cat >makefile <<'@EOF'
# makefile for standard procedural databases
# Author:  Eric Haines, 3D/Eye, Inc.
#   Works on HP 300 and 800 machines, Apollo DN10000 - might even work for you

CC=cc -O
INC=def.h lib.h
BASELIB=-lm

all:		balls gears mount rings teapot tetra tree

lib.o:		$(INC) lib.c
		$(CC) -c lib.c

balls:		lib.o balls.o
		$(CC) -o balls balls.o lib.o $(BASELIB)

balls.o:	$(INC) balls.c
		$(CC) -c balls.c

gears:		lib.o gears.o
		$(CC) -o gears gears.o lib.o $(BASELIB)

gears.o:	$(INC) gears.c
		$(CC) -c gears.c

mount:		lib.o mount.o
		$(CC) -o mount mount.o lib.o $(BASELIB)

mount.o:	$(INC) mount.c
		$(CC) -c mount.c

rings:		lib.o rings.o
		$(CC) -o rings rings.o lib.o $(BASELIB)

rings.o:	$(INC) rings.c
		$(CC) -c rings.c

teapot:		lib.o teapot.o
		$(CC) -o teapot teapot.o lib.o $(BASELIB)

teapot.o:	$(INC) teapot.c
		$(CC) -c teapot.c

tetra:		lib.o tetra.o
		$(CC) -o tetra tetra.o lib.o $(BASELIB)

tetra.o:	$(INC) tetra.c
		$(CC) -c tetra.c

tree:		lib.o tree.o
		$(CC) -o tree tree.o lib.o $(BASELIB)

tree.o:		$(INC) tree.c
		$(CC) -c tree.c

clean:
		rm -f balls gears mount rings teapot tetra tree
		rm -f balls.o gears.o mount.o rings.o teapot.o tetra.o tree.o
		rm -f lib.o
@EOF

chmod 644 makefile

echo x - patchlevel.h
cat >patchlevel.h <<'@EOF'
#define	PATCHLEVEL	0
@EOF

chmod 644 patchlevel.h

echo x - def.h
cat >def.h <<'@EOF'
/*
 * def.h contains some useful definitions for "C" programs.
 *
 * Author:  Eric Haines, 3D/Eye, Inc.
 */

/* exit codes - define as you wish */
#define	EXIT_SUCCESS	0
#define	EXIT_FAIL	1

#define EPSILON		5.0e-6

#ifndef FALSE
#define	FALSE	0
#endif

#ifndef NULL
#define	NULL	0
#endif

#ifndef TRUE
#define	TRUE	1
#endif

#ifndef PI
#define	PI	3.141592653589793
#endif

typedef	double	COORD3[3] ;
typedef	double	COORD4[4] ;

#define	X	0
#define	Y	1
#define	Z	2
#define	W	3

typedef	double		MATRIX[4][4] ;	/* row major form */

#define ABSOLUTE(A)		( (A) < 0 ? -(A) : (A) )
#define	FRACTION(A)		( (A) - (int)(A) )
#define	IS_VAL_ALMOST_ZERO(A,E)	( ABSOLUTE(A) <= (E) )
#define	MAX(A,B)		( (A) > (B) ? (A) : (B) )
#define	MIN(A,B)		( (A) < (B) ? (A) : (B) )
#define SQR(A)			( (A) * (A) )

#define ADD2_COORD3(r,a)	{ (r)[X] += (a)[X]; (r)[Y] += (a)[Y];\
				  (r)[Z] += (a)[Z]; }
#define ADD3_COORD3(r,a,b)	{ (r)[X] = (a)[X] + (b)[X];\
				  (r)[Y] = (a)[Y] + (b)[Y];\
				  (r)[Z] = (a)[Z] + (b)[Z]; }
#define COPY_COORD3(r,a)	{ (r)[X] = (a)[X];\
				  (r)[Y] = (a)[Y];\
				  (r)[Z] = (a)[Z];}
#define COPY_COORD4(r,a)	{ (r)[X] = (a)[X];\
				  (r)[Y] = (a)[Y];\
				  (r)[Z] = (a)[Z];\
				  (r)[W] = (a)[W]; }
#define CROSS(r,a,b)		{ (r)[X] = (a)[Y] * (b)[Z] - (a)[Z] * (b)[Y];\
				  (r)[Y] = (a)[Z] * (b)[X] - (a)[X] * (b)[Z];\
				  (r)[Z] = (a)[X] * (b)[Y] - (a)[Y] * (b)[X]; }
#define DOT_PRODUCT(a,b)	( (a)[X] * (b)[X] +\
				  (a)[Y] * (b)[Y] +\
				  (a)[Z] * (b)[Z] )
#define DOT4(a,b)		( (a)[X] * (b)[X] +\
				  (a)[Y] * (b)[Y] +\
				  (a)[Z] * (b)[Z] +\
				  (a)[W] * (b)[W] )
#define IS_COORD3_ALMOST_ZERO(a,E)	(\
				     IS_VAL_ALMOST_ZERO( (a)[X], (E) )\
				  && IS_VAL_ALMOST_ZERO( (a)[Y], (E) )\
				  && IS_VAL_ALMOST_ZERO( (a)[Z], (E) ) )
#define MAX_COORD3(a)		( MAX( MAX( (a)[X], (a)[Y] ), (a)[Z] ) )
#define SET_COORD3(r,A,B,C)	{ (r)[X] = (A); (r)[Y] = (B); (r)[Z] = (C); }
#define SET_COORD4(r,A,B,C,D)	{ (r)[X] = (A); (r)[Y] = (B); (r)[Z] = (C);\
				  (r)[W] = (D); }
#define SUB2_COORD3(r,a)	{ (r)[X] -= (a)[X]; (r)[Y] -= (a)[Y];\
				  (r)[Z] -= (a)[Z]; }
#define SUB3_COORD3(r,a,b)	{ (r)[X] = (a)[X] - (b)[X];\
				  (r)[Y] = (a)[Y] - (b)[Y];\
				  (r)[Z] = (a)[Z] - (b)[Z]; }
@EOF

chmod 644 def.h

echo x - lib.h
cat >lib.h <<'@EOF'
/*
 * lib.h - vector library definitions
 *
 * Author:  Eric Haines, 3D/Eye, Inc.
 */

#define	X_AXIS	0
#define	Y_AXIS	1
#define	Z_AXIS	2

/* Output library definitions */
#define OUTPUT_CURVES		0	/* true curve output */
#define OUTPUT_PATCHES		1	/* polygonal patches output */

#define OUTPUT_RESOLUTION	4	/* amount of polygonalization */

int	lib_get_size() ;
double	lib_normalize_vector() ;
double	lib_gauss_rand() ;
@EOF

chmod 644 lib.h

echo x - lib.c
cat >lib.c <<'@EOF'
/*
 * lib.c - a library of common operations, and object output routines.
 *
 * Author:  Eric Haines, 3D/Eye, Inc.
 */

#include <stdio.h>
#include <math.h>
#include "def.h"
#include "lib.h"

/*
 * Take the first argument on the command line as a size factor.
 * No change if no/illegal argument found.
 */
int	lib_get_size( argc, argv, p_size )
int	argc ;
char	*argv[] ;
int	*p_size ;
{
int	num ;

    if ( argc >= 2 ) {
	if ( argc > 2 ) {
	    /* too many arguments */
	    return( FALSE ) ;
	}
	if ( sscanf( argv[1], "%d", &num ) == 1 ) {
	    if ( num > 0 ) {
		*p_size = num ;
	    } else {
		/* size <= 0 */
		return( FALSE ) ;
	    }
	} else {
	    /* unknown argument */
	    return( FALSE ) ;
	}
    }
    return( TRUE ) ;
}

/*
 * Normalize the vector (X,Y,Z) so that X*X + Y*Y + Z*Z = 1.
 *
 * The normalization divisor is returned.  If the divisor is zero, no
 * normalization occurs.
 *
 */
double	lib_normalize_vector( cvec )
COORD3	cvec ;
{
double	divisor ;

    divisor = sqrt( (double)DOT_PRODUCT( cvec, cvec ) ) ;

    if ( divisor != 0.0 ) {
	cvec[X] /= divisor ;
	cvec[Y] /= divisor ;
	cvec[Z] /= divisor ;
    }

    return( divisor ) ;
}

/*
 * Set all matrix elements to zero.
 */
lib_zero_matrix( mx )
MATRIX	mx ;
{
int	i, j ;

    for ( i = 0 ; i < 4 ; ++i ) {
	for ( j = 0 ; j < 4 ; ++j ) {
	    mx[i][j] = 0.0 ;
	}
    }
}

/*
 * Create identity matrix.
 */
lib_create_identity_matrix( mx )
MATRIX	mx ;
{
int	i ;

    lib_zero_matrix( mx ) ;
    for ( i = 0 ; i < 4 ; ++i ) {
	mx[i][i] = 1.0 ;
    }
}

/*
 * Create a rotation matrix along the given axis by the given angle in radians.
 */
lib_create_rotate_matrix( mx, axis, angle )
MATRIX	mx ;
int	axis ;
double	angle ;
{
double	cosine, sine ;

    lib_zero_matrix( mx ) ;

    cosine = cos( (double)angle ) ;
    sine = sin( (double)angle ) ;

    switch ( axis ) {
	case X_AXIS:
	    mx[0][0] = 1.0 ;
	    mx[1][1] = mx[2][2] = cosine ;
	    mx[1][2] = sine ;
	    mx[2][1] = -sine ;
	    break ;
	case Y_AXIS:
	    mx[1][1] = 1.0 ;
	    mx[0][0] = mx[2][2] = cosine ;
	    mx[2][0] = sine ;
	    mx[0][2] = -sine ;
	    break ;
	case Z_AXIS:
	    mx[2][2] = 1.0 ;
	    mx[0][0] = mx[1][1] = cosine ;
	    mx[0][1] = sine ;
	    mx[1][0] = -sine ;
	    break ;
    }
    mx[3][3] = 1.0 ;
}

/*
 * Create a rotation matrix along the given axis by the given angle in radians.
 * The axis is a set of direction cosines.
 */
lib_create_axis_rotate_matrix( mx, axis, angle )
MATRIX	mx ;
COORD3	axis ;
double	angle ;
{
double	cosine, one_minus_cosine, sine ;

    lib_zero_matrix( mx ) ;

    cosine = cos( (double)angle ) ;
    sine = sin( (double)angle ) ;
    one_minus_cosine = 1.0 - cosine ;

    mx[0][0] = SQR(axis[X]) + (1.0 - SQR(axis[X])) * cosine ;
    mx[0][1] = axis[X] * axis[Y] * one_minus_cosine + axis[Z] * sine ;
    mx[0][2] = axis[X] * axis[Z] * one_minus_cosine - axis[Y] * sine ;

    mx[1][0] = axis[X] * axis[Y] * one_minus_cosine - axis[Z] * sine ;
    mx[1][1] = SQR(axis[Y]) + (1.0 - SQR(axis[Y])) * cosine ;
    mx[1][2] = axis[Y] * axis[Z] * one_minus_cosine + axis[X] * sine ;

    mx[2][0] = axis[X] * axis[Z] * one_minus_cosine + axis[Y] * sine ;
    mx[2][1] = axis[Y] * axis[Z] * one_minus_cosine - axis[X] * sine ;
    mx[2][2] = SQR(axis[Z]) + (1.0 - SQR(axis[Z])) * cosine ;

    mx[3][3] = 1.0 ;
}

/*
 * Multiply a vector by a matrix.
 */
lib_transform_vector( vres, vec, mx )
COORD3	vres ;
COORD3	vec ;
MATRIX	mx ;
{
    vres[X] =
	vec[X]*mx[0][0] + vec[Y]*mx[1][0] + vec[Z]*mx[2][0] ;
    vres[Y] =
	vec[X]*mx[0][1] + vec[Y]*mx[1][1] + vec[Z]*mx[2][1] ;
    vres[Z] =
	vec[X]*mx[0][2] + vec[Y]*mx[1][2] + vec[Z]*mx[2][2] ;
}

/*
 * Multiply a point by a matrix.
 */
lib_transform_point( vres, vec, mx )
COORD3	vres ;
COORD3	vec ;
MATRIX	mx ;
{
    vres[X] =
	vec[X]*mx[0][0] + vec[Y]*mx[1][0] + vec[Z]*mx[2][0] + mx[3][0] ;
    vres[Y] =
	vec[X]*mx[0][1] + vec[Y]*mx[1][1] + vec[Z]*mx[2][1] + mx[3][1] ;
    vres[Z] =
	vec[X]*mx[0][2] + vec[Y]*mx[1][2] + vec[Z]*mx[2][2] + mx[3][2] ;
}

/*
 * Multiply a 4 element vector by a matrix.
 */
lib_transform_coord( vres, vec, mx )
COORD4	vres ;
COORD4	vec ;
MATRIX	mx ;
{
    vres[X] =
	vec[X]*mx[0][0] + vec[Y]*mx[1][0] + vec[Z]*mx[2][0] + vec[W]*mx[3][0] ;
    vres[Y] =
	vec[X]*mx[0][1] + vec[Y]*mx[1][1] + vec[Z]*mx[2][1] + vec[W]*mx[3][1] ;
    vres[Z] =
	vec[X]*mx[0][2] + vec[Y]*mx[1][2] + vec[Z]*mx[2][2] + vec[W]*mx[3][2] ;
    vres[W] =
	vec[X]*mx[0][3] + vec[Y]*mx[1][3] + vec[Z]*mx[2][3] + vec[W]*mx[3][3] ;
}

/*
 * Compute transpose of matrix.
 */
lib_transpose_matrix( mxres, mx )
MATRIX	mxres ;
MATRIX	mx ;
{
int	i, j ;

    for ( i = 0 ; i < 4 ; ++i ) {
	for ( j = 0 ; j < 4 ; ++j ) {
	    mxres[j][i] = mx[i][j] ;
	}
    }
}

/*
 * Multiply two 4x4 matrices.
 */
lib_matrix_multiply( mxres, mx1, mx2 )
MATRIX	mxres ;
MATRIX	mx1 ;
MATRIX	mx2 ;
{
int	i, j ;

    for ( i = 0 ; i < 4 ; i++ ) {
	for ( j = 0 ; j < 4 ; j++ ) {
	    mxres[i][j] = mx1[i][0]*mx2[0][j] +
			  mx1[i][1]*mx2[1][j] +
			  mx1[i][2]*mx2[2][j] +
			  mx1[i][3]*mx2[3][j] ;
	}
    }
}

/*
 * Rotate a vector pointing towards the major-axis faces (i.e. the major-axis
 * component of the vector is defined as the largest value) 90 degrees to
 * another cube face.  Mod_face is a face number.
 *
 * If the routine is called six times, with mod_face=0..5, the vector will be
 * rotated to each face of a cube.  Rotations are:
 *	mod_face = 0 mod 3, +Z axis rotate
 *	mod_face = 1 mod 3, +X axis rotate
 *	mod_face = 2 mod 3, -Y axis rotate
 */
lib_rotate_cube_face( vec, major_axis, mod_face )
COORD3	vec ;
int	major_axis ;
int	mod_face ;
{
double	swap ;

    mod_face = (mod_face+major_axis) % 3 ;
    if ( mod_face == 0 ) {
	swap   = vec[X] ;
	vec[X] = -vec[Y] ;
	vec[Y] = swap ;
    }
    else if ( mod_face == 1 ) {
	swap   = vec[Y] ;
	vec[Y] = -vec[Z] ;
	vec[Z] = swap ;
    }
    else {
	swap   = vec[X] ;
	vec[X] = -vec[Z] ;
	vec[Z] = swap ;
    }
}

/*
 * Portable gaussian random number generator (from "Numerical Recipes", GASDEV)
 * Returns a uniform random deviate between 0.0 and 1.0.  'iseed' must be
 * less than M1 to avoid repetition, and less than (2**31-C1)/A1 [= 300718]
 * to avoid overflow.
 */
#define	M1	134456
#define	IA1	8121
#define	IC1	28411
#define	RM1	1.0/M1

double	lib_gauss_rand(iseed)
int	iseed ;
{
int	ix1, ix2 ;
double	fac, r, v1, v2 ;

    ix2 = iseed ;

    do {
	ix1 = (IC1+ix2*IA1) % M1 ;
	ix2 = (IC1+ix1*IA1) % M1 ;
	v1 = ix1 * 2.0 * RM1 - 1.0 ;
	v2 = ix2 * 2.0 * RM1 - 1.0 ;
	r = v1*v1 + v2*v2 ;
    } while ( r >= 1.0 ) ;

    fac = sqrt( (double)( -2.0 * log( (double)r ) / r ) ) ;
    return( v1 * fac ) ;
}


/* OUTPUT ROUTINES
 *
 * Files are output as lines of text.  For each entity, the first line
 * defines its type.  The rest of the first line and possibly other lines
 * contain further information about the entity.  Entities include:
 *
 * "v"  - viewing vectors and angles
 * "l"  - positional light location
 * "b"  - background color
 * "f"  - object material properties
 * "c"  - cone or cylinder primitive
 * "s"  - sphere primitive
 * "p"  - polygon primitive
 * "pp" - polygonal patch primitive
 */

/*
 * Output viewpoint location.  The parameters are:
 *   From:  the eye location.
 *   At:  a position to be at the center of the image.  A.k.a. "lookat"
 *   Up:  a vector defining which direction is up.
 *
 * Note that no assumptions are made about normalizing the data (e.g. the
 * from-at distance does not have to be 1).  Also, vectors are not
 * required to be perpendicular to each other.
 *
 * For all databases some viewing parameters are always the same:
 *
 * Viewing angle is defined as from the center of top pixel row to bottom
 *	pixel row and left column to right column.
 * Yon is "at infinity."
 * Resolution is always 512 x 512.
 */
lib_output_viewpoint( from, at, up, angle, hither, resx, resy )
COORD3	from ;
COORD3	at ;
COORD3	up ;
double	angle ;
double	hither ;
int	resx ;
int	resy ;
{
    printf( "v\n" ) ;
    printf( "from %g %g %g\n", from[X], from[Y], from[Z] ) ;
    printf( "at %g %g %g\n", at[X], at[Y], at[Z] ) ;
    printf( "up %g %g %g\n", up[X], up[Y], up[Z] ) ;
    printf( "angle %g\n", angle ) ;
    printf( "hither %g\n", hither ) ;
    printf( "resolution %d %d\n", resx, resy ) ;
}

/*
 * Output light.  A light is defined by position.  All lights have the same
 * intensity.
 *
 */
lib_output_light( center_pt )
COORD3 	center_pt ;
{
    printf( "l %g %g %g\n", center_pt[X], center_pt[Y], center_pt[Z] ) ;
}

/*
 * Output background color.  A color is simply RGB (monitor dependent, but
 * that's life).  The format is:
 *	"b" red green blue
 */
lib_output_background_color( color )
COORD3 	color ;
{
    printf( "b %g %g %g\n", color[X], color[Y], color[Z] ) ;
}

/*
 * Output a color and shading parameters for the object in the format:
 *	"f" red green blue Kd Ks Shine T index_of_refraction
 *
 * Kd is the diffuse component, Ks the specular, Shine is the Phong cosine
 * power for highlights, T is transmittance (fraction of light passed per
 * unit).  0 <= Kd <= 1 and 0 <= Ks <= 1, though it is not required that
 * Kd + Ks == 1.
 *
 * The fill color is used to color the objects following it until a new color
 * is assigned or the file ends.
 */
lib_output_color( color, kd, ks, shine, t, i_of_r )
COORD3 	color ;
double	kd ;
double	ks ;
double	shine ;
double	t ;
double	i_of_r ;
{
    printf( "f %g %g %g %g %g %g %g %g\n", color[X], color[Y], color[Z],
						kd, ks, shine, t, i_of_r ) ;
}

/*
 * Output cylinder or cone.  A cylinder is defined as having a radius and an
 * axis defined by two points, which also define the top and bottom edge of the
 * cylinder.  A cone is defined similarly, the difference being that the apex
 * and base radii are different.  The apex radius is defined as being smaller
 * than the base radius.  Note that the surface exists without endcaps.
 *
 * If format=OUTPUT_CURVES, output the cylinder/cone in format:
 *	"c"
 *	base[X] base[Y] base[Z] base_radius
 *	apex[X] apex[Y] apex[Z] apex_radius
 *
 * If the format=OUTPUT_PATCHES, the surface is polygonalized and output.
 * (4*OUTPUT_RESOLUTION) polygons are output as rectangles by
 * lib_output_polypatch.
 */
lib_output_cylcone( base_pt, apex_pt, format )
COORD4	base_pt ;
COORD4	apex_pt ;
int	format ;
{
int	num_pol ;
double	angle, divisor ;
COORD3	axis, dir ;
COORD3	norm_axis, start_norm, start_radius[4] ;
COORD3	lip_norm[4], lip_pt[4] ;
MATRIX	mx ;

    if ( format == OUTPUT_CURVES ) {
	printf( "c\n" ) ;
	printf( "%g %g %g %g\n",
		base_pt[X], base_pt[Y], base_pt[Z], base_pt[W] ) ;
	printf( "%g %g %g %g\n",
		apex_pt[X], apex_pt[Y], apex_pt[Z], apex_pt[W] ) ;
    }
    else {
	SUB3_COORD3( axis, apex_pt, base_pt ) ;
	COPY_COORD3( norm_axis, axis ) ;
	(void)lib_normalize_vector( norm_axis ) ;

	SET_COORD3( dir, 0.0, 0.0, 1.0 ) ;
	CROSS( start_norm, axis, dir ) ;

	divisor = lib_normalize_vector( start_norm ) ;
	if ( ABSOLUTE( divisor ) < EPSILON ) {
	    SET_COORD3( dir, 1.0, 0.0, 0.0 ) ;
	    CROSS( start_norm, axis, dir ) ;
	    (void)lib_normalize_vector( start_norm ) ;
	}

	start_radius[0][X] = start_norm[X] * base_pt[W] ;
	start_radius[0][Y] = start_norm[Y] * base_pt[W] ;
	start_radius[0][Z] = start_norm[Z] * base_pt[W] ;
	ADD3_COORD3( lip_pt[0], base_pt, start_radius[0] ) ;

	start_radius[1][X] = start_norm[X] * apex_pt[W] ;
	start_radius[1][Y] = start_norm[Y] * apex_pt[W] ;
	start_radius[1][Z] = start_norm[Z] * apex_pt[W] ;
	ADD3_COORD3( lip_pt[1], apex_pt, start_radius[1] ) ;

	COPY_COORD3( lip_norm[0], start_norm ) ;
	COPY_COORD3( lip_norm[1], start_norm ) ;

	for ( num_pol = 0 ; num_pol < 4*OUTPUT_RESOLUTION ; ++num_pol ) {
	    COPY_COORD3( lip_pt[3], lip_pt[0] ) ;
	    COPY_COORD3( lip_pt[2], lip_pt[1] ) ;
	    COPY_COORD3( lip_norm[3], lip_norm[0] ) ;
	    COPY_COORD3( lip_norm[2], lip_norm[1] ) ;

	    angle = 2.0 * PI *
		    (double)( num_pol+1 ) / (double)( 4*OUTPUT_RESOLUTION ) ;
	    lib_create_axis_rotate_matrix( mx, norm_axis, angle ) ;

	    lib_transform_vector( lip_pt[0], start_radius[0], mx ) ;
	    ADD2_COORD3( lip_pt[0], base_pt ) ;
	    lib_transform_vector( lip_pt[1], start_radius[1], mx ) ;
	    ADD2_COORD3( lip_pt[1], apex_pt ) ;

	    lib_transform_vector( lip_norm[0], start_norm, mx ) ;
	    COPY_COORD3( lip_norm[1], lip_norm[0] ) ;

	    lib_output_polypatch( 4, lip_pt, lip_norm ) ;
	}
    }
}

/*
 * Output sphere.  A sphere is defined by a radius and center position.
 *
 * If format=OUTPUT_CURVES, output the sphere in format:
 *	"s" center[X] center[Y] center[Z] radius
 *
 * If the format=OUTPUT_PATCHES, the sphere is polygonalized and output.
 * The sphere is polygonalized by splitting it into 6 faces (of a cube
 * projected onto the sphere) and dividing these faces by equally spaced
 * great circles.  OUTPUT_RESOLUTION affects the number of great circles.
 * (6*2*OUTPUT_RESOLUTION*OUTPUT_RESOLUTION) polygons are output as triangles
 * using lib_output_polypatch.
 */
lib_output_sphere( center_pt, format )
COORD4	center_pt ;
int	format ;
{
int	num_face, num_edge, num_tri, num_vert, u_pol, v_pol ;
double	angle ;
COORD3	edge_norm[3], edge_pt[3] ;
COORD3	x_axis[OUTPUT_RESOLUTION+1], y_axis[OUTPUT_RESOLUTION+1] ;
COORD3	pt[OUTPUT_RESOLUTION+1][OUTPUT_RESOLUTION+1] ;
COORD3	mid_axis ;
MATRIX	rot_mx ;

    if ( format == OUTPUT_CURVES ) {
	printf( "s %g %g %g %g\n",
		center_pt[X], center_pt[Y], center_pt[Z], center_pt[W] ) ;
    }
    else {
	/* calculate axes used to find grid points */
	for ( num_edge = 0 ; num_edge <= OUTPUT_RESOLUTION ; ++num_edge ) {
	    angle = (PI/4.0) * (2.0*(double)num_edge/OUTPUT_RESOLUTION - 1.0) ;

	    mid_axis[X] = 1.0 ; mid_axis[Y] = 0.0 ; mid_axis[Z] = 0.0 ;
	    lib_create_rotate_matrix( rot_mx, Y_AXIS, angle ) ;
	    lib_transform_vector( x_axis[num_edge], mid_axis, rot_mx ) ;

	    mid_axis[X] = 0.0 ; mid_axis[Y] = 1.0 ; mid_axis[Z] = 0.0 ;
	    lib_create_rotate_matrix( rot_mx, X_AXIS, angle ) ;
	    lib_transform_vector( y_axis[num_edge], mid_axis, rot_mx ) ;
	}

	/* set up grid of points on +Z sphere surface */
	for ( u_pol = 0 ; u_pol <= OUTPUT_RESOLUTION ; ++u_pol ) {
	    for ( v_pol = 0 ; v_pol <= OUTPUT_RESOLUTION ; ++v_pol ) {
		CROSS( pt[u_pol][v_pol], x_axis[u_pol], y_axis[v_pol] ) ;
		(void)lib_normalize_vector( pt[u_pol][v_pol] ) ;
	    }
	}
	for ( num_face = 0 ; num_face < 6 ; ++num_face ) {
	    /* transform points to cube face */
	    for ( u_pol = 0 ; u_pol <= OUTPUT_RESOLUTION ; ++u_pol ) {
		for ( v_pol = 0 ; v_pol <= OUTPUT_RESOLUTION ; ++v_pol ) {
		    lib_rotate_cube_face( pt[u_pol][v_pol]
					, Z_AXIS
					, num_face
					) ;
		}
	    }
	    /* output grid */
	    for ( u_pol = 0 ; u_pol < OUTPUT_RESOLUTION ; ++u_pol ) {
		for ( v_pol = 0 ; v_pol < OUTPUT_RESOLUTION ; ++v_pol ) {
		    for ( num_tri = 0 ; num_tri < 2 ; ++num_tri ) {
			for ( num_edge = 0 ; num_edge < 3 ; ++num_edge ) {
			    num_vert = (num_tri*2 + num_edge) % 4 ;
			    if ( num_vert == 0 ) {
				COPY_COORD3( edge_pt[num_edge]
					   , pt[u_pol][v_pol] ) ;
			    }
			    else if ( num_vert == 1 ) {
				COPY_COORD3( edge_pt[num_edge]
					   , pt[u_pol][v_pol+1] ) ;
			    }
			    else if ( num_vert == 2 ) {
				COPY_COORD3( edge_pt[num_edge]
					   , pt[u_pol+1][v_pol+1] ) ;
			    }
			    else {
				COPY_COORD3( edge_pt[num_edge]
					   , pt[u_pol+1][v_pol] ) ;
			    }
			    COPY_COORD3( edge_norm[num_edge]
				       , edge_pt[num_edge] ) ;
			    edge_pt[num_edge][X] =
				    edge_pt[num_edge][X] * center_pt[W] +
							   center_pt[X] ;
			    edge_pt[num_edge][Y] =
				    edge_pt[num_edge][Y] * center_pt[W] +
							   center_pt[Y] ;
			    edge_pt[num_edge][Z] =
				    edge_pt[num_edge][Z] * center_pt[W] +
							   center_pt[Z] ;

			}
			lib_output_polypatch( 3, edge_pt, edge_norm ) ;
		    }
		}
	    }
	}
    }
}

/*
 * Output polygon.  A polygon is defined by a set of vertices.  With these
 * databases, a polygon is defined to have all points coplanar.  A polygon has
 * only one side, with the order of the vertices being counterclockwise as you
 * face the polygon (right-handed coordinate system).
 *
 * The output format is always:
 *	"p" total_vertices
 *	vert1[X] vert1[Y] vert1[Z]
 *	[etc. for total_vertices polygons]
 *
 */
lib_output_polygon( tot_vert, vert )
int	tot_vert ;
COORD3 	vert[] ;
{
int	num_vert ;

    printf( "p %d\n", tot_vert ) ;

    for ( num_vert = 0 ; num_vert < tot_vert ; ++num_vert ) {
	printf( "%g %g %g\n", vert[num_vert][X]
			    , vert[num_vert][Y]
			    , vert[num_vert][Z] ) ;
    }
}

/*
 * Output polygonal patch.  A patch is defined by a set of vertices and their
 * normals.  With these databases, a patch is defined to have all points
 * coplanar.  A patch has only one side, with the order of the vertices being
 * counterclockwise as you face the patch (right-handed coordinate system).
 *
 * The output format is always:
 *	"pp" total_vertices
 *	vert1[X] vert1[Y] vert1[Z] norm1[X] norm1[Y] norm1[Z]
 *	[etc. for total_vertices polygonal patches]
 *
 */
lib_output_polypatch( tot_vert, vert, norm )
int	tot_vert ;
COORD3	vert[] ;
COORD3	norm[] ;
{
int	num_vert ;

    printf( "pp %d\n", tot_vert ) ;

    for ( num_vert = 0 ; num_vert < tot_vert ; ++num_vert ) {
	printf( "%g %g %g %g %g %g\n", vert[num_vert][X]
				     , vert[num_vert][Y]
				     , vert[num_vert][Z]
				     , norm[num_vert][X]
				     , norm[num_vert][Y]
				     , norm[num_vert][Z] ) ;
    }
}
@EOF

chmod 644 lib.c

echo x - balls.c
cat >balls.c <<'@EOF'
/*
 * balls.c - Also known as "sphereflake".  Create a set of shiny spheres, with
 *	each sphere blooming sets of 9 more spheres with 1/3rd radius.
 *	A square floor polygon is added.  Three light sources.
 *
 * Author:  Eric Haines, 3D/Eye, Inc.
 *
 * SizeFactor determines the number of objects output.
 *	Total spheres = sum of n=0,SF of (9**SF).
 *
 *	SizeFactor	# spheres	# squares
 *	     1		    10		     1
 *	     2		    91		     1
 *	     3		   820		     1
 *
 *	     4		  7381		     1
 */

#include <stdio.h>
#include <math.h>
#include "def.h"
#include "lib.h"

static	int	SizeFactor = 4 ;

#define	OUTPUT_FORMAT		OUTPUT_CURVES

/* radius of sphere enclosing all balls (a scale factor) */
#define	BV_RADIUS		1.0

main(argc,argv)
int	argc ;
char	*argv[] ;
{
COORD3	back_color, obj_color ;
COORD3	from, at, up ;
COORD3	ground[4], bvec, light, objset[9] ;
COORD4	center_pt, direction ;

    if ( !lib_get_size( argc, argv, &SizeFactor ) ) {
	fprintf( stderr, "usage: %s [size]\n", *argv ) ;
	exit(EXIT_FAIL) ;
    }

    /* output viewpoint */
    SET_COORD3( from, 2.1, 1.3, 1.7 ) ;
    SET_COORD3( at, 0.0, 0.0, 0.0 ) ;
    SET_COORD3( up, 0.0, 0.0, 1.0 ) ;
    lib_output_viewpoint( from, at, up, 45.0, 1.0, 512, 512 ) ;

    /* output background color - UNC sky blue */
    SET_COORD3( back_color, 0.078, 0.361, 0.753 ) ;
    lib_output_background_color( back_color ) ;

    /* output light sources */
    SET_COORD3( light, 4.0, 3.0, 2.0 ) ;
    lib_output_light( light ) ;
    SET_COORD3( light, 1.0, -4.0, 4.0 ) ;
    lib_output_light( light ) ;
    SET_COORD3( light, -3.0, 1.0, 5.0 ) ;
    lib_output_light( light ) ;

    /* output floor polygon - beige */
    SET_COORD3( back_color, 1.0, 0.75, 0.33 ) ;
    lib_output_color( back_color, 1.0, 0.0, 0.0, 0.0, 0.0 ) ;
    SET_COORD3( bvec, 12.0 * BV_RADIUS, 12.0 * BV_RADIUS, -BV_RADIUS / 2.0 ) ;
    SET_COORD3( ground[0],  bvec[X],  bvec[Y], bvec[Z] ) ;
    SET_COORD3( ground[1], -bvec[X],  bvec[Y], bvec[Z] ) ;
    SET_COORD3( ground[2], -bvec[X], -bvec[Y], bvec[Z] ) ;
    SET_COORD3( ground[3],  bvec[X], -bvec[Y], bvec[Z] ) ;
    lib_output_polygon( 4, ground ) ;

    /* set up object color - off white */
    SET_COORD3( obj_color, 1.0, 0.9, 0.7 ) ;
    lib_output_color( obj_color, 0.5, 0.5, 3.0, 0.0, 0.0 ) ;

    /* create set of spawned points */
    create_objset( objset ) ;

    /* compute and output object */
    SET_COORD4( center_pt, 0.0, 0.0, 0.0, BV_RADIUS / 2.0 ) ;
    SET_COORD4( direction, 0.0, 0.0, 1.0, 1.0/3.0 ) ;
    output_object( SizeFactor, center_pt, direction, objset ) ;

    exit(EXIT_SUCCESS) ;
}

/* Create the set of 9 vectors needed to generate the sphere set. */
create_objset( objset )
COORD3	objset[9] ;
{
int	num_set, num_vert ;
double	dist ;
COORD3	axis, temp_pt, trio_dir[3] ;
MATRIX	mx ;

    dist = 1.0 / sqrt( (double)2.0 ) ;

    SET_COORD3( trio_dir[0], dist, dist,   0.0 ) ;
    SET_COORD3( trio_dir[1], dist,  0.0, -dist ) ;
    SET_COORD3( trio_dir[2],  0.0, dist, -dist ) ;

    SET_COORD3( axis, 1.0, -1.0, 0.0 ) ;
    (void)lib_normalize_vector( axis ) ;
    lib_create_axis_rotate_matrix(
	    mx
	  , axis
	  , asin( (double) ( 2.0 / sqrt( (double)6.0 ) ) ) ) ;

    for ( num_vert = 0 ; num_vert < 3 ; ++num_vert ) {
	lib_transform_vector( temp_pt, trio_dir[num_vert], mx ) ;
	COPY_COORD3( trio_dir[num_vert], temp_pt ) ;
    }

    for ( num_set = 0 ; num_set < 3 ; ++num_set ) {
	lib_create_rotate_matrix( mx, Z_AXIS, num_set*2.0*PI/3.0 ) ;
	for ( num_vert = 0 ; num_vert < 3 ; ++num_vert ) {
	    lib_transform_vector( objset[num_set*3+num_vert],
		    trio_dir[num_vert], mx ) ;
	}
    }
}


/*
 * Output the parent sphere, then output the children of the sphere.
 */
output_object( depth, center, direction, objset )
int	depth ;
COORD4	center ;
COORD4	direction ;
COORD3	objset[9] ;
{
int	num_vert ;
double	angle, scale ;
COORD3	axis, z_axis ;
COORD4	child_pt, child_dir ;
MATRIX	mx ;

    /* output sphere at location & radius defined by center */
    lib_output_sphere( center, OUTPUT_FORMAT ) ;

    /* check if children should be generated */
    if ( depth > 0 ) {
	--depth ;

	/* rotation matrix to new axis from +Z axis */
	if ( direction[Z] >= 1.0 ) {
	    /* identity matrix */
	    lib_create_identity_matrix( mx ) ;
	}
	else if ( direction[Z] <= -1.0 ) {
	    lib_create_rotate_matrix( mx, Y_AXIS, PI ) ;
	}
	else {
	    SET_COORD3( z_axis, 0.0, 0.0, 1.0 ) ;
	    CROSS( axis, z_axis, direction ) ;
	    (void)lib_normalize_vector( axis ) ;
	    angle = acos( (double)DOT_PRODUCT( z_axis, direction ) ) ;
	    lib_create_axis_rotate_matrix( mx, axis, angle ) ;
	}

	/* scale down location of new spheres */
	scale = center[W] * (1.0 + direction[W] ) ;

	for ( num_vert = 0 ; num_vert < 9 ; ++num_vert ) {
	    lib_transform_vector( child_pt, objset[num_vert], mx ) ;
	    child_pt[X] = child_pt[X] * scale + center[X] ;
	    child_pt[Y] = child_pt[Y] * scale + center[Y] ;
	    child_pt[Z] = child_pt[Z] * scale + center[Z] ;
	    /* scale down radius */
	    child_pt[W] = center[W] * direction[W] ;
	    SUB3_COORD3( child_dir, child_pt, center ) ;
	    child_dir[X] /= scale ;
	    child_dir[Y] /= scale ;
	    child_dir[Z] /= scale ;
	    child_dir[W] = direction[W] ;
	    output_object( depth, child_pt, child_dir, objset ) ;
	}
    }
}
@EOF

chmod 644 balls.c

echo x - gears.c
cat >gears.c <<'@EOF'
/*
 * gears.c - Create a set of gears.  Each gear face has 144 vertices, and
 *	contains concavities.  Note that the first 3 vertices of all polygons
 *	define the two edges of a convex section of the polygon.  Background
 *	square floor is reflective.  Five light sources.  
 *
 * Author:  Eric Haines, 3D/Eye, Inc.
 *
 * SizeFactor determines the number of polygons output.
 *	Total gears = SF**3:  concave polygons = 2 * SF**3
 *	rectangles = 4*TEETH * SF**3
 *
 *	SizeFactor	# gears	   # gear faces    # rectangles
 *	     1		    1		  2		144
 *	     2		    8		 16	       1152
 *	     3		   27		 54	       3888
 *
 *	     4		   64		128	       9216
 */

#include <stdio.h>
#include <math.h>
#include "def.h"
#include "lib.h"

static	int	SizeFactor = 4 ;

/* define number of teeth on a gear - must be a multiple of 4 */
#define	TEETH			36
/* define ratio of radius taken up by teeth and the gear thickness */
/* the outer radius is made slightly smaller that the full radius to create
 * a finite separation between intermeshing gears.  This gets rid of the bug
 * of having two surfaces occupy exactly the same space.  Note that if these
 * are changed, the gears may interpenetrate.
 */
#define	OUTER_EDGE_RATIO	0.995
#define	INNER_EDGE_RATIO	0.9
#define	EDGE_DIFF		( 1.0 - INNER_EDGE_RATIO )
/* ratio of width of gear to thickness */
#define	DEPTH_RATIO		0.1

main(argc,argv)
int	argc ;
char	*argv[] ;
{
int	ix, iy, iz ;
double	angle, color_scale, outer_radius, thickness ;
COORD3	back_color, gear_color ;
COORD3	from, at, up ;
COORD3	ground[4], light, offset, zero_pt, center ;

    if ( !lib_get_size( argc, argv, &SizeFactor ) ) {
	fprintf( stderr, "usage: %s [size]\n", *argv ) ;
	exit(EXIT_FAIL) ;
    }

    /* output viewpoint */
    SET_COORD3( from, -1.1, -2.1, 2.6 ) ;
    SET_COORD3( at, 0.0, 0.0, 0.0 ) ;
    SET_COORD3( up, 0.0, 0.0, 1.0 ) ;
    lib_output_viewpoint( from, at, up, 45.0, 1.0, 512, 512 ) ;

    /* output background color - UNC sky blue */
    SET_COORD3( back_color, 0.078, 0.361, 0.753 ) ;
    lib_output_background_color( back_color ) ;

    /* output light sources */
    SET_COORD3( light, 2.0, 4.0, 4.0 ) ;
    lib_output_light( light ) ;
    SET_COORD3( light, -2.0, 4.0, 3.0 ) ;
    lib_output_light( light ) ;
    SET_COORD3( light, 2.0, -2.5, 2.5 ) ;
    lib_output_light( light ) ;
    SET_COORD3( light, -1.0, -4.0, 2.0 ) ;
    lib_output_light( light ) ;
    /* just behind the eye */
    SET_COORD3( light, -1.111, -2.121, 2.626 ) ;
    lib_output_light( light ) ;

    /* output ground polygon - off white */
    SET_COORD3( back_color, 1.0, 0.85, 0.7 ) ;
    lib_output_color( back_color, 0.75, 0.25, 25.0, 0.0, 0.0 ) ;
    SET_COORD3( ground[0], 2.0, 2.0, 0.0 ) ;
    SET_COORD3( ground[1], -2.0, 2.0, 0.0 ) ;
    SET_COORD3( ground[2], -2.0, -2.0, 0.0 ) ;
    SET_COORD3( ground[3], 2.0, -2.0, 0.0 ) ;
    lib_output_polygon( 4, ground ) ;

    outer_radius = 1.0 /
	( (double)SizeFactor - (double)(SizeFactor-1) * EDGE_DIFF / 2.0 ) ;
    /* calculate first gear center */
    zero_pt[X] = zero_pt[Y] = -1.0 + outer_radius ;
    zero_pt[Z] = 1.0 ;
    /* calculate offset */
    offset[X] = offset[Y] = outer_radius * ( 2.0 - EDGE_DIFF ) ;
    offset[Z] = -1.0 / (double)SizeFactor ;

    /* create gears */
    for ( iz = 0 ; iz < SizeFactor ; ++iz ) {
	center[Z] = zero_pt[Z] + (double)iz * offset[Z] ;
	for ( iy = 0 ; iy < SizeFactor ; ++iy ) {
	    center[Y] = zero_pt[Y] + (double)iy * offset[Y] ;
	    for ( ix = 0 ; ix < SizeFactor ; ++ix ) {
		center[X] = zero_pt[X] + (double)ix * offset[X] ;

		/* output pseudo-random gear color */
		SET_COORD3( gear_color
			  , 0.01 + FRACTION( (double)(ix*3+iy*2+iz+1)*5.0/7.0 )
			  , 0.01 + FRACTION( (double)(iy*3+iz*2+ix+1)*3.0/7.0 )
			  , 0.01 + FRACTION( (double)(iz*3+ix*2+iy+1)*2.0/7.0 )
			  ) ;
		color_scale = MAX_COORD3( gear_color ) ;
		gear_color[X] /= color_scale ;
		gear_color[Y] /= color_scale ;
		gear_color[Z] /= color_scale ;
		if ( ( ix*4 + iy*2 + iz ) % 5 == 0 ) {
		    lib_output_color( gear_color,
						0.75, 0.25, 50.0, 0.95, 1.1 ) ;
		}
		else {
		    lib_output_color( gear_color,
						1.0, 0.0, 0.0, 0.0, 0.0 ) ;
		}

		/* output gear */
		angle = PI * (double)( (ix+iy+iz) % 2 ) / (double)(TEETH) ;
		thickness =
		    MIN( DEPTH_RATIO, 1.0 / ( 2.0 * (double)SizeFactor ) ) ;
		create_gear( center,
			     angle,
			     OUTER_EDGE_RATIO * outer_radius,
			     (1.0 - EDGE_DIFF) * outer_radius,
			     thickness ) ;
	    }
	}
    }
    exit(EXIT_SUCCESS) ;
}

/* Create gear */
create_gear( center, offset_angle, outer_radius, inner_radius, thickness )
COORD3	center ;
double	offset_angle ;
double	outer_radius ;
double	inner_radius ;
double	thickness ;
{
int	next_side, num_side, num_teeth ;
double	gear_angle, tooth_angle ;
COORD3	side_pts[4], gear_pts[4*TEETH], outer_pt, inner_pt ;

    SET_COORD3( outer_pt, outer_radius, 0.0, 0.0 ) ;
    SET_COORD3( inner_pt, inner_radius, 0.0, 0.0 ) ;

    tooth_angle = 2.0 * PI / (double)TEETH ;

    /* output gear top */
    for ( num_teeth = 0 ; num_teeth < TEETH ; ++num_teeth ) {
	gear_angle = offset_angle +
		2.0 * PI * (double)num_teeth / (double)TEETH ;
	create_tooth( gear_angle
		    , tooth_angle
		    , center
		    , outer_pt
		    , inner_pt
		    , (COORD3 *)gear_pts[num_teeth*4]
		    ) ;
    }
    lib_output_polygon( 4*TEETH, gear_pts ) ;

    /* output teeth */
    for ( num_side = 0 ; num_side < 4 * TEETH ; ++num_side ) {
	next_side = ( num_side + 1 ) % ( 4 * TEETH ) ;
	COPY_COORD3( side_pts[0], gear_pts[num_side] ) ;
	COPY_COORD3( side_pts[1], gear_pts[num_side] ) ;
	side_pts[1][Z] -= thickness ;
	COPY_COORD3( side_pts[2], gear_pts[next_side] ) ;
	side_pts[2][Z] -= thickness ;
	COPY_COORD3( side_pts[3], gear_pts[next_side] ) ;
	lib_output_polygon( 4, side_pts ) ;
    }

    /* output gear bottom */
    outer_pt[Z] = inner_pt[Z] = -thickness ;
    for ( num_teeth = 0 ; num_teeth < TEETH ; ++num_teeth ) {
	gear_angle = offset_angle -
		2.0 * PI * (double)num_teeth / (double)TEETH ;
	create_tooth( gear_angle
		    , -tooth_angle
		    , center
		    , outer_pt
		    , inner_pt
		    , (COORD3 *)gear_pts[num_teeth*4]
		    ) ;
    }
    lib_output_polygon( 4*TEETH, gear_pts ) ;
}

/* Create gear tooth */
create_tooth( gear_angle, tooth_angle, center, outer_pt, inner_pt, edge_pts )
double	gear_angle ;
double	tooth_angle ;
COORD3	center ;
COORD3	outer_pt ;
COORD3	inner_pt ;
COORD3	edge_pts[] ;
{
MATRIX	mx ;

    lib_create_rotate_matrix( mx
			    , Z_AXIS
			    , gear_angle - 0.19 * tooth_angle ) ;
    lib_transform_coord( edge_pts[0], outer_pt, mx ) ;
    ADD2_COORD3( edge_pts[0], center ) ;
    lib_create_rotate_matrix( mx
			    , Z_AXIS
			    , gear_angle + 0.19 * tooth_angle ) ;
    lib_transform_coord( edge_pts[1], outer_pt, mx ) ;
    ADD2_COORD3( edge_pts[1], center ) ;
    lib_create_rotate_matrix( mx
			    , Z_AXIS
			    , gear_angle + 0.3 * tooth_angle ) ;
    lib_transform_coord( edge_pts[2], inner_pt, mx ) ;
    ADD2_COORD3( edge_pts[2], center ) ;
    lib_create_rotate_matrix( mx
			    , Z_AXIS
			    , gear_angle + 0.7 * tooth_angle ) ;
    lib_transform_coord( edge_pts[3], inner_pt, mx ) ;
    ADD2_COORD3( edge_pts[3], center ) ;
}
@EOF

chmod 644 gears.c

echo x - mount.c
cat >mount.c <<'@EOF'
/*
 * mount.c - creates a fractal mountain, using Carpenter's method with a
 *	different extension to square grids.  A pyramid of 4 glass spheres
 *	is added in front of the mountain.  One light source.
 *
 *	NOTE: the hashing function used to generate the database originally is
 *	faulty.  The function causes repetition to occur within the fractal
 *	mountain (obviously not very fractal behavior!).  A new hashing
 *	function is included immediately after the old one:  merely define
 *	NEW_HASH if you want to use a better hashing function.  To perform ray
 *	tracing comparison tests you should still use the old, faulty database
 *	(it may have repetition, but it's still a good test image).
 *     
 * Author:  Eric Haines, 3D/Eye, Inc.
 *
 * SizeFactor determines the number of objects output.
 *	Total triangular polygons = 2 * (4**SizeFactor)
 *
 *	SizeFactor	# triangles	# spheres
 *	     1		     8		     4
 *	     2		    32		     4
 *	     3		   128		     4
 *
 *	     6		  8192		     4
 */

#include <stdio.h>
#include <math.h>
#include "def.h"
#include "lib.h"

static	int	SizeFactor = 6 ;

#define	OUTPUT_FORMAT		OUTPUT_CURVES

/* to use the corrected hashing function, uncomment this next line */
/* #define NEW_HASH */

/* fractal dimension - affects variance of z.  Between 2 and 3 */
#define	FRACTAL_DIMENSION	2.2
/* change MOUNTAIN_NO to get a different mountain */
#define	MOUNTAIN_NO		21

/* lower left corner and width of mountain definitions */
#define	X_CORNER	-1.0
#define	Y_CORNER	-1.0
#define	WIDTH		 2.0

#ifndef NEW_HASH

/* Hashing function to get a seed for the random number generator. */
/* This is the old, buggy hashing function - use it if you wish to
 * obtain the same image as in the November 1987 IEEE CG&A article. */
#define	hash_rand(A,B,C)	( ( (((unsigned long)(A))<<(23-(C))) +	\
				    (((unsigned long)(B))<<(15-(C)))	\
				  + (((unsigned long)(A))<<(7-(C))) ) & 0xffff)
#else

/* New, corrected hashing function.  Use for a true fractal mountain */
/* 134456 is M1 in routine lib_gauss_rand() */
#define	hash_rand(A,B,C)	( ( C <= 15 ) ?				\
				  ( ABSOLUTE(				\
				    (((unsigned long)(A))<<(31-(C))) +	\
				    (((unsigned long)(B))<<(15-(C))) )	\
				  % 134456 ) :				\
				  ( ABSOLUTE(				\
				    (((unsigned long)(A))<<(31-(C))) +	\
				    (((unsigned long)(B))>>((C)-15)) )	\
				  % 134456 ) )
				  )
#endif

static  double	Roughness ;

main(argc,argv)
int	argc ;
char	*argv[] ;
{
int	num_pts ;
double	ratio ;
COORD3	back_color, obj_color ;
COORD3	from, at, up, light ;
COORD4	center ;

    if ( !lib_get_size( argc, argv, &SizeFactor ) ) {
	fprintf( stderr, "usage: %s [size]\n", *argv ) ;
	exit(EXIT_FAIL) ;
    }

    /* output viewpoint */
    SET_COORD3( from, -1.6, 1.6, 1.7 ) ;
    SET_COORD3( at, 0.0, 0.0, 0.0 ) ;
    SET_COORD3( up, 0.0, 0.0, 1.0 ) ;
    lib_output_viewpoint( from, at, up, 45.0, 0.5, 512, 512 ) ;

    /* output background color - UNC sky blue */
    SET_COORD3( back_color, 0.078, 0.361, 0.753 ) ;
    lib_output_background_color( back_color ) ;

    /* output light sources */
    SET_COORD3( light, -100.0, -100.0, 100.0 ) ;
    lib_output_light( light ) ;

    /* set up crystal sphere color - clear white */
    SET_COORD3( obj_color, 1.0, 1.0, 1.0 ) ;
    lib_output_color( obj_color, 0.1, 0.9, 100.0, 0.9, 1.5 ) ;

    /* output crystal spheres */
    SET_COORD4( center, -0.8, 0.8, 1.00, 0.17 ) ;
    create_spheres( center ) ;

    /* set up mountain color - grey */
    SET_COORD3( obj_color, 0.5, 0.45, 0.35 ) ;
    lib_output_color( obj_color, 1.0, 0.0, 0.0, 0.0, 0.0 ) ;

    /* grow mountain */
    num_pts = 1<<SizeFactor ;
    ratio = 2.0 /
	    exp( (double)( log( (double)2.0 ) / (FRACTAL_DIMENSION-1.0) ) ) ;
    Roughness = sqrt( (double)( SQR(ratio) - 1.0 ) ) ;
    grow_mountain( (double)num_pts, num_pts, 0, 0, 0.0, 0.0, 0.0, 0.0 ) ;

    exit(EXIT_SUCCESS) ;
}

/* create a pyramid of crystal spheres */
create_spheres( center )
COORD4	center ;
{
int	i ;
double	angle ;
COORD3	axis, pt, new_pt ;
COORD4	sphere ;
MATRIX	mx ;

    SET_COORD3( axis, 1.0, 1.0, 0.0 ) ;
    (void)lib_normalize_vector( axis ) ;
    angle = acos( (double)( -1.0/3.0 ) ) ;

    /* set center of pyramid */
    SET_COORD3( pt, 0.0, 0.0, center[W] * sqrt( (double)( 3.0/2.0 ) ) ) ;

    COPY_COORD4( sphere, center ) ;
    ADD2_COORD3( sphere, pt ) ;
    lib_output_sphere( sphere, OUTPUT_FORMAT ) ;

    lib_create_axis_rotate_matrix( mx, axis, angle ) ;
    lib_transform_vector( new_pt, pt, mx ) ;

    for ( i = 0 ; i < 3 ; ++i ) {
	lib_create_rotate_matrix( mx, Z_AXIS, (double)i * 2.0 * PI / 3.0 ) ;
	lib_transform_vector( sphere, new_pt, mx ) ;
	ADD2_COORD3( sphere, center ) ;
	lib_output_sphere( sphere, OUTPUT_FORMAT ) ;
    }
}

/*
 * Build mountain section.  If at width > 1, split quadrilateral into four
 * parts.  Else if at width == 1, output quadrilateral as two triangles.
 */
grow_mountain( fnum_pts, width, ll_x, ll_y, ll_fz, lr_fz, ur_fz, ul_fz )
double	fnum_pts ;
int	width ;
int	ll_x ;
int	ll_y ;
double	ll_fz ;
double	lr_fz ;
double	ur_fz ;
double	ul_fz ;
{
int	half_width, iz, num_tri, num_tri_vert, num_vert ;
double	l_fx, r_fx, l_fy, u_fy ;
double	lower_fz, right_fz, upper_fz, left_fz, middle_fz ;
double	rise_height, hside_length ;
COORD3	tri_vert[3] ;

    if ( width == 1 ) {
	/* calculate x and y coordinates of corners */
	l_fx = X_CORNER + (double)ll_x * WIDTH / fnum_pts ;
	r_fx = X_CORNER + (double)(ll_x+1) * WIDTH / fnum_pts ;
	l_fy = Y_CORNER + (double)ll_y * WIDTH / fnum_pts ;
	u_fy = Y_CORNER + (double)(ll_y+1) * WIDTH / fnum_pts ;

	/* output two triangles for section */
	for ( num_tri = 0 ; num_tri < 2 ; ++num_tri ) {
	    for ( num_vert = 0 ; num_vert < 3 ; ++num_vert ) {
		num_tri_vert = ( num_vert + num_tri * 2 ) % 4 ;
		switch ( num_tri_vert ) {
		case 0:
		    SET_COORD3( tri_vert[num_vert], l_fx, l_fy, ll_fz ) ;
		    break ;
		case 1:
		    SET_COORD3( tri_vert[num_vert], r_fx, l_fy, lr_fz ) ;
		    break ;
		case 2:
		    SET_COORD3( tri_vert[num_vert], r_fx, u_fy, ur_fz ) ;
		    break ;
		case 3:
		    SET_COORD3( tri_vert[num_vert], l_fx, u_fy, ul_fz ) ;
		    break ;
		}
	    }
	    lib_output_polygon( 3, tri_vert ) ;
	}
    }

    else {
	/* subdivide edges and move in z direction */
	half_width = width>>1 ;
	hside_length = (double)half_width * WIDTH / fnum_pts ;
	rise_height = hside_length * Roughness ;

	/* for each midpoint, find z */
	iz = MOUNTAIN_NO + hash_rand( ll_x + half_width, ll_y, SizeFactor ) ;
	lower_fz = ( ll_fz + lr_fz ) / 2.0 +
					rise_height * lib_gauss_rand( iz ) ;
	iz = MOUNTAIN_NO +
		hash_rand( ll_x + width, ll_y + half_width, SizeFactor ) ;
	right_fz = ( lr_fz + ur_fz ) / 2.0 +
					rise_height * lib_gauss_rand( iz ) ;
	iz = MOUNTAIN_NO +
		hash_rand( ll_x + half_width, ll_y + width, SizeFactor ) ;
	upper_fz = ( ur_fz + ul_fz ) / 2.0 +
					rise_height * lib_gauss_rand( iz ) ;
	iz = MOUNTAIN_NO + hash_rand( ll_x, ll_y + half_width, SizeFactor ) ;
	left_fz  = ( ul_fz + ll_fz ) / 2.0 +
					rise_height * lib_gauss_rand( iz ) ;
	iz = MOUNTAIN_NO +
	    hash_rand( ll_x + half_width, ll_y + half_width, SizeFactor ) ;
	middle_fz = ( ll_fz + lr_fz + ur_fz + ul_fz ) / 4.0 +
			    1.4142136 * rise_height * lib_gauss_rand( iz ) ;

	/* check subsections for subdivision or output */
	grow_mountain( fnum_pts, half_width, ll_x, ll_y,
				    ll_fz, lower_fz, middle_fz, left_fz ) ;
	grow_mountain( fnum_pts, half_width, ll_x+half_width, ll_y,
				    lower_fz, lr_fz, right_fz, middle_fz ) ;
	grow_mountain( fnum_pts, half_width, ll_x+half_width, ll_y+half_width,
				    middle_fz, right_fz, ur_fz, upper_fz ) ;
	grow_mountain( fnum_pts, half_width, ll_x, ll_y+half_width,
				    left_fz, middle_fz, upper_fz, ul_fz ) ;
    }
}
@EOF

chmod 644 mount.c

echo x - rings.c
cat >rings.c <<'@EOF'
/*
 * rings.c - Create objects with 6 pentagonal rings which connect the midpoints
 *	of the edges of a dodecahedron.  A pyramid of these objects is formed,
 *	which the viewer looks upon from the point.  A plane is placed behind
 *	the pyramid for shadows.  Three light sources.
 *
 * Author:  Eric Haines, 3D/Eye, Inc.
 *
 * SizeFactor determines the number of objects output.
 *	Each object has 30 cylinders and 30 spheres.
 *	Total objects = SF*SF + (SF-1)*(SF-1) + ... + 1 plus 1 backdrop square.
 *	formula for # of spheres or cylinders = 5*SF*(SF+1)*(2*SF+1)
 *
 *	SizeFactor	# spheres	# cylinders	# squares
 *	     1		    30		      30	     1
 *	     2		   150		     150	     1
 *	     3		   420		     420	     1
 *
 *	     7		  4200		    4200	     1
 */

#include <stdio.h>
#include <math.h>
#include "def.h"
#include "lib.h"

static	int	SizeFactor = 7 ;

#define	OUTPUT_FORMAT		OUTPUT_CURVES

/* if spread out is > 1, succeeding layers spread out more */
#define	SPREAD_OUT		1

main(argc,argv)
int	argc ;
char	*argv[] ;
{
int	prev_elem, num_elem, num_depth, num_objx, num_objz ;
double	radius, spread, y_diff, xz_diff ;
COORD4	base_pt, apex_pt ;
COORD3	from, at, up ;
COORD3	wvec, light ;
COORD3	back_color, ring_color[6] ;
COORD3	wall[4], offset, dodec[30] ;

    if ( !lib_get_size( argc, argv, &SizeFactor ) ) {
	fprintf( stderr, "usage: %s [size]\n", *argv ) ;
	exit(EXIT_FAIL) ;
    }

    radius = 0.07412 ;	/* cone and sphere radius */

    /* calculate spread of objects */
    spread = 1 / sin( (double)( PI/8.0 ) ) ;
    if ( SPREAD_OUT <= spread ) {
	y_diff = spread / SPREAD_OUT ;
	xz_diff = 1.0 ;
    }
    else {
	y_diff = 1.0 ;
	xz_diff = SPREAD_OUT / spread ;
    }

    /* output viewpoint */
    SET_COORD3( from, -1.0, -spread, 0.5 ) ;
    SET_COORD3( at, from[X], from[Y] + 1.0, from[Z] ) ;
    SET_COORD3( up, 0.0, 0.0, 1.0 ) ;
    lib_output_viewpoint( from, at, up, 45.0, 1.0, 512, 512 ) ;

    /* output background color - UNC sky blue */
    /* note that the background color should never be seen */
    SET_COORD3( back_color, 0.078, 0.361, 0.753 ) ;
    lib_output_background_color( back_color ) ;

    /* output light source */
    SET_COORD3( light, 3.0, -spread, 3.0 ) ;
    lib_output_light( light ) ;
    SET_COORD3( light, -4.0, -spread, 1.0 ) ;
    lib_output_light( light ) ;
    SET_COORD3( light, 2.0, -spread, -4.0 ) ;
    lib_output_light( light ) ;

    /* output wall polygon - white */
    SET_COORD3( back_color, 1.0, 1.0, 1.0 ) ;
    lib_output_color( back_color, 1.0, 0.0, 0.0, 0.0, 0.0 ) ;
    /* just spans 45 degree view + 1% */
    wvec[Y] = y_diff * ( SizeFactor + 1 ) ;
    wvec[X] = wvec[Z] = 1.01 * ( wvec[Y] - from[Y] ) * tan( PI / 8.0 ) ;
    SET_COORD3( wall[0],  wvec[X]+from[X], wvec[Y],  wvec[Z]+from[Z] ) ;
    SET_COORD3( wall[1], -wvec[X]+from[X], wvec[Y],  wvec[Z]+from[Z] ) ;
    SET_COORD3( wall[2], -wvec[X]+from[X], wvec[Y], -wvec[Z]+from[Z] ) ;
    SET_COORD3( wall[3],  wvec[X]+from[X], wvec[Y], -wvec[Z]+from[Z] ) ;
    lib_output_polygon( 4, wall ) ;

    /* set up ring colors - RGB and complements */
    SET_COORD3( ring_color[0], 1.0, 0.0, 0.0 ) ;
    SET_COORD3( ring_color[1], 0.0, 1.0, 0.0 ) ;
    SET_COORD3( ring_color[2], 0.0, 0.0, 1.0 ) ;
    SET_COORD3( ring_color[3], 0.0, 1.0, 1.0 ) ;
    SET_COORD3( ring_color[4], 1.0, 0.0, 1.0 ) ;
    SET_COORD3( ring_color[5], 1.0, 1.0, 0.0 ) ;

    create_dodec( radius, dodec ) ;
    /* radius of osculating cylinders and spheres (no derivation given) */
    base_pt[W] = apex_pt[W] = radius ;

    for ( num_depth = 0 ; num_depth < SizeFactor ; ++num_depth ) {
	offset[Y] = y_diff * (double)(num_depth+1) ;
	for ( num_objz = 0 ; num_objz <= num_depth ; ++num_objz ) {
	    offset[Z] = xz_diff * (double)(2*num_objz - num_depth) ;
	    for ( num_objx = 0 ; num_objx <= num_depth ; ++num_objx ) {
		offset[X] = xz_diff * (double)(2*num_objx - num_depth) ;
		for ( num_elem = 0 ; num_elem < 30 ; ++num_elem ) {
		    COPY_COORD3( base_pt, dodec[num_elem] ) ;
		    ADD2_COORD3( base_pt, offset ) ;
		    if ( num_elem%5 == 0 ) {
			prev_elem = num_elem + 4 ;
			/* new ring beginning - output color */
			lib_output_color( ring_color[num_elem/5]
					, 0.5
					, 0.5, 3.0
					, 0.0, 0.0
					) ;
		    }
		    else {
			prev_elem = num_elem - 1 ;
		    }
		    COPY_COORD3( apex_pt, dodec[prev_elem] ) ;
		    ADD2_COORD3( apex_pt, offset ) ;

		    lib_output_cylcone( base_pt, apex_pt, OUTPUT_FORMAT ) ;
		    lib_output_sphere( base_pt, OUTPUT_FORMAT ) ;
		}
	    }
	}
    }
    exit(EXIT_SUCCESS) ;
}

/* Create the set of 30 points needed to generate the rings */
create_dodec( minor_radius, vertex )
double	minor_radius ;
COORD3	vertex[30] ;
{
int	num_vertex, num_pentagon ;
double	scale, x_rotation, z_rotation ;
COORD3	temp_vertex ;
MATRIX	x_matrix, z_matrix ;

    /* scale object to fit in a sphere of radius 1 */
    scale = 1.0 / ( 1.0 + minor_radius ) ;

    /*
     * define one pentagon as on the XY plane, with points starting along +X
     * and N fifths of the way around the Z axis.
     */
    for ( num_vertex = 0 ; num_vertex < 5 ; ++num_vertex ) {
	vertex[num_vertex][X] = scale * cos((double)num_vertex * 2.0*PI/5.0 ) ;
	vertex[num_vertex][Y] = scale * sin((double)num_vertex * 2.0*PI/5.0 ) ;
	vertex[num_vertex][Z] = 0.0 ;
	vertex[num_vertex][W] = 1.0 ;
    }

    /*
     * find the rotation angle (in radians) along the X axis:
     * angle between two adjacent dodecahedron faces.
     */
    x_rotation = 2.0 *
		acos( cos( (double)(PI/3.0) ) / sin( (double)(PI/5.0) ) ) ;
    lib_create_rotate_matrix( x_matrix, X_AXIS, x_rotation ) ;

    /*
     * Find each of the other 5 pentagons:  rotate along the X axis,
     * then rotate on the Z axis.
     */
    for ( num_pentagon = 1 ; num_pentagon < 6 ; ++num_pentagon ) {
	/*
	 * find the rotation angle (in radians) along the Z axis:
	 * 1/10th plus N fifths of the way around * 2 * PI.
	 */
	z_rotation = PI*( 2.0*(double)(num_pentagon-1)+1.0 ) / 5.0 ;
	lib_create_rotate_matrix( z_matrix, Z_AXIS, z_rotation ) ;

	for ( num_vertex = 0 ; num_vertex < 5 ; ++num_vertex ) {

	    lib_transform_point( temp_vertex
			       , vertex[num_vertex]
			       , x_matrix
			       ) ;

	    lib_transform_point( vertex[5*num_pentagon+num_vertex]
			       , temp_vertex
			       , z_matrix
			       ) ;
	}
    }
}
@EOF

chmod 644 rings.c

echo x - teapot.c
cat >teapot.c <<'@EOF'
/*
 * teapot.c - Tessellate the famous teapot into triangular patches, and sit the
 *	object on top of a checkerboard which is meshed to the same degree.  In
 *	other words, the number of squares on the checkerboard is the same as
 *	the amount of meshing for each of the 32 patches on the teapot, with
 *	each mesh quadrilateral on the teapot further split into two triangles.
 *	Two light sources.
 *
 *	Note that the teapot should always be rendered as a double sided object
 *	(since for some patches both sides can be seen).  (See IEEE CG&A
 *	January 1987 for a history of the teapot.  Note that their vertex list
 *	has duplicates, e.g. 93 and 271, as well as unused vertices, e.g.
 *	205,206,216,223).  Also, note that the bottom (the last four patches)
 *	is not flat - blame Frank Crow, not me.
 *
 * Author:  Eric Haines, 3D/Eye, Inc.
 *
 * SizeFactor determines the number of objects output.
 *	Total patches = 32*2*n*n - 8*n     [degenerates are deleted]
 *	Total squares = n*n
 *
 *	SizeFactor	# patches	# squares
 *	     1		    56		     1		[has more degenerates!]
 *	     2		   240		     4
 *	     3		   552		     9
 *
 *	    12		  9120		   144
 */

#include <stdio.h>
#include <math.h>
#include "def.h"
#include "lib.h"

static	int	SizeFactor = 12 ;

/* comment out the next line to generate all patches except the bottom,
 * i.e. the original Newell teapot
 */
#define BOTTOM

#ifdef	BOTTOM
#define	NUM_PATCHES	32
#else
#define	NUM_PATCHES	28
#endif

int	check_for_cusp() ;

main(argc,argv)
int	argc ;
char	*argv[] ;
{
COORD3	back_color ;
COORD3	light ;
COORD3	from, at, up ;

    if ( !lib_get_size( argc, argv, &SizeFactor ) ) {
	fprintf( stderr, "usage: %s [size]\n", *argv ) ;
	exit(EXIT_FAIL) ;
    }

    if ( SizeFactor == 1 ) {
	fprintf( stderr,
	  "warning: a size of 1 is not supported - use at your own risk\n" ) ;
    }

    /* output viewpoint */
    SET_COORD3( from, 4.86, 7.2, 5.4 ) ;
    SET_COORD3( at, 0.0, 0.0, 0.0 ) ;
    SET_COORD3( up, 0.0, 0.0, 1.0 ) ;
    lib_output_viewpoint( from, at, up, 45.0, 1.0, 512, 512 ) ;

    /* output background color - UNC sky blue */
    SET_COORD3( back_color, 0.078, 0.361, 0.753 ) ;
    lib_output_background_color( back_color ) ;

    /* output light sources */
    SET_COORD3( light, -3.1, 9.8, 12.1 ) ;
    lib_output_light( light ) ;
    SET_COORD3( light, 11.3, 5.1, 8.8 ) ;
    lib_output_light( light ) ;

    output_checkerboard() ;
    output_teapot() ;
    exit(EXIT_SUCCESS) ;
}

output_checkerboard()
{
int	sstep, tstep ;
COORD3	vert[4] ;
COORD3	obj_color ;

    SET_COORD3( obj_color, 1.0, 1.0, 1.0 ) ;
    lib_output_color( obj_color, 0.5, 0.5, 3.0, 0.0, 0.0 ) ;
    for ( sstep = 0 ; sstep < SizeFactor ; sstep++ ) {
	for ( tstep = 0 ; tstep < SizeFactor ; tstep++ ) {
	    if ( ( sstep + tstep ) % 2 ) {
		loc_to_square( sstep, tstep, vert ) ;
		lib_output_polygon( 4, vert ) ;
	    }
	}
    }

    SET_COORD3( obj_color, 0.5, 0.5, 0.5 ) ;
    lib_output_color( obj_color, 0.5, 0.5, 3.0, 0.0, 0.0 ) ;
    for ( sstep = 0 ; sstep < SizeFactor ; sstep++ ) {
	for ( tstep = 0 ; tstep < SizeFactor ; tstep++ ) {
	    if ( !(( sstep + tstep ) % 2) ) {
		loc_to_square( sstep, tstep, vert ) ;
		lib_output_polygon( 4, vert ) ;
	    }
	}
    }
}

loc_to_square( sstep, tstep, vert )
int	sstep ;
int	tstep ;
COORD3	vert[4] ;
{
int	i ;

    for ( i = 0 ; i < 4 ; i++ ) {
	/* vertex 0 & 3 are x.low, 1 & 2 are x.high */
	vert[i][X] = 4.0 * ( 2.0 * (double)(sstep + (i%3 ? 1 : 0) ) /
		(double)SizeFactor - 1.0 ) ;

	/* vertex 0 & 1 are y.low, 2 & 3 are y.high */
	vert[i][Y] = 4.0 * ( 2.0 * (double)(tstep + (i/2 ? 1 : 0) ) /
		(double)SizeFactor - 1.0 ) ;

	vert[i][Z] = 0.0 ;
    }
}

static	int	Patches[32][4][4] = {
/* rim */
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
3, 16, 17, 18, 7, 19, 20, 21, 11, 22, 23, 24, 15, 25, 26, 27,
18, 28, 29, 30, 21, 31, 32, 33, 24, 34, 35, 36, 27, 37, 38, 39,
30, 40, 41, 0, 33, 42, 43, 4, 36, 44, 45, 8, 39, 46, 47, 12,
/* body */
12, 13, 14, 15, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59,
15, 25, 26, 27, 51, 60, 61, 62, 55, 63, 64, 65, 59, 66, 67, 68,
27, 37, 38, 39, 62, 69, 70, 71, 65, 72, 73, 74, 68, 75, 76, 77,
39, 46, 47, 12, 71, 78, 79, 48, 74, 80, 81, 52, 77, 82, 83, 56,
56, 57, 58, 59, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95,
59, 66, 67, 68, 87, 96, 97, 98, 91, 99, 100, 101, 95, 102, 103, 104,
68, 75, 76, 77, 98, 105, 106, 107, 101, 108, 109, 110, 104, 111, 112, 113,
77, 82, 83, 56, 107, 114, 115, 84, 110, 116, 117, 88, 113, 118, 119, 92,
/* handle */
120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135,
123, 136, 137, 120, 127, 138, 139, 124, 131, 140, 141, 128, 135, 142, 143, 132,
132, 133, 134, 135, 144, 145, 146, 147, 148, 149, 150, 151, 68, 152, 153, 154,
135, 142, 143, 132, 147, 155, 156, 144, 151, 157, 158, 148, 154, 159, 160, 68,
/* spout */
161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176,
164, 177, 178, 161, 168, 179, 180, 165, 172, 181, 182, 169, 176, 183, 184, 173,
173, 174, 175, 176, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196,
176, 183, 184, 173, 188, 197, 198, 185, 192, 199, 200, 189, 196, 201, 202, 193,
/* lid */
203, 203, 203, 203, 204, 205, 206, 207, 208, 208, 208, 208, 209, 210, 211, 212,
203, 203, 203, 203, 207, 213, 214, 215, 208, 208, 208, 208, 212, 216, 217, 218,
203, 203, 203, 203, 215, 219, 220, 221, 208, 208, 208, 208, 218, 222, 223, 224,
203, 203, 203, 203, 221, 225, 226, 204, 208, 208, 208, 208, 224, 227, 228, 209,
209, 210, 211, 212, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240,
212, 216, 217, 218, 232, 241, 242, 243, 236, 244, 245, 246, 240, 247, 248, 249,
218, 222, 223, 224, 243, 250, 251, 252, 246, 253, 254, 255, 249, 256, 257, 258,
224, 227, 228, 209, 252, 259, 260, 229, 255, 261, 262, 233, 258, 263, 264, 237,
/* bottom */
265, 265, 265, 265, 266, 267, 268, 269, 270, 271, 272, 273, 92, 119, 118, 113,
265, 265, 265, 265, 269, 274, 275, 276, 273, 277, 278, 279, 113, 112, 111, 104,
265, 265, 265, 265, 276, 280, 281, 282, 279, 283, 284, 285, 104, 103, 102, 95,
265, 265, 265, 265, 282, 286, 287, 266, 285, 288, 289, 270, 95, 94, 93, 92
} ;

static	COORD3	Verts[290] = {
 1.4, 0, 2.4,
 1.4, -0.784, 2.4,
 0.784, -1.4, 2.4,
 0, -1.4, 2.4,
 1.3375, 0, 2.53125,
 1.3375, -0.749, 2.53125,
 0.749, -1.3375, 2.53125,
 0, -1.3375, 2.53125,
 1.4375, 0, 2.53125,
 1.4375, -0.805, 2.53125,
 0.805, -1.4375, 2.53125,
 0, -1.4375, 2.53125,
 1.5, 0, 2.4,
 1.5, -0.84, 2.4,
 0.84, -1.5, 2.4,
 0, -1.5, 2.4,
 -0.784, -1.4, 2.4,
 -1.4, -0.784, 2.4,
 -1.4, 0, 2.4,
 -0.749, -1.3375, 2.53125,
 -1.3375, -0.749, 2.53125,
 -1.3375, 0, 2.53125,
 -0.805, -1.4375, 2.53125,
 -1.4375, -0.805, 2.53125,
 -1.4375, 0, 2.53125,
 -0.84, -1.5, 2.4,
 -1.5, -0.84, 2.4,
 -1.5, 0, 2.4,
 -1.4, 0.784, 2.4,
 -0.784, 1.4, 2.4,
 0, 1.4, 2.4,
 -1.3375, 0.749, 2.53125,
 -0.749, 1.3375, 2.53125,
 0, 1.3375, 2.53125,
 -1.4375, 0.805, 2.53125,
 -0.805, 1.4375, 2.53125,
 0, 1.4375, 2.53125,
 -1.5, 0.84, 2.4,
 -0.84, 1.5, 2.4,
 0, 1.5, 2.4,
 0.784, 1.4, 2.4,
 1.4, 0.784, 2.4,
 0.749, 1.3375, 2.53125,
 1.3375, 0.749, 2.53125,
 0.805, 1.4375, 2.53125,
 1.4375, 0.805, 2.53125,
 0.84, 1.5, 2.4,
 1.5, 0.84, 2.4,
 1.75, 0, 1.875,
 1.75, -0.98, 1.875,
 0.98, -1.75, 1.875,
 0, -1.75, 1.875,
 2, 0, 1.35,
 2, -1.12, 1.35,
 1.12, -2, 1.35,
 0, -2, 1.35,
 2, 0, 0.9,
 2, -1.12, 0.9,
 1.12, -2, 0.9,
 0, -2, 0.9,
 -0.98, -1.75, 1.875,
 -1.75, -0.98, 1.875,
 -1.75, 0, 1.875,
 -1.12, -2, 1.35,
 -2, -1.12, 1.35,
 -2, 0, 1.35,
 -1.12, -2, 0.9,
 -2, -1.12, 0.9,
 -2, 0, 0.9,
 -1.75, 0.98, 1.875,
 -0.98, 1.75, 1.875,
 0, 1.75, 1.875,
 -2, 1.12, 1.35,
 -1.12, 2, 1.35,
 0, 2, 1.35,
 -2, 1.12, 0.9,
 -1.12, 2, 0.9,
 0, 2, 0.9,
 0.98, 1.75, 1.875,
 1.75, 0.98, 1.875,
 1.12, 2, 1.35,
 2, 1.12, 1.35,
 1.12, 2, 0.9,
 2, 1.12, 0.9,
 2, 0, 0.45,
 2, -1.12, 0.45,
 1.12, -2, 0.45,
 0, -2, 0.45,
 1.5, 0, 0.225,
 1.5, -0.84, 0.225,
 0.84, -1.5, 0.225,
 0, -1.5, 0.225,
 1.5, 0, 0.15,
 1.5, -0.84, 0.15,
 0.84, -1.5, 0.15,
 0, -1.5, 0.15,
 -1.12, -2, 0.45,
 -2, -1.12, 0.45,
 -2, 0, 0.45,
 -0.84, -1.5, 0.225,
 -1.5, -0.84, 0.225,
 -1.5, 0, 0.225,
 -0.84, -1.5, 0.15,
 -1.5, -0.84, 0.15,
 -1.5, 0, 0.15,
 -2, 1.12, 0.45,
 -1.12, 2, 0.45,
 0, 2, 0.45,
 -1.5, 0.84, 0.225,
 -0.84, 1.5, 0.225,
 0, 1.5, 0.225,
 -1.5, 0.84, 0.15,
 -0.84, 1.5, 0.15,
 0, 1.5, 0.15,
 1.12, 2, 0.45,
 2, 1.12, 0.45,
 0.84, 1.5, 0.225,
 1.5, 0.84, 0.225,
 0.84, 1.5, 0.15,
 1.5, 0.84, 0.15,
 -1.6, 0, 2.025,
 -1.6, -0.3, 2.025,
 -1.5, -0.3, 2.25,
 -1.5, 0, 2.25,
 -2.3, 0, 2.025,
 -2.3, -0.3, 2.025,
 -2.5, -0.3, 2.25,
 -2.5, 0, 2.25,
 -2.7, 0, 2.025,
 -2.7, -0.3, 2.025,
 -3, -0.3, 2.25,
 -3, 0, 2.25,
 -2.7, 0, 1.8,
 -2.7, -0.3, 1.8,
 -3, -0.3, 1.8,
 -3, 0, 1.8,
 -1.5, 0.3, 2.25,
 -1.6, 0.3, 2.025,
 -2.5, 0.3, 2.25,
 -2.3, 0.3, 2.025,
 -3, 0.3, 2.25,
 -2.7, 0.3, 2.025,
 -3, 0.3, 1.8,
 -2.7, 0.3, 1.8,
 -2.7, 0, 1.575,
 -2.7, -0.3, 1.575,
 -3, -0.3, 1.35,
 -3, 0, 1.35,
 -2.5, 0, 1.125,
 -2.5, -0.3, 1.125,
 -2.65, -0.3, 0.9375,
 -2.65, 0, 0.9375,
 -2, -0.3, 0.9,
 -1.9, -0.3, 0.6,
 -1.9, 0, 0.6,
 -3, 0.3, 1.35,
 -2.7, 0.3, 1.575,
 -2.65, 0.3, 0.9375,
 -2.5, 0.3, 1.125,
 -1.9, 0.3, 0.6,
 -2, 0.3, 0.9,
 1.7, 0, 1.425,
 1.7, -0.66, 1.425,
 1.7, -0.66, 0.6,
 1.7, 0, 0.6,
 2.6, 0, 1.425,
 2.6, -0.66, 1.425,
 3.1, -0.66, 0.825,
 3.1, 0, 0.825,
 2.3, 0, 2.1,
 2.3, -0.25, 2.1,
 2.4, -0.25, 2.025,
 2.4, 0, 2.025,
 2.7, 0, 2.4,
 2.7, -0.25, 2.4,
 3.3, -0.25, 2.4,
 3.3, 0, 2.4,
 1.7, 0.66, 0.6,
 1.7, 0.66, 1.425,
 3.1, 0.66, 0.825,
 2.6, 0.66, 1.425,
 2.4, 0.25, 2.025,
 2.3, 0.25, 2.1,
 3.3, 0.25, 2.4,
 2.7, 0.25, 2.4,
 2.8, 0, 2.475,
 2.8, -0.25, 2.475,
 3.525, -0.25, 2.49375,
 3.525, 0, 2.49375,
 2.9, 0, 2.475,
 2.9, -0.15, 2.475,
 3.45, -0.15, 2.5125,
 3.45, 0, 2.5125,
 2.8, 0, 2.4,
 2.8, -0.15, 2.4,
 3.2, -0.15, 2.4,
 3.2, 0, 2.4,
 3.525, 0.25, 2.49375,
 2.8, 0.25, 2.475,
 3.45, 0.15, 2.5125,
 2.9, 0.15, 2.475,
 3.2, 0.15, 2.4,
 2.8, 0.15, 2.4,
 0, 0, 3.15,
 0.8, 0, 3.15,
 0.8, -0.45, 3.15,
 0.45, -0.8, 3.15,
 0, -0.8, 3.15,
 0, 0, 2.85,
 0.2, 0, 2.7,
 0.2, -0.112, 2.7,
 0.112, -0.2, 2.7,
 0, -0.2, 2.7,
 -0.45, -0.8, 3.15,
 -0.8, -0.45, 3.15,
 -0.8, 0, 3.15,
 -0.112, -0.2, 2.7,
 -0.2, -0.112, 2.7,
 -0.2, 0, 2.7,
 -0.8, 0.45, 3.15,
 -0.45, 0.8, 3.15,
 0, 0.8, 3.15,
 -0.2, 0.112, 2.7,
 -0.112, 0.2, 2.7,
 0, 0.2, 2.7,
 0.45, 0.8, 3.15,
 0.8, 0.45, 3.15,
 0.112, 0.2, 2.7,
 0.2, 0.112, 2.7,
 0.4, 0, 2.55,
 0.4, -0.224, 2.55,
 0.224, -0.4, 2.55,
 0, -0.4, 2.55,
 1.3, 0, 2.55,
 1.3, -0.728, 2.55,
 0.728, -1.3, 2.55,
 0, -1.3, 2.55,
 1.3, 0, 2.4,
 1.3, -0.728, 2.4,
 0.728, -1.3, 2.4,
 0, -1.3, 2.4,
 -0.224, -0.4, 2.55,
 -0.4, -0.224, 2.55,
 -0.4, 0, 2.55,
 -0.728, -1.3, 2.55,
 -1.3, -0.728, 2.55,
 -1.3, 0, 2.55,
 -0.728, -1.3, 2.4,
 -1.3, -0.728, 2.4,
 -1.3, 0, 2.4,
 -0.4, 0.224, 2.55,
 -0.224, 0.4, 2.55,
 0, 0.4, 2.55,
 -1.3, 0.728, 2.55,
 -0.728, 1.3, 2.55,
 0, 1.3, 2.55,
 -1.3, 0.728, 2.4,
 -0.728, 1.3, 2.4,
 0, 1.3, 2.4,
 0.224, 0.4, 2.55,
 0.4, 0.224, 2.55,
 0.728, 1.3, 2.55,
 1.3, 0.728, 2.55,
 0.728, 1.3, 2.4,
 1.3, 0.728, 2.4,
 0, 0, 0,
 1.425, 0, 0,
 1.425, 0.798, 0,
 0.798, 1.425, 0,
 0, 1.425, 0,
 1.5, 0, 0.075,
 1.5, 0.84, 0.075,
 0.84, 1.5, 0.075,
 0, 1.5, 0.075,
 -0.798, 1.425, 0,
 -1.425, 0.798, 0,
 -1.425, 0, 0,
 -0.84, 1.5, 0.075,
 -1.5, 0.84, 0.075,
 -1.5, 0, 0.075,
 -1.425, -0.798, 0,
 -0.798, -1.425, 0,
 0, -1.425, 0,
 -1.5, -0.84, 0.075,
 -0.84, -1.5, 0.075,
 0, -1.5, 0.075,
 0.798, -1.425, 0,
 1.425, -0.798, 0,
 0.84, -1.5, 0.075,
 1.5, -0.84, 0.075
} ;

/* Compute points on each spline surface of teapot by brute force.
 * Forward differencing would be faster, but this is compact & simple.
 */
output_teapot()
{
/* bezier form */
static MATRIX ms = { -1.0,  3.0, -3.0,  1.0,
		      3.0, -6.0,  3.0,  0.0,
		     -3.0,  3.0,  0.0,  0.0,
		      1.0,  0.0,  0.0,  0.0 } ;
int	surf, i, r, c, sstep, tstep, num_tri, num_vert, num_tri_vert ;
double	s[3], t[3] ;
COORD3	vert[4], norm[4] ;
COORD3	obj_color ;
MATRIX	mst, g, mgm[3], tmtx ;

    SET_COORD3( obj_color, 1.0, 0.5, 0.1 ) ;
    lib_output_color( obj_color, 0.75, 0.25, 20.0, 0.0, 0.0 ) ;

    lib_transpose_matrix( mst, ms ) ;

    for ( surf = 0 ; surf < NUM_PATCHES ; surf++ ) {

	/* get M * G * M matrix for x,y,z */
	for ( i = 0 ; i < 3 ; i++ ) {
	    /* get control patches */
	    for ( r = 0 ; r < 4 ; r++ ) {
		for ( c = 0 ; c < 4 ; c++ ) {
		    g[r][c] = Verts[Patches[surf][r][c]][i] ;
		}
	    }

	    lib_matrix_multiply( tmtx, ms, g ) ;
	    lib_matrix_multiply( mgm[i], tmtx, mst ) ;
	}

	/* step along, get points, and output */
	for ( sstep = 0 ; sstep < SizeFactor ; sstep++ ) {
	    for ( tstep = 0 ; tstep < SizeFactor ; tstep++ ) {
		for ( num_tri = 0 ; num_tri < 2 ; num_tri++ ) {
		    for ( num_vert = 0 ; num_vert < 3 ; num_vert++ ) {
			num_tri_vert = ( num_vert + num_tri * 2 ) % 4 ;
			/* trickiness: add 1 to sstep if 1 or 2 */
			s[num_vert] = (double)(sstep + (num_tri_vert/2 ? 1:0) )
				    / (double)SizeFactor ;
			/* trickiness: add 1 to tstep if 2 or 3 */
			t[num_vert] = (double)(tstep + (num_tri_vert%3 ? 1:0) )
				    / (double)SizeFactor ;
		    }
		    points_from_basis( 3, s, t, mgm,
			    vert, norm ) ;
		    /* don't output degenerate polygons */
		    if ( check_for_cusp( 3, vert, norm ) ) {
			lib_output_polypatch( 3, vert, norm ) ;
		    }
		}
	    }
	}
    }
}

points_from_basis( tot_vert, s, t, mgm, vert, norm )
int	tot_vert ;
double	s[] ;
double	t[] ;
MATRIX	mgm[3] ;
COORD3	vert[] ;
COORD3	norm[] ;
{
int	i, num_vert, p ;
double	sval, tval, dsval, dtval, sxyz, txyz ;
COORD3	sdir, tdir ;
COORD4	sp, tp, dsp, dtp, tcoord ;

    for ( num_vert = 0 ; num_vert < tot_vert ; num_vert++ ) {

	sxyz = s[num_vert] ;
	txyz = t[num_vert] ;

	/* get power vectors and their derivatives */
	for ( p = 4, sval = tval = 1.0 ; p-- ; ) {
	    sp[p] = sval ;
	    tp[p] = tval ;
	    sval *= sxyz ;
	    tval *= txyz ;

	    if ( p == 3 ) {
		dsp[p] = dtp[p] = 0.0 ;
		dsval = dtval = 1.0 ;
	    } else {
		dsp[p] = dsval * (double)(3-p) ;
		dtp[p] = dtval * (double)(3-p) ;
		dsval *= sxyz ;
		dtval *= txyz ;
	    }
	}

	/* do for x,y,z */
	for ( i = 0 ; i < 3 ; i++ ) {
	    /* multiply power vectors times matrix to get value */
	    lib_transform_coord( tcoord, sp, mgm[i] ) ;
	    vert[num_vert][i] = DOT4( tcoord, tp ) ;

	    /* get s and t tangent vectors */
	    lib_transform_coord( tcoord, dsp, mgm[i] ) ;
	    sdir[i] = DOT4( tcoord, tp ) ;

	    lib_transform_coord( tcoord, sp, mgm[i] ) ;
	    tdir[i] = DOT4( tcoord, dtp ) ;
	}

	/* find normal */
	CROSS( norm[num_vert], tdir, sdir ) ;
	(void)lib_normalize_vector( norm[num_vert] ) ;
    }
}

/* at the center of the lid's handle and at bottom are cusp points -
 * check if normal is (0 0 0), if so, check that polygon is not degenerate.
 * If degenerate, return FALSE, else set normal.
 */
int	check_for_cusp( tot_vert, vert, norm )
int	tot_vert ;
COORD3	vert[] ;
COORD3	norm[] ;
{
int	count, i, nv ;

    for ( count = 0, i = tot_vert ; i-- ; ) {
	/* check if vertex is at cusp */
	if ( IS_VAL_ALMOST_ZERO( vert[i][X], 0.0001 ) &&
	     IS_VAL_ALMOST_ZERO( vert[i][Y], 0.0001 ) ) {
	    count++ ;
	    nv = i ;
	}
    }

    if ( count > 1 ) {
	/* degenerate */
	return( FALSE ) ;
    }
    if ( count == 1 ) {
	/* check if point is somewhere above the middle of the teapot */
	if ( vert[nv][Z] > 1.5 ) {
	    /* cusp at lid */
	    SET_COORD3( norm[nv], 0.0, 0.0, 1.0 ) ;
	} else {
	    /* cusp at bottom */
	    SET_COORD3( norm[nv], 0.0, 0.0, -1.0 ) ;
	}
    }
    return( TRUE ) ;
}
@EOF

chmod 644 teapot.c

echo x - tetra.c
cat >tetra.c <<'@EOF'
/*
 * tetra.c - Create a tetrahedral pyramid.  This environment is based on the
 *	scene used by Glassner ("Space Subdivision for Fast Ray Tracing," IEEE
 *	CG&A, October 1984) and Kay & Kajiya ("Ray Tracing Complex Scenes,"
 *	SIGGRAPH '86 Proceedings) for testing their ray tracers.
 *	One light source.
 *
 * Author:  Eric Haines, 3D/Eye, Inc.
 *
 * Note:  the view and light positions are the same (after transformation to
 *	a different set of world coordinates) as used by Kay & Kajiya,
 *	courtesy of Tim Kay.  For some reason, the number of shadow rays
 *	generated is different (Kay gets 34K, I get 46K).  One light source.
 *
 * SizeFactor determines the number of polygons output.
 *	Total triangular polygons = 4**SF
 *
 *	SizeFactor	# triangles
 *	     1		     4
 *	     2		    16
 *	     3		    64
 *
 *	     6		  4096
 */

#include <stdio.h>
#include <math.h>
#include "def.h"
#include "lib.h"

static	int	SizeFactor = 6 ;

main(argc,argv)
int	argc ;
char	*argv[] ;
{
COORD3  back_color, tetra_color, light ;
COORD3  from, at, up ;
COORD4  center_pt ;

    if ( !lib_get_size( argc, argv, &SizeFactor ) ) {
	fprintf( stderr, "usage: %s [size]\n", *argv ) ;
	exit(EXIT_FAIL) ;
    }

    /* output viewpoint */
    SET_COORD3( from, 1.022846, -3.177154, -2.174512 ) ;
    SET_COORD3( at, -0.004103, -0.004103, 0.216539 ) ;
    SET_COORD3( up, -0.816497, -0.816497, 0.816497 ) ;
    lib_output_viewpoint( from, at, up, 45.0, 1.0, 512, 512 ) ;

    /* output background color - UNC sky blue */
    SET_COORD3( back_color, 0.078, 0.361, 0.753 ) ;
    lib_output_background_color( back_color ) ;

    /* output light source */
    SET_COORD3( light, 1.876066, -18.123936, -5.000422 ) ;
    lib_output_light( light ) ;

    /* output tetrahedron color - red */
    SET_COORD3( tetra_color, 1.0, 0.2, 0.2 ) ;
    lib_output_color( tetra_color, 1.0, 0.0, 0.0, 0.0, 0.0 ) ;

    /* compute and output tetrahedral object */
    SET_COORD4( center_pt, 0.0, 0.0, 0.0, 1.0 ) ;
    create_tetra( SizeFactor, center_pt ) ;

    exit(EXIT_SUCCESS) ;
}


/* Create tetrahedrons recursively */
create_tetra( depth, center )
int	depth ;
COORD4	center ;
{
int	num_face, num_vert, swap, vert_ord[3] ;
int	x_dir, y_dir, z_dir ;
COORD3	face_pt[3], obj_pt[4] ;
COORD4	sub_center ;

    if ( depth <= 1 ) {
	/* Output tetrahedron */

	/* find opposite corners of a cube which form a tetrahedron */
	for ( num_vert = 0, x_dir = -1 ; x_dir <= 1 ; x_dir += 2 ) {
	    for ( y_dir = -1 ; y_dir <= 1 ; y_dir += 2 ) {
		for ( z_dir = -1 ; z_dir <= 1 ; z_dir += 2 ) {
		    if ( x_dir*y_dir*z_dir == 1 ) {
			obj_pt[num_vert][X] =
					center[X] + (double)x_dir * center[W] ;
			obj_pt[num_vert][Y] =
					center[Y] + (double)y_dir * center[W] ;
			obj_pt[num_vert][Z] =
					center[Z] + (double)z_dir * center[W] ;
			++num_vert ;
		    }
		}
	    }
	}

	/* find faces and output */
	for ( num_face = 0 ; num_face < 4 ; ++num_face ) {
	    /* output order:
	     *   face 0:  points 0 1 2
	     *   face 1:  points 3 2 1
	     *   face 2:  points 2 3 0
	     *   face 3:  points 1 0 3
	     */
	    for ( num_vert = 0 ; num_vert < 3 ; ++num_vert ) {
		vert_ord[num_vert] = (num_face + num_vert) % 4 ;
	    }
	    if ( num_face%2 == 1 ) {
		swap = vert_ord[0] ;
		vert_ord[0] = vert_ord[2] ;
		vert_ord[2] = swap ;
	    }

	    for ( num_vert = 0 ; num_vert < 3 ; ++num_vert ) {
		COPY_COORD3( face_pt[num_vert], obj_pt[vert_ord[num_vert]] ) ;
	    }
	    lib_output_polygon( 3, face_pt ) ;
	}
    }

    else {
	/* Create sub-tetrahedra */

	/* find opposite corners of a cube to form sub-tetrahedra */
	for ( x_dir = -1 ; x_dir <= 1 ; x_dir += 2 ) {
	    for ( y_dir = -1 ; y_dir <= 1 ; y_dir += 2 ) {
		for ( z_dir = -1 ; z_dir <= 1 ; z_dir += 2 ) {
		    if ( x_dir*y_dir*z_dir == 1 ) {
			sub_center[X] =
				center[X] + (double)x_dir * center[W] / 2.0 ;
			sub_center[Y] =
				center[Y] + (double)y_dir * center[W] / 2.0 ;
			sub_center[Z] =
				center[Z] + (double)z_dir * center[W] / 2.0 ;
			sub_center[W] = center[W] / 2.0 ;

			create_tetra( depth-1, sub_center ) ;
		    }
		}
	    }
	}
    }
}
@EOF

chmod 644 tetra.c

echo x - tree.c
cat >tree.c <<'@EOF'
/*
 * tree.c - Creates a tree using Aono & Kunii's generation method.
 *	(See IEEE CG&A May 1984).  A square polygon is placed beneath the
 *	tree to act as a field.  Seven light sources.
 *
 * Author:  Eric Haines, 3D/Eye, Inc.
 *
 * SizeFactor determines the number of objects output.
 *	Total objects = 2**(SF+1)-1 cones and spheres + 1 square polygon.
 *
 *	SizeFactor	# spheres	   # cones	# squares
 *	     1		     3		       3	     1
 *	     2		     7		       7	     1
 *	     3		    15		      15	     1
 *
 *	    11		  4095		    4095	     1
 */

#include <stdio.h>
#include <math.h>
#include "def.h"
#include "lib.h"

static	int	SizeFactor = 11 ;

#define	OUTPUT_FORMAT		OUTPUT_CURVES

/* the following affect the shape of the tree */
#define	BR_ANGLE_0		40.0
#define	BR_ANGLE_1		25.0
#define	BR_CONTR_0		0.65
#define	BR_CONTR_1		0.70
#define	BR_DIAMETER		0.67
#define	DIV_ANGLE		140.0
#define	WIDTH_HEIGHTH_RATIO	0.15

static	MATRIX	Rst_mx[2] ;

main(argc,argv)
int	argc ;
char	*argv[] ;
{
COORD3	field[4] ;
COORD3	from, at, up ;
COORD3	light ;
COORD3	back_color, tree_color ;

    if ( !lib_get_size( argc, argv, &SizeFactor ) ) {
	fprintf( stderr, "usage: %s [size]\n", *argv ) ;
	exit(EXIT_FAIL) ;
    }

    /* output viewpoint */
    SET_COORD3( from, 4.5, 0.4, 2.0 ) ;
    SET_COORD3( at, 0.0, 0.0, 1.5 ) ;
    SET_COORD3( up, 0.0, 0.0, 1.0 ) ;
    lib_output_viewpoint( from, at, up, 45.0, 1.0, 512, 512 ) ;

    /* output background color - UNC sky blue */
    SET_COORD3( back_color, 0.078, 0.361, 0.753 ) ;
    lib_output_background_color( back_color ) ;

    /* output light source */
    SET_COORD3( light, -5.0, 5.0, 50.0 ) ;
    lib_output_light( light ) ;
    SET_COORD3( light, 30.0, -30.0, 30.0 ) ;
    lib_output_light( light ) ;
    SET_COORD3( light, -40.0, -30.0, 20.0 ) ;
    lib_output_light( light ) ;
    SET_COORD3( light, 10.0, 30.0, 40.0 ) ;
    lib_output_light( light ) ;
    SET_COORD3( light, -30.0, 40.0, 10.0 ) ;
    lib_output_light( light ) ;
    SET_COORD3( light, 50.0, 25.0, 20.0 ) ;
    lib_output_light( light ) ;
    SET_COORD3( light, -10.0, -60.0, 30.0 ) ;
    lib_output_light( light ) ;

    /* output field polygon - green */
    SET_COORD3( back_color, 0.2, 0.7, 0.2 ) ;
    lib_output_color( back_color, 1.0, 0.0, 0.0, 0.0, 0.0 ) ;
    SET_COORD3( field[0],  50.0,  50.0, 0.0 ) ;
    SET_COORD3( field[1], -50.0,  50.0, 0.0 ) ;
    SET_COORD3( field[2], -50.0, -50.0, 0.0 ) ;
    SET_COORD3( field[3],  50.0, -50.0, 0.0 ) ;
    lib_output_polygon( 4, field ) ;

    /* set up tree color - brown */
    SET_COORD3( tree_color, 0.55, 0.4, 0.2 ) ;
    lib_output_color( tree_color, 1.0, 0.0, 0.0, 0.0, 0.0 ) ;

    /* create tree */
    create_tree() ;

    exit(EXIT_SUCCESS) ;
}

/*
 * Set up matrices for growth of each branch with respect to the
 * parent branch, then grow each branch.
 */
create_tree()
{
int	i ;
double	branch_angle, branch_contraction, divergence ;
MATRIX	ident_mx, temp1_mx, temp2_mx, tempr_mx, tempst_mx ;

    for ( i = 0 ; i < 2 ; ++i ) {
	if ( i == 0 ) {
	    branch_angle = BR_ANGLE_0 ;
	    divergence = 90.0 ;
	    branch_contraction = BR_CONTR_0 ;
	}
	else {
	    branch_angle = BR_ANGLE_1 ;
	    divergence = DIV_ANGLE + 90.0 ;
	    branch_contraction = BR_CONTR_1 ;
	}

	/* rotate along X axis by branching angle */
	lib_create_rotate_matrix( temp1_mx, X_AXIS, branch_angle*PI/180.0 ) ;

	/* rotate along Z axis by divergence angle */
	lib_create_rotate_matrix( temp2_mx, Z_AXIS, divergence*PI/180.0 ) ;

	lib_matrix_multiply( tempr_mx, temp1_mx, temp2_mx ) ;

	/* include translation of branch, scaled */
	lib_create_identity_matrix( tempst_mx ) ;
	tempst_mx[0][0] = tempst_mx[1][1] = tempst_mx[2][2] =
							branch_contraction ;
	tempst_mx[3][2] = 1.0 ;

	/* concatenate */
	lib_matrix_multiply( Rst_mx[i], tempr_mx, tempst_mx ) ;
    }

    /* set up initial matrix */
    lib_create_identity_matrix( ident_mx ) ;
    grow_tree( ident_mx, 1.0, SizeFactor ) ;
}

/* grow tree branches recursively */
grow_tree( cur_mx, scale, depth )
MATRIX	cur_mx ;
double	scale ;
int	depth ;
{
int	i ;
COORD3	vec ;
COORD4	apex, base ;
MATRIX	new_mx ;

    /* output branch */
    SET_COORD3( vec, 0.0, 0.0, 0.0 ) ;
    lib_transform_point( base, vec, cur_mx ) ;
    base[W] = scale * WIDTH_HEIGHTH_RATIO ;

    SET_COORD3( vec, 0.0, 0.0, 1.0 ) ;
    lib_transform_point( apex, vec, cur_mx ) ;
    apex[W] = base[W] * BR_DIAMETER ;

    lib_output_cylcone( base, apex, OUTPUT_FORMAT ) ;
    lib_output_sphere( apex, OUTPUT_FORMAT ) ;

    if ( depth > 0 ) {
	--depth ;

	for ( i = 0 ; i < 2 ; ++i ) {
	    lib_matrix_multiply( new_mx, Rst_mx[i], cur_mx ) ;
	    grow_tree( new_mx, scale * BR_DIAMETER, depth ) ;
	}
    }
}
@EOF

chmod 644 tree.c

exit 0
