/*
 * lib.c - a library of vector operations, a random number generator, and
 *     object output routines.
 *
 * Version:  3.1a
 * Author:  Eric Haines, 3D/Eye, Inc.
 *
 * Modified: 1 October 1992
 *           Alexander R. Enzmann
 *
 *           I made quite a few changes in order to support multiple raytracers,
 *           with the hopes that this library would become even more useful than
 *           it already is.
 *
 * Modified: 17 Jan 1993
 *           Eduard [esp] Schwan
 *           Removed unused local variables & returned value in lib_output_color
 *
 * Modified: 17 Mar 1993
 *           Eduard [esp] Schwan
 *           Changed POV-Ray refs to OUTPUT_POVRAY_10 & OUTPUT_POVRAY_20
 *           Changed POV-Ray 1.5 refs to 2.0
 *           Passed bg_color to display routines
 *             (unfortunate side-effect is you should now call
 *             lib_output_background_color BEFORE lib_output_viewpoint.
 *             This may not be the best approach - please review!)
 *           Added OUTPUT_RAWTRI output, for creating Raw output for Raw2POV.
 *           Added separator comment lines between fns for reaability
 *
 * Modified: 30 Mar 1993
 *           Eduard [esp] Schwan
 *           Made poly static arrays dynamic
 *           Created lib_init/lib_shutdown routines
 *           Removed unused local variables
 *           Added creation of enclosing sphere of bg color for POV-Ray
 *
 * Modified: 12 April 1993
 *           Alexander R. Enzmann
 *           Updated to include Eric Haines' version 3.1 SPD data types
 *
 * Modified: 16 April 1993
 *           Eduard [esp] Schwan
 *           Sprinkled MULTITASK macros in loops to let cooperative
 *           multitasking OS systems (Macs) breathe, and added a
 *           string version # & accessor routine for reporting.
 *
 * Modified: 27 July 1993
 *           David Hook
 *           Added converter for "art" ray tracer.
 *
 * Modified: 5 August 1993
 *           Eduard [esp] Schwan
 *           Fixed filename bug in lib_get_opts,
 *           cleaned up Mac #pragma statements, made USAGE a fn,
 *           added '/' as command line switch prefix for MSDOS dudes.
 *           Added POV 2.0 cone and disc primitive syntax, removed
 *           redundant "object { }" wrapper around POV 2.0 objects,
 *           added POV 2.0 "background {}" statement.  Also cleaned
 *           up indenting by adding tab_() routines.  Added A Enzmann's
 *           fixes to the POV 1.0 cylinder/cone routines.
 *
 * Modified: 31 August 1993
 *           Eric Haines
 *           Minor syntax things (hey, I like 4 space indents).
 *
 */

/*-----------------------------------------------------------------*/
/* include section */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "lib.h"
#include "drv.h"


/*-----------------------------------------------------------------*/
/* defines/constants section */

#define VBUFFER_SIZE	1024
#define POLYEND_SIZE	512

#define VERT(i, a) ((verts[vbuffer[i]])[a])

/*-----------------------------------------------------------------*/
/* Here are some local variables that are used to control things like
   the current output file, current texture, ... */
static FILE *outfile      = stdout;
static char *texture_name = NULL;
static int  texture_count = 0;
static double texture_ior = 1.0;
static int  format        = OUTPUT_NFF;
static int  u_resolution  = OUTPUT_RESOLUTION;
static int  v_resolution  = OUTPUT_RESOLUTION;
static COORD3 bk_color = {0.0, 0.0, 0.0};
static COORD3 fg_color = {0.0, 0.0, 0.0};
static double view_bounds[2][3];
static int view_init_flag = 0;
static char *lib_version_str = LIB_VERSION;

surface_ptr lib_surfaces = NULL;
object_ptr lib_objects = NULL;
light_ptr lib_lights = NULL;
viewpoint view = {
	{0, 0, -10},
	{0, 0, 0},
	{0, 1, 0},
	45, 1, 1.0e-3, 10, 128, 128,
	{ {1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1} } };

/* Polygon stack for making PLG files */
static object_ptr polygon_stack = NULL;

/* Storage for polygon indices */
static unsigned int *vbuffer = NULL;
static int *poly_end = NULL;

/* Globals to determine which axes can be used to split the polygon */
static int axis1 = 0;
static int axis2 = 1;

/* Globals for tracking indentation level of output file */
static int	tab_width = 4;
static int	tab_level = 0;


/*-----------------------------------------------------------------*/
static void show_usage PARAMS((void))
{
    /* localize the usage strings to be expanded only once for space savings */
#if defined(applec) || defined(THINK_C)
    /* and don't write to stdout on Macs, which don't have console I/O, and  */
    /* won't ever get this error anyway, since parms are auto-generated.     */
#else

    fprintf(stderr, "usage [-s size] [-r format] [-c|t]\n");
    fprintf(stderr, "-s size - input size of database\n");
    fprintf(stderr, "-r format - input database format to output\n");
    fprintf(stderr, "	0   Output direct to the screen (sys dependent)\n");
    fprintf(stderr, "	1   NFF - MTV\n");
    fprintf(stderr, "	2   POV-Ray 1.0\n");
    fprintf(stderr, "	3   Polyray v1.4, v1.5\n");
    fprintf(stderr, "	4   Vivid 2.0\n");
    fprintf(stderr, "	5   QRT 1.5\n");
    fprintf(stderr, "	6   Rayshade\n");
    fprintf(stderr, "	7   POV-Ray 2.0 (format is subject to change)\n");
    fprintf(stderr, "	8   RTrace 8.0.0\n");
    fprintf(stderr, "	9   PLG format for use with rend386\n");
    fprintf(stderr, "	10  Raw triangle output\n");
    fprintf(stderr, "	11  art 2.3\n");
    fprintf(stderr, "-c - output true curved descriptions\n");
    fprintf(stderr, "-t - output tessellated triangle descriptions\n");

#endif
} /* show_usage */


/*-----------------------------------------------------------------*/
/*
 * Command line option parser
 *
 * -s size - input size of database
 * -r format - input database format to output
 *	0   Output direct to the screen (sys dependent)
 *	1   NFF - MTV
 *	2   POV-Ray 1.0
 *	3   Polyray v1.4, v1.5
 *	4   Vivid 2.0
 *	5   QRT 1.5
 *	6   Rayshade
 *	7   POV-Ray 2.0 (format is subject to change)
 *	8   RTrace 8.0.0
 *	9   PLG format for use with "rend386"
 *	10  Raw triangle output
 *	11  art 2.3
 * -c - output true curved descriptions
 * -t - output tessellated triangle descriptions
 * -f file - DXF file to convert/display (used for showdxf.c only)
 *
 * TRUE returned if bad command line detected
 * some of these are useless for the various routines - we're being a bit
 * lazy here...
 */
int	lib_get_opts( argc, argv, p_size, p_rdr, p_curve )
int	argc ;
char	*argv[] ;
int	*p_size, *p_rdr, *p_curve ;
{
int num_arg ;
int val ;

    num_arg = 0 ;

    while ( ++num_arg < argc ) {
	if ( (*argv[num_arg] == '-') || (*argv[num_arg] == '/') ) {
	    switch( argv[num_arg][1] ) {
		case 's':	/* size selection */
		    if ( ++num_arg < argc ) {
			sscanf( argv[num_arg], "%d", &val ) ;
			if ( val < 1 ) {
			    fprintf( stderr,
				    "bad size value %d given\n",val);
			    show_usage();
			    return( TRUE ) ;
			}
			*p_size = val ;
		    } else {
			fprintf( stderr, "not enough args for -s option\n" ) ;
			show_usage();
			return( TRUE ) ;
		    }
		    break ;
		case 'c':	/* true curve output */
		    *p_curve = OUTPUT_CURVES ;
		    break ;
		case 't':	/* tessellated curve output */
		    *p_curve = OUTPUT_PATCHES ;
		    break ;
		case 'r':	/* renderer selection */
		    if ( ++num_arg < argc ) {
			sscanf( argv[num_arg], "%d", &val ) ;
			if ( val < OUTPUT_VIDEO || val >= OUTPUT_DELAYED ) {
			    fprintf( stderr,
				    "bad renderer value %d given\n",val);
			    show_usage();
			    return( TRUE ) ;
			}
			*p_rdr = val ;
		    } else {
			fprintf( stderr, "not enough args for -r option\n" ) ;
			show_usage();
			return( TRUE ) ;
		    }
		    break ;
		default:
		    fprintf( stderr, "unknown argument -%c\n",
			    argv[num_arg][1] ) ;
		    show_usage();
		    return( TRUE ) ;
	    }
	} else {
	    fprintf( stderr, "unknown argument %s\n",
		    argv[num_arg] ) ;
	    show_usage();
	    return( TRUE ) ;
	}
    }
    return( FALSE ) ;
}


/*-----------------------------------------------------------------*/
static void tab_indent PARAMS((void))
{
    int	k;
    /* Q&D way to do it... */
    for (k=0; k<tab_width*tab_level; k++)
	putc(' ', outfile);
} /* tab_printf */


/*-----------------------------------------------------------------*/
static void tab_inc PARAMS((void))
{
    tab_level++;
} /* tab_inc */


/*-----------------------------------------------------------------*/
static void tab_dec PARAMS((void))
{
    tab_level--;
    if (tab_level < 0)
	tab_level = 0;
} /* tab_dec */


/*-----------------------------------------------------------------*/
/*
 * Library initialization/teardown functions
 */
/*-----------------------------------------------------------------*/

static void
lib_initialize PARAMS((void))
{
    vbuffer = (unsigned int*)malloc(VBUFFER_SIZE * sizeof(unsigned int));
    poly_end = (int*)malloc(POLYEND_SIZE * sizeof(int));
    if (!vbuffer || !poly_end) {
	fprintf(stderr, "Error(lib_initialize): Can't allocate memory.\n");
	exit(1);
    }
} /* lib_initialize */


/*-----------------------------------------------------------------*/
static void
lib_shutdown PARAMS((void))
{
    if (vbuffer) {
	free(vbuffer);
	vbuffer = NULL;
    }
    if (poly_end) {
	free(poly_end);
	poly_end = NULL;
    }
} /* lib_shutdown */


/* Library info functions */

/*-----------------------------------------------------------------*/
char *
lib_get_version_str PARAMS((void))
{
    return lib_version_str;
} /* lib_get_version_str */


/*-----------------------------------------------------------------*/
/*
 * Routines to set/reset the various output parameters
 */
/*-----------------------------------------------------------------*/
void
lib_set_output_file(new_outfile)
    FILE *new_outfile;
{
    if (new_outfile == NULL)
	outfile = stdout;
    else
	outfile = new_outfile;
}

/*-----------------------------------------------------------------*/
void
lib_set_default_texture(default_texture)
    char *default_texture;
{
    texture_name = default_texture;
}

/*-----------------------------------------------------------------*/
void
lib_set_raytracer(default_tracer)
    int default_tracer;
{
    if (default_tracer < OUTPUT_VIDEO ||
        default_tracer > OUTPUT_DELAYED) {
	fprintf(stderr, "Unknown renderer index: %d\n", default_tracer);
	exit(1);
    }
    format = default_tracer;
}

/*-----------------------------------------------------------------*/
void
lib_set_polygonalization(u_steps, v_steps)
    int u_steps, v_steps;
{
    if (u_steps > 0 && v_steps > 0) {
	u_resolution = u_steps;
	v_resolution = v_steps;
    }
}

/*-----------------------------------------------------------------*/
static void
lookup_surface_stats(index, tcount, tior)
    int index, *tcount;
    double *tior;
{
    surface_ptr temp_ptr = lib_surfaces;

    while (temp_ptr != NULL && temp_ptr->surf_index != index)
	temp_ptr = temp_ptr->next;
    *tior = temp_ptr->ior;
    if (*tior < 1.0) *tior = 1.0;
    *tcount = temp_ptr->surf_index;
}

/* OUTPUT ROUTINES */

/*-----------------------------------------------------------------*/
/*
 * 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.
 *   Fov:    Vertical field of view of the camera
 *   Aspect: Aspect ratio of horizontal fov to vertical fov
 *   Hither: Minimum distance to any ray-surface intersection
 *   Resx:   X resolution of resulting image
 *   Resy:   Y resolution of resulting image
 *
 * 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."
 */
void
lib_output_viewpoint(from, at, up,
		     fov_angle, aspect_ratio, hither,
		     resx, resy)
    COORD3 from, at, up;
    double fov_angle, aspect_ratio, hither;
    int    resx, resy;
{
    COORD4 viewvec, rightvec;
    double frustrumheight, frustrumwidth;

    switch (format) {
	case OUTPUT_DELAYED:
	case OUTPUT_VIDEO:
	case OUTPUT_PLG:
	    /* Save the various view parameters */
	    COPY_COORD3(view.from, from);
	    COPY_COORD3(view.at, at);
	    COPY_COORD3(view.up, up);
	    view.angle  = fov_angle;
	    view.hither = hither;
	    view.resx   = resx;
	    view.resy   = resy;
	    view.aspect = aspect_ratio;

	    /* Make the 3D clipping box for this view */
	    view_bounds[0][0] = -view.resx/2;
	    view_bounds[1][0] =  view.resx/2;
	    view_bounds[0][1] = -view.resy/2;
	    view_bounds[1][1] =  view.resy/2;
	    view_bounds[0][2] = view.hither;
	    view_bounds[1][2] = 1.0e10;

	    /* Generate the perspective view matrix */
	    lib_create_view_matrix(view.tx, view.from, view.at,
				   view.up, view.resx, view.resy,
				   view.angle, view.aspect);

	    /* Turn on graphics using system dependent video routines */
	    if (format == OUTPUT_VIDEO) {
	        display_init(view.resx, view.resy, bk_color);
	        view_init_flag = 1;
	    }
	    break;

        case OUTPUT_NFF:
	    fprintf(outfile, "v\n");
	    fprintf(outfile, "from %g %g %g\n", from[X], from[Y], from[Z]);
	    fprintf(outfile, "at %g %g %g\n", at[X], at[Y], at[Z]);
	    fprintf(outfile, "up %g %g %g\n", up[X], up[Y], up[Z]);
	    fprintf(outfile, "angle %g\n", fov_angle);
	    fprintf(outfile, "hither %g\n", hither);
	    fprintf(outfile, "resolution %d %d\n", resx, resy);
	    break;

	case OUTPUT_POVRAY_10:
	case OUTPUT_POVRAY_20:
	    /* Let's get a set of vectors that are all at right angles to each
	       other that describe the view given. */
	    lib_normalize_vector(up);
	    SUB3_COORD3(viewvec, at, from);
	    lib_normalize_vector(viewvec);
	    CROSS(rightvec, up, viewvec);
	    lib_normalize_vector(rightvec);
	    CROSS(up, viewvec, rightvec);
	    lib_normalize_vector(up);

	    /* Calculate the height of the view frustrum in world coordinates.
	       and then scale the right and up vectors appropriately. */
	    frustrumheight = 2.0 * tan(PI * fov_angle / 360.0);
	    frustrumwidth = aspect_ratio * frustrumheight;
	    up[X] *= frustrumheight;
	    up[Y] *= frustrumheight;
	    up[Z] *= frustrumheight;
	    rightvec[X] *= frustrumwidth;
	    rightvec[Y] *= frustrumwidth;
	    rightvec[Z] *= frustrumwidth;

	    tab_indent();
	    fprintf(outfile, "camera {\n");
	    tab_inc();

	    if (format == OUTPUT_POVRAY_10) {
		tab_indent();
		fprintf(outfile, "location <%g %g %g>\n",
			from[X], from[Y], from[Z]);
		tab_indent();
		fprintf(outfile, "direction <%g %g %g>\n",
			viewvec[X], viewvec[Y], viewvec[Z]);
		tab_indent();
		fprintf(outfile, "right <%g %g %g>\n",
			rightvec[X], rightvec[Y], rightvec[Z]);
		tab_indent();
		fprintf(outfile, "up <%g %g %g>\n",
			up[X], up[Y], up[Z]);
		}
	    else {
		tab_indent();
		fprintf(outfile, "location <%g, %g, %g>\n",
			from[X], from[Y], from[Z]);
		tab_indent();
		fprintf(outfile, "direction <%g, %g, %g>\n",
			viewvec[X], viewvec[Y], viewvec[Z]);
		tab_indent();
		fprintf(outfile, "right <%g, %g, %g>\n",
			rightvec[X], rightvec[Y], rightvec[Z]);
		tab_indent();
		fprintf(outfile, "up <%g, %g, %g>\n",
			up[X], up[Y], up[Z]);
	    }
	    tab_dec();
	    fprintf(outfile, "} // camera\n");
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_POLYRAY:
	    tab_indent();
	    fprintf(outfile, "viewpoint {\n");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "from <%g, %g, %g>\n", from[X], from[Y], from[Z]);
	    tab_indent();
	    fprintf(outfile, "at <%g, %g, %g>\n", at[X], at[Y], at[Z]);
	    tab_indent();
	    fprintf(outfile, "up <%g, %g, %g>\n", up[X], up[Y], up[Z]);
	    tab_indent();
	    fprintf(outfile, "angle %g\n", fov_angle);
	    tab_indent();
	    fprintf(outfile, "aspect %g\n", aspect_ratio);
	    tab_indent();
	    fprintf(outfile, "hither %g\n", hither);
	    tab_indent();
	    fprintf(outfile, "resolution %d, %d\n", resx, resy);

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "}\n");
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_VIVID:
	    tab_indent();
	    fprintf(outfile, "studio {\n");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "from %g %g %g\n", from[X], from[Y], from[Z]);
	    tab_indent();
	    fprintf(outfile, "at %g %g %g\n", at[X], at[Y], at[Z]);
	    tab_indent();
	    fprintf(outfile, "up %g %g %g\n", up[X], up[Y], up[Z]);
	    tab_indent();
	    fprintf(outfile, "angle %g\n", fov_angle);
	    tab_indent();
	    fprintf(outfile, "aspect %g\n", aspect_ratio);
	    tab_indent();
	    fprintf(outfile, "resolution %d %d\n", resx, resy);
	    tab_indent();
	    fprintf(outfile, "no_exp_trans\n");

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "}\n");
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_QRT:
	    tab_indent();
	    fprintf(outfile, "OBSERVER = (\n");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "loc = (%g,%g,%g),\n", from[X], from[Y], from[Z]);
	    tab_indent();
	    fprintf(outfile, "lookat = (%g,%g,%g),\n", at[X], at[Y], at[Z]);
	    tab_indent();
	    fprintf(outfile, "up = (%g,%g,%g)\n", up[X], up[Y], up[Z]);
	    tab_dec();
	    tab_indent();
	    fprintf(outfile, ")\n");

	    tab_indent();
	    fprintf(outfile, "FOC_LENGTH = %g\n",
		    35.0 / tan(PI * fov_angle / 360.0));
	    tab_indent();
	    fprintf(outfile, "DEFAULT (\n");
	    tab_inc();
	    tab_indent();
	    fprintf(outfile, "aspect = %g,\n", 6.0 * aspect_ratio / 7.0);
	    tab_indent();
	    fprintf(outfile, "x_res = %d, y_res = %d\n", resx, resy);
	    tab_dec();
	    tab_indent();
	    fprintf(outfile, ")\n");

	    /* QRT insists on having the output file as part of the data text */
	    tab_indent();
	    fprintf(outfile, "FILE_NAME = qrt.tga\n");
	    break;

	case OUTPUT_RAYSHADE:
	    fprintf(outfile, "eyep %g %g %g\n", from[X], from[Y], from[Z]);
	    fprintf(outfile, "lookp %g %g %g\n", at[X], at[Y], at[Z]);
	    fprintf(outfile, "up %g %g %g\n", up[X], up[Y], up[Z]);
	    fprintf(outfile, "fov %g %g\n", aspect_ratio * fov_angle,
		     fov_angle);
	    fprintf(outfile, "screen %d %d\n", resx, resy);
	    fprintf(outfile, "sample 1 nojitter\n");
	    break;

	case OUTPUT_RTRACE:
	    fprintf(outfile, "View\n");
	    fprintf(outfile, "%g %g %g\n", from[X], from[Y], from[Z]);
	    fprintf(outfile, "%g %g %g\n", at[X], at[Y], at[Z]);
	    fprintf(outfile, "%g %g %g\n", up[X], up[Y], up[Z]);
	    fprintf(outfile, "%g %g\n", aspect_ratio * fov_angle/2,
		     fov_angle/2);
	    break;

	case OUTPUT_ART:
	    fprintf(outfile, "maxhitlevel 4\n");
	    fprintf(outfile, "screensize 0.0, 0.0\n");
	    fprintf(outfile, "fieldofview %f\n", fov_angle);
	    fprintf(outfile, "up(%f, %f, %f)\n", up[X], up[Y], up[Z]);
	    fprintf(outfile, "lookat(%f, %f, %f, ", from[X], from[Y], from[Z]);
	    fprintf(outfile, "%f, %f, %f, 0.0)\n", at[X], at[Y], at[Z]);
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_RAWTRI:
	    break;
    }
}

/*-----------------------------------------------------------------*/
/*
 * Output light.  A light is defined by position.  All lights have the same
 * intensity.
 *
 */
void
lib_output_light(center_pt)
    COORD4 center_pt;
{
    double lscale;
    light_ptr new_light;

    if (center_pt[W] != 0.0)
	lscale = center_pt[W];
    else
	lscale = 1.0;

    switch (format) {
	case OUTPUT_DELAYED:
	    new_light = (light_ptr)malloc(sizeof(struct light_struct));
	    if (new_light == NULL)
	       /* Quietly fail & return */
	       return;
	    COPY_COORD4(new_light->center_pt, center_pt);
	    new_light->center_pt[W] = lscale;
	    new_light->next = lib_lights;
	    lib_lights = new_light;
	    break;

	case OUTPUT_VIDEO:
	case OUTPUT_PLG:
	    /* Not currently doing anything with lights */
	    break;

	case OUTPUT_NFF:
	    fprintf(outfile, "l %g %g %g\n",
		    center_pt[X], center_pt[Y], center_pt[Z]);
	    break;

	case OUTPUT_POVRAY_10:
	    tab_indent();
	    fprintf(outfile, "object {\n");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "light_source {\n");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "<%g %g %g>",
		    center_pt[X], center_pt[Y], center_pt[Z]);
	    fprintf(outfile, " color red %g green %g blue %g\n",
		    lscale, lscale, lscale);

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "} // light\n");

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "} // object\n");

	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_POVRAY_20:
	    tab_indent();
	    fprintf(outfile, "light_source {\n");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "<%g, %g, %g>",
		    center_pt[X], center_pt[Y], center_pt[Z]);
	    fprintf(outfile, " color red %g green %g blue %g\n",
		    lscale, lscale, lscale);

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "} // light\n");

	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_POLYRAY:
	    tab_indent();
	    fprintf(outfile, "light <%g, %g, %g>, <%g, %g, %g>\n",
		    lscale, lscale, lscale,
		    center_pt[X], center_pt[Y], center_pt[Z]);
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_VIVID:
	    tab_indent();
	    fprintf(outfile, "light {type point position %g %g %g",
		    center_pt[X], center_pt[Y], center_pt[Z]);
	    fprintf(outfile, " color %g %g %g }\n",
		    lscale, lscale, lscale);
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_QRT:
	    tab_indent();
	    fprintf(outfile, "LAMP ( loc = (%g,%g,%g), dist = 0, radius = 1,",
		    center_pt[X], center_pt[Y], center_pt[Z]);
	    fprintf(outfile, " amb = (%g,%g,%g) )\n",
		    lscale, lscale, lscale);
	    break;

	case OUTPUT_RAYSHADE:
	    fprintf(outfile, "light %g point %g %g %g\n",
		    center_pt[W], center_pt[X], center_pt[Y], center_pt[Z]);
	    break;

	case OUTPUT_RTRACE:
	    fprintf(outfile, "1 %g %g %g %g %g %g\n",
	    center_pt[X], center_pt[Y], center_pt[Z],
	    center_pt[W], center_pt[W], center_pt[W]);

	case OUTPUT_ART:
	    tab_indent();
	    fprintf(outfile, "light \n{");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "location(%g, %g, %g)  colour 0.5, 0.5, 0.5\n",
		   center_pt[X], center_pt[Y], center_pt[Z]);

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "}\n");
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_RAWTRI:
	    break;
    }
}

/*-----------------------------------------------------------------*/
/*
 * Output background color.  A color is simply RGB (monitor dependent, but
 * that's life).
 */
void
lib_output_background_color(color)
    COORD3 color;
{
    switch (format) {
	case OUTPUT_VIDEO:
	case OUTPUT_DELAYED:
	case OUTPUT_PLG:
	    COPY_COORD3(bk_color, color);
	    break;

	case OUTPUT_NFF:
	    fprintf(outfile, "b %g %g %g\n", color[X], color[Y], color[Z]);
	    break;

	case OUTPUT_POVRAY_10:
	    /* POV-Ray 1.0 does not support a background color */
	    /* Instead, create arbitrarily large enclosing sphere of that
	     * color */
	    tab_indent();
	    fprintf(outfile, "// background color:\n");

	    tab_indent();
	    fprintf(outfile, "object {\n");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "sphere { <0 0 0> 9000  ");
	    fprintf(outfile,
	    "texture { ambient 1 diffuse 0 color red %g green %g blue %g } }\n",
		       color[X], color[Y], color[Z]);
	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "} // object - background\n");
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_POVRAY_20:
	    tab_indent();
	    fprintf(outfile, "background { color red %g green %g blue %g }\n",
		   color[X], color[Y], color[Z]);
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_POLYRAY:
	    tab_indent();
	    fprintf(outfile, "background <%g, %g, %g>\n",
		    color[X], color[Y], color[Z]);
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_VIVID:
	    /* Vivid insists on putting the background into the studio */
	    tab_indent();
	    fprintf(outfile, "studio { background %g %g %g }\n",
		    color[X], color[Y], color[Z]);
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_QRT:
	    tab_indent();
	    fprintf(outfile, "SKY ( horiz = (%g,%g,%g), zenith = (%g,%g,%g),",
		    color[X], color[Y], color[Z],
		    color[X], color[Y], color[Z]);
	    fprintf(outfile, " dither = 0 )\n");
	    break;

	case OUTPUT_RAYSHADE:
	    fprintf(outfile, "background %g %g %g\n",
		    color[X], color[Y], color[Z]);
	    break;

	case OUTPUT_RTRACE:
	    fprintf(outfile, "Colors\n");
	    fprintf(outfile, "%g %g %g\n", color[X], color[Y], color[Z]);
	    fprintf(outfile, "0 0 0\n");
	    break;

	case OUTPUT_RAWTRI:
	    break;

	case OUTPUT_ART:
	    tab_indent();
	    fprintf(outfile, "background %g, %g, %g\n",
		    color[X], color[Y], color[Z]);
	    fprintf(outfile, "\n");
	    break;
    }
}

/*-----------------------------------------------------------------*/
static char *
create_surface_name(name, val)
    char *name;
    int val;
{
    char *txname;

    if (name != NULL)
	return name;

    txname = (char *)malloc(7*sizeof(char));
    if (txname == NULL)
	return NULL;
    sprintf(txname, "txt%03d", val);
    txname[6] = '\0';
    return txname;
}

/*-----------------------------------------------------------------*/
/*
 * Output color and shading parameters for all following objects
 *
 * For POV-Ray and Polyray, a character string will be returned that
 * identified this texture.  The default texture will be updated with
 * the name generated by this function.
 *
 * Meaning of the color and shading parameters:
 *    name   = name that this surface can be referenced by...
 *    color  = surface color
 *    ka     = ambient component
 *    kd     = diffuse component
 *    ks     = amount contributed from the reflected direction
 *    shine  = contribution from specular highlights
 *    ang    = angle at which the specular highlight falls to 50% of maximum
 *    t      = amount from the refracted direction
 *    i_of_r = index of refraction of the surface
 *
 */
char *
lib_output_color(name, color, ka, kd, ks, shine, ang, kt, i_of_r)
    char *name;
    COORD3 color;
    double ka, kd, ks, shine, ang, kt, i_of_r;
{
    surface_ptr new_surf;
    char *txname = NULL;
    double phong_pow;

    /* Increment the number of surface types we know about */
    ++texture_count;
    texture_ior = i_of_r;

    /* Calculate the Phong coefficient */
    phong_pow = PI * ang / 180.0;
    if (phong_pow <= 0.0)
	phong_pow = 100000.0;
    else if (phong_pow >= (PI/4.0))
	phong_pow = 0.000001;
    else
	phong_pow = -(log(2.0) / log(cos(2.0 * phong_pow)));

    switch (format) {
	case OUTPUT_DELAYED:
	    new_surf = (surface_ptr)malloc(sizeof(struct surface_struct));
	    if (new_surf == NULL)
	       /* Quietly fail */
		return NULL;
	    new_surf->surf_name = create_surface_name(name, texture_count);
	    new_surf->surf_index = texture_count;
	    COPY_COORD3(new_surf->color, color);
	    new_surf->ka = ka;
	    new_surf->kd = kd;
	    new_surf->ks = ks;
	    new_surf->shine = shine;
	    new_surf->ang = ang;
	    new_surf->kt = kt;
	    new_surf->ior = i_of_r;
	    new_surf->next = lib_surfaces;
	    lib_surfaces = new_surf;
	    break;

	case OUTPUT_PLG:
	case OUTPUT_VIDEO:
	    COPY_COORD3(fg_color, color);
	    break;

	case OUTPUT_NFF:
	    fprintf(outfile, "f %g %g %g %g %g %g %g %g\n",
		color[X], color[Y], color[Z], kd, ks, phong_pow, kt, i_of_r);
	    break;

	case OUTPUT_POVRAY_10:
	    txname = create_surface_name(name, texture_count);
	    tab_indent();
	    fprintf(outfile, "#declare %s = texture {\n", txname);
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "color red %g green %g blue %g",
		    color[X], color[Y], color[Z]);
	    if (kt > 0)
		fprintf(outfile, " alpha %g", kt);
	    fprintf(outfile, "\n");

	    tab_indent();
	    fprintf(outfile, "ambient %g\n", ka);

	    tab_indent();
	    fprintf(outfile, "diffuse %g\n", kd);

	    if (shine != 0) {
		tab_indent();
		fprintf(outfile, "phong %g phong_size %g\n", shine, phong_pow);
	    }

	    if (ks != 0) {
		tab_indent();
		fprintf(outfile, "reflection %g\n", ks);
	    }

	    if (kt != 0) {
		tab_indent();
		fprintf(outfile, "refraction 1.0 ior %g\n", i_of_r);
	    }

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "} // texture %s\n", txname);
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_POVRAY_20:
	    txname = create_surface_name(name, texture_count);
	    tab_indent();
	    fprintf(outfile, "#declare %s = texture {\n", txname);
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "pigment {\n");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "color red %g green %g blue %g",
		    color[X], color[Y], color[Z]);
	    if (kt > 0)
		fprintf(outfile, " filter %g", kt);
	    fprintf(outfile, "\n");

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "} // pigment\n");

	    tab_indent();
	    fprintf(outfile, "// normal {\n");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "// bumps, ripples, etc.\n");

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "// } // normal\n");

	    tab_indent();
	    fprintf(outfile, "finish {\n");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "ambient %g\n", ka);

	    tab_indent();
	    fprintf(outfile, "diffuse %g\n", kd);

	    if (shine != 0) {
		tab_indent();
		fprintf(outfile, "phong %g  phong_size %g\n", shine, phong_pow);
	    }

	    if (ks != 0) {
		tab_indent();
		fprintf(outfile, "reflection %g\n", ks);
	    }

	    if (kt != 0) {
		tab_indent();
		fprintf(outfile, "refraction 1.0 ior %g\n", i_of_r);
	    }

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "} // finish\n");

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "} // texture %s\n", txname);
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_POLYRAY:
	    txname = create_surface_name(name, texture_count);
	    tab_indent();
	    fprintf(outfile, "define %s\n", txname);

	    tab_indent();
	    fprintf(outfile, "texture {\n");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "surface {\n");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "ambient <%g, %g, %g>, %g\n",
		   color[X], color[Y], color[Z], ka);

	    tab_indent();
	    fprintf(outfile, "diffuse <%g, %g, %g>, %g\n",
		   color[X], color[Y], color[Z], kd);

	    if (shine != 0) {
		tab_indent();
		fprintf(outfile, "specular white, %g\n", shine);
		tab_indent();
		fprintf(outfile, "microfacet Phong %g\n", ang);
	    }

	    if (ks != 0) {
		tab_indent();
		fprintf(outfile, "reflection white, %g\n", ks);
	    }

	    if (kt != 0) {
		tab_indent();
		fprintf(outfile, "transmission white, %g, %g\n", kt, i_of_r);
	    }

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "}\n");

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "}\n");
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_VIVID:
	    tab_indent();
	    fprintf(outfile, "surface {\n");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "ambient %g %g %g\n",
		    ka * color[X], ka * color[Y], ka * color[Z], ka);

	    tab_indent();
	    fprintf(outfile, "diffuse %g %g %g\n",
		    kd * color[X], kd * color[Y], kd * color[Z], ka);

	    if (shine != 0) {
		tab_indent();
		fprintf(outfile, "shine %g %g %g %g\n",
		       phong_pow, shine, shine, shine);
	    }
	    if (ks != 0) {
		tab_indent();
		fprintf(outfile, "specular %g %g %g\n", ks, ks, ks);
	    }
	    if (kt != 0) {
		tab_indent();
		fprintf(outfile, "transparent %g %g %g\n",
		     kt * color[X], kt * color[Y], kt * color[Z]);
		tab_indent();
		fprintf(outfile, "ior %g\n", i_of_r);
	    }

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, "}\n");
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_QRT:
	    tab_indent();
	    fprintf(outfile, "DEFAULT (\n");
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "amb = (%g,%g,%g),\n",
		    ka * color[X], ka * color[Y], ka * color[Z], ka);
	    tab_indent();
	    fprintf(outfile, "diff = (%g,%g,%g),\n",
		    kd * color[X], kd * color[Y], kd * color[Z], ka);
	    tab_indent();
	    fprintf(outfile, "reflect = %g, sreflect = %g,\n",
		    shine, phong_pow);
	    tab_indent();
	    fprintf(outfile, "mirror = (%g,%g,%g),\n",
		    ks * color[X], ks * color[Y], ks * color[Z]);
	    tab_indent();
	    fprintf(outfile, "trans = (%g,%g,%g), index = %g,\n",
		    kt * color[X], kt * color[Y], kt * color[Z], i_of_r);
	    tab_indent();
	    fprintf(outfile, "dither = 0\n");

	    tab_dec();
	    tab_indent();
	    fprintf(outfile, ")\n");
	    fprintf(outfile, "\n");
	    break;

	case OUTPUT_RAYSHADE:
	    txname = create_surface_name(name, texture_count);
	    tab_indent();
	    fprintf(outfile, "surface %s\n", txname);
	    tab_inc();

	    tab_indent();
	    fprintf(outfile, "ambient %g %g %g\n",
		   ka * color[X], ka * color[Y], ka * color[Z]);
	    tab_indent();
	    fprintf(outfile, "diffuse %g %g %g\n",
		   kd * color[X], kd * color[Y], kd * color[Z]);

	    if (shine != 0) {
		tab_indent();
		fprintf(outfile, "specular %g %g %g\n", shine, shine, shine);
		tab_indent();
		fprintf(outfile, "specpow %g\n", phong_pow);
	    }

	    if (ks != 0) {
		if (shine == 0.0) {
		    /* If there is no Phong highlighting, but there is
		       reflectivity, then we need to define the color of
		       specular reflections */
		    tab_indent();
		    fprintf(outfile, "specular 1.0 1.0 1.0\n");
		    tab_indent();
		    fprintf(outfile, "specpow 0.0\n");
		}
		tab_indent();
		fprintf(outfile, "reflect %g\n", ks);
	    }

	    if (kt != 0) {
		tab_indent();
		fprintf(outfile, "transp %g index %g\n", kt, i_of_r);
	    }

	    tab_dec();
	    tab_indent();
	    break;

	case OUTPUT_RTRACE:
	    if (shine > 0 && ks == 0.0) ks = shine;
	    fprintf(outfile, "1 %g %g %g %g %g %g %g %g %g %g 0 %g %g %g\n",
		    color[X], color[Y], color[Z],
		    kd, kd, kd,
		    ks, ks, ks,
		    (phong_pow > 100.0 ? 100.0 : phong_pow),
		    kt, kt, kt);
	    break;

	case OUTPUT_RAWTRI:
	    txname = create_surface_name(name, texture_count);
	    break;

	case OUTPUT_ART:
	    tab_indent();
	    fprintf(outfile, "colour %f, %f, %f\n",
		    color[X], color[Y], color[Z]);
	    tab_indent();
	    fprintf(outfile, "ambient %f, %f, %f\n",
		    color[X] * 0.05, color[Y] * 0.05, color[Z] * 0.05);

	    if (ks != 0.0) {
		tab_indent();
		fprintf(outfile, "material %f, %f, %f, %f\n",
			    i_of_r, kd, ks, phong_pow);
	    } else {
		tab_indent();
		fprintf(outfile, "material %f, %f, 0.0, 0.0\n",
			    i_of_r, kd);
	    }

	    tab_indent();
	    fprintf(outfile, "reflectance %f\n", ks);
	    tab_indent();
	    fprintf(outfile, "transparency %f\n", kt);
	    fprintf(outfile, "\n");
	    break;
    }

    /* Stash away the current texture name */
    texture_name = txname;

    return txname;
}

/*-----------------------------------------------------------------*/
static void
lib_output_polygon_cylcone(base_pt, apex_pt)
    COORD4 base_pt, apex_pt;
{
    double angle, delta_angle, divisor;
    COORD3 axis, dir, norm_axis, start_norm;
    COORD3 norm[4], vert[4], start_radius[4];
    MATRIX mx;
    int    i;

    SUB3_COORD3(axis, apex_pt, base_pt);
    COPY_COORD3(norm_axis, axis);
    lib_normalize_vector(norm_axis);

    SET_COORD3(dir, 0.0, 0.0, 1.0);
    CROSS(start_norm, axis, dir);
    start_norm[W] = 0.0;

    divisor = lib_normalize_vector(start_norm);
    if (ABSOLUTE(divisor) < EPSILON) {
	SET_COORD3(dir, 1.0, 0.0, 0.0);
	CROSS(start_norm, axis, dir);
	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(vert[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(vert[1], apex_pt, start_radius[1]);

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

    delta_angle = 2.0 * PI / (double)(2*u_resolution);
    for (i=1,angle=delta_angle;i<=2*u_resolution;++i,angle+=delta_angle) {
	lib_create_axis_rotate_matrix(mx, norm_axis, angle);
	lib_transform_vector(vert[2], start_radius[1], mx);
	ADD2_COORD3(vert[2], apex_pt);
	lib_transform_vector(norm[2], start_norm, mx);
	lib_output_polypatch(3, vert, norm);
	COPY_COORD3(vert[1], vert[2]);
	COPY_COORD3(norm[1], norm[2]);
	lib_transform_vector(vert[2], start_radius[0], mx);
	ADD2_COORD3(vert[2], base_pt);
	lib_output_polypatch(3, vert, norm);

	COPY_COORD3(vert[0], vert[2]);
	COPY_COORD3(norm[0], norm[2]);

	PLATFORM_MULTITASK();
    }
}

/*-----------------------------------------------------------------*/
/*
 * 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_POLYGONS, the surface is polygonalized and output.
 * (4*OUTPUT_RESOLUTION) polygons are output as rectangles by
 * lib_output_polypatch.
 */
void
lib_output_cylcone(base_pt, apex_pt, curve_format)
    COORD4 base_pt, apex_pt;
    int curve_format;
{
    object_ptr new_object;
    COORD4  axis;
    double  len, cottheta, xang, yang;

    if (format == OUTPUT_DELAYED) {
	/* Save all the pertinent information */
	new_object = (object_ptr)malloc(sizeof(struct object_struct));
	if (new_object == NULL)
	   /* Quietly fail */
	   return;
	new_object->object_type  = CONE_OBJ;
	new_object->curve_format = curve_format;
	new_object->surf_index   = texture_count;
	COPY_COORD4(new_object->object_data.cone.apex_pt, apex_pt);
	COPY_COORD4(new_object->object_data.cone.base_pt, base_pt);
	new_object->next_object = lib_objects;
	lib_objects = new_object;

    } else if (curve_format == OUTPUT_CURVES) {
	switch (format) {
	    case OUTPUT_VIDEO:
	    case OUTPUT_PLG:
		lib_output_polygon_cylcone(base_pt, apex_pt);
		break;

	    case OUTPUT_NFF:
		fprintf(outfile, "c\n" ) ;
		fprintf(outfile, "%g %g %g %g\n",
			base_pt[X], base_pt[Y], base_pt[Z], base_pt[W]);
		fprintf(outfile, "%g %g %g %g\n",
			apex_pt[X], apex_pt[Y], apex_pt[Z], apex_pt[W]);
		break;

	    case OUTPUT_POVRAY_10:
		/*
		Since POV-Ray uses infinite primitives, we will start
		with a cone aligned with the z-axis (QCone_Z) and figure
		out how to clip and scale it to match what we want
		*/
		if (apex_pt[W] < base_pt[W]) {
		    /* Put the bigger end at the top */
		    COPY_COORD4(axis, base_pt);
		    COPY_COORD4(base_pt, apex_pt);
		    COPY_COORD4(apex_pt, axis);
		    }
		/* Find the axis and axis length */
		SUB3_COORD3(axis, apex_pt, base_pt);
		len = lib_normalize_vector(axis);
		if (len < EPSILON) {
		   /* Degenerate cone/cylinder */
		   fprintf(outfile, "// degenerate cone/cylinder!  Ignored...\n");
		   break;
		   }
		if (ABSOLUTE(apex_pt[W] - base_pt[W]) < EPSILON) {
		   /* Treat this thing as a cylinder */
		   cottheta = len;
		   tab_indent();
		   fprintf(outfile, "object {\n");
		   tab_inc();

		   tab_indent();
		   fprintf(outfile, "quadric { <1 1 0> <0 0 0> <0 0 0> -1 } // cylinder\n");

		   tab_indent();
		   fprintf(outfile, "clipped_by {\n");
		   tab_inc();

		   tab_indent();
		   fprintf(outfile, "intersection {\n");
		   tab_inc();

		   tab_indent();
		   fprintf(outfile, "plane { <0 0 -1> 0 }\n");
		   tab_indent();
		   fprintf(outfile, "plane { <0 0  1> 1 }\n");

		   tab_dec();
		   tab_indent();
		   fprintf(outfile, "} // intersection\n");

		   tab_dec();
		   tab_indent();
		   fprintf(outfile, "} // clip\n");

		   tab_indent();
		   fprintf(outfile, "scale <%g %g 1>\n", base_pt[W], base_pt[W]);
		   }
		else {
		   /* Determine alignment */
		   cottheta = len / (apex_pt[W] - base_pt[W]);
		   tab_indent();
		   fprintf(outfile, "object {\n");
		   tab_inc();

		   tab_indent();
		   fprintf(outfile, "quadric{ <1 1 -1> <0 0 0> <0 0 0> 0 } // cone\n");

		   tab_indent();
		   fprintf(outfile, "clipped_by {\n");
		   tab_inc();

		   tab_indent();
		   fprintf(outfile, "intersection {\n");
		   tab_inc();

		   tab_indent();
		   fprintf(outfile, "plane { <0 0 -1> %g}\n", -base_pt[W]);
		   tab_indent();
		   fprintf(outfile, "plane { <0 0  1> %g}\n", apex_pt[W]);

		   tab_dec();
		   tab_indent();
		   fprintf(outfile, "} // intersection\n");

		   tab_dec();
		   tab_indent();
		   fprintf(outfile, "} // clip\n");

		   tab_indent();
		   fprintf(outfile, "translate <0 0 %g>\n", -base_pt[W]);
		   }

		tab_indent();
		fprintf(outfile, "scale <1 1 %g>\n", cottheta);

		len = sqrt(axis[X] * axis[X] + axis[Z] * axis[Z]);
		xang = -180.0 * asin(axis[Y]) / PI;
		if (len < EPSILON)
		  yang = 0.0;
		else
		  yang = 180.0 * acos(axis[Z] / len) / PI;
		if (axis[X] < 0)
		   yang = -yang;
		tab_indent();
		fprintf(outfile, "rotate <%g %g 0>\n", xang, yang);
		tab_indent();
		fprintf(outfile, "translate <%g %g %g>\n",
			base_pt[X], base_pt[Y], base_pt[Z]);
		if (texture_name != NULL) {
		   tab_indent();
		   fprintf(outfile, "texture { %s }\n", texture_name);
		   }

		tab_dec();
		tab_indent();
		fprintf(outfile, "} // object\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_POVRAY_20:
		/* of course if apex_pt[W] ~= base_pt[W], could do cylinder */
		tab_indent();
		fprintf(outfile, "cone {\n");
		tab_inc();

		tab_indent();
		fprintf(outfile, "<%g, %g, %g>, %g,\n",
			apex_pt[X], apex_pt[Y], apex_pt[Z], apex_pt[W]);
		tab_indent();
		fprintf(outfile, "<%g, %g, %g>, %g open\n",
			base_pt[X], base_pt[Y], base_pt[Z], base_pt[W]);
		if (texture_name != NULL) {
		   tab_indent();
		   fprintf(outfile, "texture { %s }\n", texture_name);
		   }

		tab_dec();
		tab_indent();
		fprintf(outfile, "}\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_POLYRAY:
		tab_indent();
		fprintf(outfile, "object {");
		tab_inc();

		tab_indent();
		if (base_pt[W] == apex_pt[W])
		   fprintf(outfile, "cylinder <%g, %g, %g>, <%g, %g, %g>, %g",
			  base_pt[X], base_pt[Y], base_pt[Z],
			  apex_pt[X], apex_pt[Y], apex_pt[Z], apex_pt[W]);
		else
		   fprintf(outfile, "cone <%g, %g, %g>, %g, <%g, %g, %g>, %g",
			  base_pt[X], base_pt[Y], base_pt[Z], base_pt[W],
			  apex_pt[X], apex_pt[Y], apex_pt[Z], apex_pt[W]);
		if (texture_name != NULL)
		   fprintf(outfile, " %s", texture_name);
		fprintf(outfile, "\n");

		tab_dec();
		tab_indent();
		fprintf(outfile, "}\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_VIVID:
		tab_indent();
		fprintf(outfile, "cone {");
		tab_inc();

		tab_indent();
		fprintf(outfile, " base %g %g %g base_radius %g\n",
			base_pt[X], base_pt[Y], base_pt[Z], base_pt[W]);
		tab_indent();
		fprintf(outfile, " apex %g %g %g apex_radius %g\n",
			apex_pt[X], apex_pt[Y], apex_pt[Z], apex_pt[W]);

		tab_dec();
		tab_indent();
		fprintf(outfile, "}\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_QRT:
		fprintf(outfile, "BEGIN_BBOX\n");
		lib_output_polygon_cylcone(base_pt, apex_pt);
		fprintf(outfile, "END_BBOX\n");
		break;

	    case OUTPUT_RAYSHADE:
		fprintf(outfile, "cone ");
		if (texture_name != NULL)
		    fprintf(outfile, "%s ", texture_name);
		fprintf(outfile, " %g %g %g %g %g %g %g %g\n",
			base_pt[W], base_pt[X], base_pt[Y], base_pt[Z],
			apex_pt[W], apex_pt[X], apex_pt[Y], apex_pt[Z]);
		break;

	    case OUTPUT_RTRACE:
		fprintf(outfile, "4 %d %g %g %g %g %g %g %g %g %g\n",
			texture_count, texture_ior,
			base_pt[X], base_pt[Y], base_pt[Z], base_pt[W],
			apex_pt[X], apex_pt[Y], apex_pt[Z], apex_pt[W]);
		break;

	    case OUTPUT_ART:
		if (base_pt[W] != apex_pt[W]) {
		    tab_indent();
		    fprintf(outfile, "cone {\n");
		    tab_inc();
		    tab_indent();
		    fprintf(outfile, "radius %f  center(%f, %f, %f)\n",
			    base_pt[W], base_pt[X], base_pt[Y], base_pt[Z]);
		    tab_indent();
		    fprintf(outfile, "radius %f  center(%f, %f, %f)\n",
			    apex_pt[W], apex_pt[X], apex_pt[Y], apex_pt[Z]);
		} else {
		    tab_indent();
		    fprintf(outfile, "cylinder {\n");
		    tab_inc();
		    tab_indent();
		    fprintf(outfile, "radius %f  center(%f, %f, %f)\n",
			    base_pt[W], base_pt[X], base_pt[Y], base_pt[Z]);
		    tab_indent();
		    fprintf(outfile, "center(%f, %f, %f)\n",
			    apex_pt[X], apex_pt[Y], apex_pt[Z]);
		}

		tab_dec();
		tab_indent();
		fprintf(outfile, "}\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_RAWTRI:
		lib_output_polygon_cylcone(base_pt, apex_pt);
		break;
	}
    } else {
	lib_output_polygon_cylcone(base_pt, apex_pt);
    }
}

/*-----------------------------------------------------------------*/
static void
disc_evaluator(trans, theta, v, r, vert)
    MATRIX trans;
    double theta, v, r;
    COORD3 vert;
{
    COORD3 tvert;

    /* Compute the position of the point */
    SET_COORD3(tvert, (r + v) * cos(theta), (r + v) * sin(theta), 0.0);
    lib_transform_vector(vert, tvert, trans);
}

/*-----------------------------------------------------------------*/
static void
lib_output_polygon_disc(center, normal, iradius, oradius)
    COORD3 center, normal;
    double iradius, oradius;
{
    double  divisor, u, v, delta_u, delta_v;
    MATRIX mx, imx;
    int i, j;
    COORD3 norm, vert[4];

    COPY_COORD3(norm, normal);
    if ((divisor = lib_normalize_vector(norm)) < EPSILON) {
	fprintf(stderr, "Bad disc normal\n");
	exit(1);
    }
    lib_create_canonical_matrix(mx, imx, center, norm);
    delta_u = 2.0 * PI / (double)(4 * u_resolution);
    delta_v = (oradius - iradius) / (double)(v_resolution);

    /* Dump out polygons */
    for (i=0,u=0.0;i<4*u_resolution;i++,u+=delta_u) {
	PLATFORM_MULTITASK();
	for (j=0,v=0.0;j<v_resolution;j++,v+=delta_v) {
	    disc_evaluator(imx, u, v, iradius, vert[0]);
	    disc_evaluator(imx, u+delta_u, v, iradius, vert[1]);
	    disc_evaluator(imx, u+delta_u, v+delta_v, iradius, vert[2]);
	    lib_output_polygon(3, vert);
	    COPY_COORD3(vert[1], vert[2]);
	    disc_evaluator(imx, u, v+delta_v, iradius, vert[2]);
	    lib_output_polygon(3, vert);
        }
    }
}

/*-----------------------------------------------------------------*/
void
lib_output_disc(center, normal, iradius, oradius, curve_format)
    COORD3 center, normal;
    double iradius, oradius;
    int curve_format;
{
    object_ptr new_object;
    COORD4  axis, base, apex;
    double  len, xang, yang;

    if (format == OUTPUT_DELAYED) {
	/* Save all the pertinent information */
	new_object = (object_ptr)malloc(sizeof(struct object_struct));
	if (new_object == NULL)
	    /* Quietly fail */
	    return;
	new_object->object_type  = DISC_OBJ;
	new_object->curve_format = curve_format;
	new_object->surf_index   = texture_count;
	COPY_COORD4(new_object->object_data.disc.center, center);
	COPY_COORD4(new_object->object_data.disc.normal, normal);
	new_object->object_data.disc.iradius = iradius;
	new_object->object_data.disc.iradius = oradius;
	new_object->next_object = lib_objects;
	lib_objects = new_object;
    } else if (curve_format == OUTPUT_CURVES) {
	switch (format) {
	    case OUTPUT_VIDEO:
	    case OUTPUT_NFF:
	    case OUTPUT_PLG:
	    case OUTPUT_VIVID:
	    case OUTPUT_RAYSHADE:
	    case OUTPUT_RAWTRI:
		lib_output_polygon_disc(center, normal, iradius, oradius);
		break;

	    case OUTPUT_POVRAY_10:
		/* A disc is a plane intersected with either one or two
		 * spheres
		 */
		COPY_COORD3(axis, normal);
		len = lib_normalize_vector(axis);
		tab_indent();
		fprintf(outfile, "object {\n");
		tab_inc();

		tab_indent();
		fprintf(outfile, "plane { <0 0 1> 1 }\n");

		tab_indent();
		fprintf(outfile, "clipped_by {\n");
		tab_inc();

		if (iradius > 0.0) {
		    tab_indent();
		    fprintf(outfile, "intersection {\n");
		    tab_inc();

		    tab_indent();
		    fprintf(outfile, "sphere { <0 0 0> %g inverse }\n",
			    iradius);
		    tab_indent();
		    fprintf(outfile, "sphere { <0 0 1> %g }\n", oradius);

		    tab_dec();
		    tab_indent();
		    fprintf(outfile, "} // intersection\n");
		}
		else {
		    tab_indent();
		    fprintf(outfile, "object { sphere { <0 0 0> %g } }\n",
			    oradius);
		}

		tab_dec();
		tab_indent();
		fprintf(outfile, "} // clip\n");

		len = sqrt(axis[X] * axis[X] + axis[Z] * axis[Z]);
		xang = -180.0 * asin(axis[Y]) / PI;
		yang = 180.0 * acos(axis[Z] / len) / PI;
		if (axis[X] < 0)
		    yang = -yang;
		tab_indent();
		fprintf(outfile, "rotate <%g %g 0>\n", xang, yang);
		tab_indent();
		fprintf(outfile, "translate <%g %g %g>\n",
			center[X], center[Y], center[Z]);

		if (texture_name != NULL) {
		    tab_indent();
		    fprintf(outfile, "texture { %s }", texture_name);
		}

		tab_dec();
		tab_indent();
		fprintf(outfile, "} // object - disc\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_POVRAY_20:
		/* disc <center> <normalVector> radius [holeRadius] */
		tab_indent();
		fprintf(outfile, "disc { <%g, %g, %g>",
			center[X], center[Y], center[Z]);
		fprintf(outfile, " <%g, %g, %g>",
			normal[X], normal[Y], normal[Z]);
		fprintf(outfile, " %g", oradius);
		if (iradius > 0.0)
		    fprintf(outfile, ", %g", iradius);
		if (texture_name != NULL)
		    fprintf(outfile, " texture { %s }", texture_name);
		fprintf(outfile, " }\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_POLYRAY:
		tab_indent();
		fprintf(outfile, "object { disc <%g, %g, %g>,",
			center[X], center[Y], center[Z]);
		fprintf(outfile, " <%g, %g, %g>,",
			normal[X], normal[Y], normal[Z]);
		if (iradius > 0.0)
		    fprintf(outfile, " %g,", iradius);
		fprintf(outfile, " %g", oradius);
		if (texture_name != NULL)
		    fprintf(outfile, " %s", texture_name);
		fprintf(outfile, " }\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_QRT:
		fprintf(outfile, "BEGIN_BBOX\n");
		lib_output_polygon_disc(center, normal, iradius, oradius);
		fprintf(outfile, "END_BBOX\n");
		break;

	    case OUTPUT_RTRACE:
		COPY_COORD3(base, center);
		base[W] = iradius;
		apex[X] = center[X] + normal[X] * EPSILON;
		apex[Y] = center[Y] + normal[Y] * EPSILON;
		apex[Z] = center[Z] + normal[Z] * EPSILON;
		apex[W] = oradius;
		lib_output_cylcone(base, apex, curve_format);
		break;
	}
    } else {
	lib_output_polygon_disc(center, normal, iradius, oradius);
    }
}


#if defined (applec)
/* Mac MPW C "prefers" to keep code segments <32k (small code model) */
#pragma segment "lib2"
#endif


/*-----------------------------------------------------------------*/
static void
sq_sphere_val(a1, a2, a3, n, e, u, v, P)
    double a1, a2, a3, n, e, u, v;
    COORD3 P;
{
    double cu, su, cv, sv;
    double icu, isu, icv, isv;

    cu = cos(u); su = sin(u);
    cv = cos(v); sv = sin(v);
    icu = SGN(cu); isu = SGN(su);
    icv = SGN(cv); isv = SGN(sv);
    cu = fabs(cu); cv = fabs(cv);
    su = fabs(su); sv = fabs(sv);
    P[X] = a1 * POW(cv, n) * POW(cu, e) * icv * icu;
    P[Y] = a2 * POW(cv, n) * POW(su, e) * icv * isu;
    P[Z] = a3 * POW(sv, n) * isv;
}

/*-----------------------------------------------------------------*/
static void
sq_sphere_norm(a1, a2, a3, n, e, u, v, N)
    double a1, a2, a3, n, e, u, v;
    COORD3 N;
{
    double cu, su, cv, sv;
    double icu, isu, icv, isv;

    cu = cos(u); su = sin(u);
    cv = cos(v); sv = sin(v);
    icu = SGN(cu); isu = SGN(su);
    icv = SGN(cv); isv = SGN(sv);

    /* May be some singularities in the values, lets catch them & put
      a fudged normal into N */
    if (e < 2 || n < 2) {
	if (ABSOLUTE(cu) < 1.0e-3 || ABSOLUTE(su) < 1.0e-3 ||
	    ABSOLUTE(cu) < 1.0e-3 || ABSOLUTE(su) < 1.0e-3) {
	   SET_COORD3(N, cu*cv, su*cv, sv);
	   lib_normalize_vector(N);
	   return;
	}
    }

    cu = fabs(cu); cv = fabs(cv);
    su = fabs(su); sv = fabs(sv);

    N[X] = a1 * POW(cv, 2-n) * POW(cu, 2-e) * icv * icu;
    N[Y] = a2 * POW(cv, 2-n) * POW(su, 2-e) * icv * isu;
    N[Z] = a3 * POW(sv, 2-n) * isv;
    lib_normalize_vector(N);
}

/*-----------------------------------------------------------------*/
void
lib_output_sq_sphere(center_pt, a1, a2, a3, n, e)
    COORD3 center_pt;
    double a1, a2, a3, n, e;
{
    object_ptr new_object;
    int i, j, u_res, v_res;
    double u, delta_u, v, delta_v;
    COORD3 verts[4], norms[4];

    if (format == OUTPUT_DELAYED) {
       /* Save all the pertinent information */
       new_object = (object_ptr)malloc(sizeof(struct object_struct));
       if (new_object == NULL)
	   /* Quietly fail */
	   return;
       new_object->object_type  = SUPERQ_OBJ;
       new_object->curve_format = OUTPUT_PATCHES;
       new_object->surf_index   = texture_count;
       COPY_COORD4(new_object->object_data.superq.center_pt, center_pt);
       new_object->object_data.superq.a1 = a1;
       new_object->object_data.superq.a2 = a2;
       new_object->object_data.superq.a3 = a3;
       new_object->object_data.superq.n  = n;
       new_object->object_data.superq.e  = e;
       new_object->next_object = lib_objects;
       lib_objects = new_object;
       return;
    }

    u_res = 4 * u_resolution;
    v_res = 4 * v_resolution;
    delta_u = 2.0 * PI / (double)u_res;
    delta_v = PI / (double)v_res;

    for (i=0,u=0.0;i<u_res;i++,u+=delta_u) {
	PLATFORM_MULTITASK();
	for (j=0,v=-PI/2.0;j<v_res;j++,v+=delta_v) {
	    if (j == 0) {
		sq_sphere_val(a1, a2, a3, n, e, u, v, verts[0]);
		sq_sphere_norm(a1, a2, a3, n, e, u, v, norms[0]);
		sq_sphere_val(a1, a2, a3, n, e, u, v+delta_v, verts[1]);
		sq_sphere_norm(a1, a2, a3, n, e, u, v+delta_v, norms[1]);
		sq_sphere_val(a1, a2, a3, n, e, u+delta_u, v+delta_v, verts[2]);
		sq_sphere_norm(a1, a2, a3, n, e, u+delta_u, v+delta_v,norms[2]);
		ADD3_COORD3(verts[0], verts[0], center_pt);
		ADD3_COORD3(verts[1], verts[1], center_pt);
		ADD3_COORD3(verts[2], verts[2], center_pt);
		lib_output_polypatch(3, verts, norms);
	    } else if (j == v_res-1) {
		sq_sphere_val(a1, a2, a3, n, e, u, v, verts[0]);
		sq_sphere_norm(a1, a2, a3, n, e, u, v, norms[0]);
		sq_sphere_val(a1, a2, a3, n, e, u, v+delta_v, verts[1]);
		sq_sphere_norm(a1, a2, a3, n, e, u, v+delta_v, norms[1]);
		sq_sphere_val(a1, a2, a3, n, e, u+delta_u, v, verts[2]);
		sq_sphere_norm(a1, a2, a3, n, e, u+delta_u, v, norms[2]);
		ADD3_COORD3(verts[0], verts[0], center_pt);
		ADD3_COORD3(verts[1], verts[1], center_pt);
		ADD3_COORD3(verts[2], verts[2], center_pt);
		lib_output_polypatch(3, verts, norms);
	    } else {
		sq_sphere_val(a1, a2, a3, n, e, u, v, verts[0]);
		sq_sphere_norm(a1, a2, a3, n, e, u, v, norms[0]);
		sq_sphere_val(a1, a2, a3, n, e, u, v+delta_v, verts[1]);
		sq_sphere_norm(a1, a2, a3, n, e, u, v+delta_v, norms[1]);
		sq_sphere_val(a1, a2, a3, n, e, u+delta_u, v+delta_v, verts[2]);
		sq_sphere_norm(a1, a2, a3, n, e, u+delta_u, v+delta_v,norms[2]);
		ADD3_COORD3(verts[0], verts[0], center_pt);
		ADD3_COORD3(verts[1], verts[1], center_pt);
		ADD3_COORD3(verts[2], verts[2], center_pt);
		lib_output_polypatch(3, verts, norms);
		COPY_COORD3(verts[1], verts[2]);
		COPY_COORD3(norms[1], norms[2]);
		sq_sphere_val(a1, a2, a3, n, e, u+delta_u, v, verts[2]);
		sq_sphere_norm(a1, a2, a3, n, e, u+delta_u, v, norms[2]);
		ADD3_COORD3(verts[2], verts[2], center_pt);
		lib_output_polypatch(3, verts, norms);
	    }
	}
    }
}

/*-----------------------------------------------------------------*/
static void
lib_output_polygon_sphere(center_pt)
    COORD4 center_pt;
{
    double  angle;
    COORD3  edge_norm[3], edge_pt[3];
    long    num_face, num_edge, num_tri, num_vert;
    COORD3  *x_axis, *y_axis, **pt;
    COORD3  mid_axis;
    MATRIX  rot_mx;
    long    u_pol, v_pol;

    /* Allocate storage for the polygon vertices */
    x_axis = (COORD3 *)malloc((u_resolution+1) * sizeof(COORD3));
    y_axis = (COORD3 *)malloc((v_resolution+1) * sizeof(COORD3));
    pt     = (COORD3 **)malloc((u_resolution+1) * sizeof(COORD3 *));
    if (x_axis == NULL || y_axis == NULL || pt == NULL) {
	fprintf(stderr, "Failed to allocate polygon data\n");
	exit(1);
    }

    for (num_edge=0;num_edge<u_resolution+1;num_edge++) {
	pt[num_edge] = (COORD3 *)malloc((v_resolution+1) * sizeof(COORD3));
	if (pt[num_edge] == NULL) {
	    fprintf(stderr, "Failed to allocate polygon data\n");
	    exit(1);
        }
    }

    /* calculate axes used to find grid points */
    for (num_edge=0;num_edge<=u_resolution;++num_edge) {
	angle = (PI/4.0) * (2.0*(double)num_edge/u_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);
    }

    for (num_edge=0;num_edge<=v_resolution;++num_edge) {
	angle = (PI/4.0) * (2.0*(double)num_edge/v_resolution - 1.0);
	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<=u_resolution;++u_pol) {
	for (v_pol=0;v_pol<=u_resolution;++v_pol) {
	    CROSS(pt[u_pol][v_pol], x_axis[u_pol], y_axis[v_pol]);
	    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<=u_resolution;++u_pol) {
	    for (v_pol=0;v_pol<=v_resolution;++v_pol) {
	        lib_rotate_cube_face(pt[u_pol][v_pol], Z_AXIS, num_face);
	    }
	}

	/* output grid */
	for (u_pol=0;u_pol<u_resolution;++u_pol) {
	    for (v_pol=0;v_pol<v_resolution;++v_pol) {
	        PLATFORM_MULTITASK();
	        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);
		}
	    }
	}
    }

    /* Release any memory used */
    for (num_edge=0;num_edge<u_resolution+1;num_edge++)
	free(pt[num_edge]);
    free(pt);
    free(y_axis);
    free(x_axis);
}

/*-----------------------------------------------------------------*/
/*
 * 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_POLYGONS, 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*u_resolution*v_resolution) polygons are output as triangles
 * using lib_output_polypatch.
 */
void
lib_output_sphere(center_pt, curve_format)
    COORD4 center_pt;
    int curve_format;
{
    object_ptr new_object;

    if (format == OUTPUT_DELAYED) {
	/* Save all the pertinent information */
	new_object = (object_ptr)malloc(sizeof(struct object_struct));
	if (new_object == NULL)
	    /* Quietly fail */
	    return;
	new_object->object_type  = SPHERE_OBJ;
	new_object->curve_format = curve_format;
	new_object->surf_index   = texture_count;
	COPY_COORD4(new_object->object_data.sphere.center_pt, center_pt);
	new_object->next_object = lib_objects;
	lib_objects = new_object;
    }
    else if (curve_format == OUTPUT_CURVES) {
	switch (format) {
	    case OUTPUT_VIDEO:
	    case OUTPUT_PLG:
		lib_output_polygon_sphere(center_pt);
		break;

	    case OUTPUT_NFF:
		fprintf(outfile, "s %g %g %g %g\n",
			center_pt[X], center_pt[Y], center_pt[Z], center_pt[W]);
		break;

	    case OUTPUT_POVRAY_10:
		tab_indent();
		fprintf(outfile, "object { sphere { <%g %g %g> %g }",
			center_pt[X], center_pt[Y], center_pt[Z], center_pt[W]);
		if (texture_name != NULL)
		    fprintf(outfile, " texture { %s }", texture_name);
		fprintf(outfile, " }\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_POVRAY_20:
		tab_indent();
		fprintf(outfile, "sphere { <%g, %g, %g>, %g",
			center_pt[X], center_pt[Y], center_pt[Z], center_pt[W]);
		if (texture_name != NULL)
		    fprintf(outfile, " texture { %s }", texture_name);
		fprintf(outfile, " }\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_POLYRAY:
		tab_indent();
		fprintf(outfile, "object { sphere <%g, %g, %g>, %g",
			center_pt[X], center_pt[Y], center_pt[Z], center_pt[W]);
		if (texture_name != NULL)
		    fprintf(outfile, " %s", texture_name);
		fprintf(outfile, " }\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_VIVID:
		tab_indent();
		fprintf(outfile, "sphere { center %g %g %g radius %g }\n",
			center_pt[X], center_pt[Y], center_pt[Z], center_pt[W]);
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_QRT:
		tab_indent();
		fprintf(outfile, "sphere ( loc = (%g, %g, %g), radius = %g )\n",
			center_pt[X], center_pt[Y], center_pt[Z], center_pt[W]);
		break;

	    case OUTPUT_RAYSHADE:
		fprintf(outfile, "sphere ");
		if (texture_name != NULL)
		    fprintf(outfile, "%s ", texture_name);
		fprintf(outfile, " %g %g %g %g\n",
			center_pt[W], center_pt[X], center_pt[Y], center_pt[Z]);
		break;

	    case OUTPUT_RTRACE:
		fprintf(outfile, "1 %d %g %g %g %g %g\n",
			texture_count, texture_ior,
			center_pt[X], center_pt[Y], center_pt[Z], center_pt[W]);
		break;

	    case OUTPUT_ART:
		tab_indent();
		fprintf(outfile, "sphere {\n");
		tab_inc();

		tab_indent();
		fprintf(outfile, "radius %f\n", center_pt[W]);
		tab_indent();
		fprintf(outfile, "center(%f, %f, %f)\n",
			center_pt[X], center_pt[Y], center_pt[Z]);

		tab_dec();
		tab_indent();
		fprintf(outfile, "}\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_RAWTRI:
		lib_output_polygon_sphere(center_pt);
		break;
	}
    } else {
	lib_output_polygon_sphere(center_pt);
    }
}


static unsigned int hfcount = 0;

/*-----------------------------------------------------------------*/
static char *
create_height_file(filename, height, width, data, type)
    char *filename;
    unsigned int height, width;
    float **data;
    int type;
{
    FILE *file;
    float v;
    unsigned int i, j;
    unsigned char r, g, b;
    unsigned char tgaheader[18];

    if (filename == NULL) {
	/* Need to create a new name for the height file */
	filename = malloc(10 * sizeof(char));
	if (filename == NULL) return NULL;
	sprintf(filename, "hf%03d.tga", hfcount++);
    }
    if ((file = fopen(filename, "wb")) == NULL) return NULL;
    if (type == 0) {
	/* Targa style height field for POV-Ray or Polyray */
	memset(tgaheader, 0, 18);
	tgaheader[2] = 2;
	tgaheader[12] = (unsigned char)(width & 0xFF);
	tgaheader[13] = (unsigned char)((width >> 8) & 0xFF);
	tgaheader[14] = (unsigned char)(height & 0xFF);
	tgaheader[15] = (unsigned char)((height >> 8) & 0xFF);
	tgaheader[16] = 24;
	tgaheader[17] = 0x20;
	fwrite(tgaheader, 18, 1, file);
	for (i=0;i<height;i++) {
	    PLATFORM_MULTITASK();
	    for (j=0;j<width;j++) {
		v = data[i][j];
		if (v < -128.0) v = -128.0;
		if (v > 127.0) v = 127.0;
		v += 128.0;
		r = v;
		v -= (float)r;
		g = (unsigned char)(256.0 * v);
		b = 0;
		fputc(b, file);
		fputc(g, file);
		fputc(r, file);
	    }
	}
    } else {
	/* Only square height fields in RayShade */
	if (height < width) width = height;
	else if (width < height) height = width;

	/* Start by storing the size as an int */
	fwrite(&height, sizeof(int), 1, file);

	/* Now store height values as native floats */
	for (i=0;i<height;i++)
	    for (j=0;j<width;j++)
		fwrite(&data[i][j], sizeof(float), 1, file);
    }
    fclose(file);

    return filename;
}

/*-----------------------------------------------------------------*/
static void
find_height_min_max(height, width, data, y0, y1)
    unsigned int height, width;
    float **data;
    float *y0, *y1;
{
    unsigned int i, j;
    float v;

    *y0 = *y1 = data[0][0];
    for (i=0;i<height;i++) {
        for (j=0;j<width;j++) {
	    v = data[i][j];
	    if (v < *y0) *y0 = v;
	    if (v > *y1) *y1 = v;
        }
    }
}

/*-----------------------------------------------------------------*/
static void
lib_output_polygon_height(height, width, data, x0, x1, y0, y1, z0, z1)
    unsigned int height, width;
    float **data;
    double x0, x1, y0, y1, z0, z1;
{
    unsigned int i, j;
    double xdelta, zdelta;
    COORD3 verts[3];

#if defined (applec)
#pragma unused (y1)
#endif /* applec */

    xdelta = (x1 - x0) / (double)(width - 1);
    zdelta = (z1 - z0) / (double)(height - 1);
    for (i=0;i<height-1;i++) {
	for (j=0;j<width-1;j++) {
	    PLATFORM_MULTITASK();
	    SET_COORD3(verts[0], x0 + j * xdelta, y0 + data[i][j],
		       z0 + i * zdelta);
	    SET_COORD3(verts[1], x0 + (j+1) * xdelta, y0 + data[i][j+1],
		       z0 + i * zdelta);
	    SET_COORD3(verts[2], x0 + (j+1) * xdelta, y0 + data[i+1][j+1],
		       z0 + (i + 1) * zdelta);
	    lib_output_polygon(3, verts);
	    COPY_COORD3(verts[1], verts[2]);
	    SET_COORD3(verts[2], x0 + j * xdelta, y0 + data[i+1][j],
		       z0 + (i + 1) * zdelta);
	    lib_output_polygon(3, verts);
	}
    }
}

/*-----------------------------------------------------------------*/
void
lib_output_height(filename, data, height, width, x0, x1, y0, y1, z0, z1)
    char *filename;
    float **data;
    unsigned int height, width;
    double x0, x1;
    double y0, y1;
    double z0, z1;
{
    object_ptr new_object;

    if (format == OUTPUT_DELAYED) {
	filename = create_height_file(filename, height, width, data, 0);
	if (filename == NULL) return;

	/* Save all the pertinent information */
	new_object = (object_ptr)malloc(sizeof(struct object_struct));
	if (new_object == NULL)
	    /* Quietly fail */
	    return;
	new_object->object_type  = HEIGHT_OBJ;
	new_object->curve_format = OUTPUT_CURVES;
	new_object->surf_index   = texture_count;
	new_object->object_data.height.width = width;
	new_object->object_data.height.height = height;
	new_object->object_data.height.data = data;
	new_object->object_data.height.filename = filename;
	new_object->object_data.height.x0 = x0;
	new_object->object_data.height.x1 = x1;
	new_object->object_data.height.y0 = y0;
	new_object->object_data.height.y1 = y1;
	new_object->object_data.height.z0 = z0;
	new_object->object_data.height.z1 = z1;
	new_object->next_object = lib_objects;
	lib_objects = new_object;
    } else {
	switch (format) {
	    case OUTPUT_VIDEO:
	    case OUTPUT_NFF:
	    case OUTPUT_PLG:
	    case OUTPUT_QRT:
	    case OUTPUT_RTRACE:
	    case OUTPUT_VIVID:
	    case OUTPUT_RAWTRI:
		lib_output_polygon_height(height, width, data,
					  x0, x1, y0, y1, z0, z1);
		break;
	    case OUTPUT_POVRAY_10:
	    case OUTPUT_POVRAY_20:
		filename = create_height_file(filename, height, width, data, 0);
		if (filename == NULL) return;

		tab_indent();
		fprintf(outfile, "object {\n");
		tab_inc();

		tab_indent();
		fprintf(outfile, "height_field { tga \"%s\" }", filename);
		if (format == OUTPUT_POVRAY_10) {
		    tab_indent();
		    fprintf(outfile, "scale <%g %g %g>\n",
			    fabs(x1 - x0), fabs(y1 - y0), fabs(z1 - z0));
		    tab_indent();
		    fprintf(outfile, "translate <%g %g %g>\n", x0, y0, z0);
		} else {
		    tab_indent();
		    fprintf(outfile, "scale <%g, %g, %g>\n",
			    fabs(x1 - x0), fabs(y1 - y0), fabs(z1 - z0));
		    tab_indent();
		    fprintf(outfile, "translate <%g, %g, %g>\n", x0, y0, z0);
		}

		if (texture_name != NULL) {
		    tab_indent();
		    fprintf(outfile, "texture { %s }", texture_name);
		}

		tab_dec();
		tab_indent();
		fprintf(outfile, "} // object - Height Field\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_POLYRAY:
		filename = create_height_file(filename, height, width, data, 0);
		if (filename == NULL) return;
		tab_indent();
		fprintf(outfile, "object { height_field \"%s\" ", filename);
		fprintf(outfile, "scale <%g, 1, %g> ",
			fabs(x1-x0), fabs(z1-z0));
		fprintf(outfile, "translate <%g, %g, %g> ", x0, y0, z0);
		if (texture_name != NULL)
		    fprintf(outfile, " %s", texture_name);
		fprintf(outfile, " }\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_RAYSHADE:
		filename = create_height_file(filename, height, width, data, 1);
		if (filename == NULL) return;
		fprintf(outfile, "heightfield ");
		if (texture_name != NULL)
		    fprintf(outfile, " %s", texture_name);
		fprintf(outfile, "\"%s\" ", filename);
		fprintf(outfile, "rotate 1 0 0 90 ");
		fprintf(outfile, "scale  %g 1 %g ",
			fabs(x1 - x0), fabs(z1 - z0));
		fprintf(outfile, "translate  %g %g %g ", x0, y0, z0);
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_ART:
		filename = create_height_file(filename, height, width, data, 1);
		if (filename == NULL) return;

		tab_indent();
		fprintf(outfile, "geometry {\n");
		tab_inc();

		tab_indent();
		fprintf(outfile, "translate(%g, %g, %g)\n", x0, y0, z0);
		tab_indent();
		fprintf(outfile, "scale(%g, 1, %g)\n",
			fabs(x1 - x0), fabs(z1 - z0));
		tab_indent();
		fprintf(outfile, "rotate(-90, x)\n");
		tab_indent();
		fprintf(outfile, "heightfield \"%s\"\n ", filename);

		tab_dec();
		tab_indent();
		fprintf(outfile, "}\n");
		fprintf(outfile, "\n");
		break;
	}
    }
}

/*-----------------------------------------------------------------*/
static void
torus_evaluator(trans, theta, phi, r0, r1, vert, norm)
    MATRIX trans;
    double theta, phi, r0, r1;
    COORD3 vert, norm;
{
    COORD3 v0, v1, tvert, tnorm;

    /* Compute the position of the point */
    SET_COORD3(tvert, (r0 + r1 * sin(theta)) * cos(phi),
		      (r0 + r1 * sin(theta)) * sin(phi),
		      r1 * cos(theta));
    /* Compute the normal at that point */
    SET_COORD3(v0, r1*cos(theta)*cos(phi),
		   r1*cos(theta)*sin(phi),
		  -r1*sin(theta));
    SET_COORD3(v1,-(r0+r1*sin(theta))*sin(phi),
		   (r0+r1*sin(theta))*cos(phi),
		   0.0);
    CROSS(tnorm, v0, v1);
    lib_normalize_vector(tnorm);
    lib_transform_point(vert, tvert, trans);
    lib_transform_vector(norm, tnorm, trans);
}

/*-----------------------------------------------------------------*/
static void
lib_output_polygon_torus(center, normal, iradius, oradius)
    COORD3 center, normal;
    double iradius, oradius;
{
    double  divisor, u, v, delta_u, delta_v;
    MATRIX mx, imx;
    int i, j;
    COORD3 vert[4], norm[4];

    if ((divisor = lib_normalize_vector(normal)) < EPSILON) {
	fprintf(stderr, "Bad torus normal\n");
	exit(1);
    }
    lib_create_canonical_matrix(mx, imx, center, normal);
    delta_u = 2.0 * PI / (double)u_resolution;
    delta_v = 2.0 * PI / (double)v_resolution;

    /* Dump out polygons */
    for (i=0,u=0.0;i<u_resolution;i++,u+=delta_u) {
	PLATFORM_MULTITASK();
	for (j=0,v=0.0;j<v_resolution;j++,v+=delta_v) {
	    torus_evaluator(imx, u, v, iradius, oradius, vert[0], norm[0]);
	    torus_evaluator(imx, u, v+delta_v, iradius, oradius,
			vert[1], norm[1]);
	    torus_evaluator(imx, u+delta_u, v+delta_v,
			iradius, oradius, vert[2], norm[2]);
	    lib_output_polypatch(3, vert, norm);
	    COPY_COORD3(vert[1], vert[2]);
	    COPY_COORD3(norm[1], norm[2]);
	    torus_evaluator(imx, u+delta_u, v, iradius, oradius,
			vert[2], norm[2]);
	    lib_output_polypatch(3, vert, norm);
	}
    }
}

/*-----------------------------------------------------------------*/
void
lib_output_torus(center, normal, iradius, oradius, curve_format)
    COORD3 center, normal;
    double iradius, oradius;
    int curve_format;
{
    object_ptr new_object;
    double len, xang, zang;

    if (format == OUTPUT_DELAYED) {
	/* Save all the pertinent information */
	new_object = (object_ptr)malloc(sizeof(struct object_struct));
	if (new_object == NULL)
	    /* Quietly fail */
	    return;
	new_object->object_type  = TORUS_OBJ;
	new_object->curve_format = curve_format;
	new_object->surf_index   = texture_count;
	COPY_COORD3(new_object->object_data.torus.center, center);
	COPY_COORD3(new_object->object_data.torus.normal, normal);
	new_object->object_data.torus.iradius = iradius;
	new_object->object_data.torus.oradius = oradius;
	new_object->next_object = lib_objects;
	lib_objects = new_object;
    } else if (curve_format == OUTPUT_CURVES) {
	switch (format) {
	    case OUTPUT_VIDEO:
	    case OUTPUT_NFF:
	    case OUTPUT_VIVID:
	    case OUTPUT_QRT:
	    case OUTPUT_POVRAY_10:
	    case OUTPUT_PLG:
	    case OUTPUT_RTRACE:
	    case OUTPUT_RAWTRI:
		lib_output_polygon_torus(center, normal, iradius, oradius);
		break;
	    case OUTPUT_POVRAY_20:
		/*
		 A torus object lies in the x-z plane.  We need to determine
		 the angles of rotation to get it lined up with "normal".
		 */
		tab_indent();
		fprintf(outfile, "torus {\n");
		tab_inc();

		tab_indent();
		fprintf(outfile, "%g, %g\n", iradius, oradius);

		(void)lib_normalize_vector(normal);
		len = sqrt(normal[X] * normal[X] + normal[Y] * normal[Y]);
		xang = 180.0 * asin(normal[Z]) / PI;
		if (len < EPSILON)
		    zang = 0.0;
		else
		    zang = -180.0 * acos(normal[Y] / len) / PI;
		if (normal[X] < 0)
		    zang = -zang;

		if (ABSOLUTE(xang) > EPSILON || ABSOLUTE(zang) > EPSILON) {
		    tab_indent();
		    fprintf(outfile, "rotate <%g, 0, %g>\n", xang, zang);
		}

		if (ABSOLUTE(center[X]) > EPSILON ||
		    ABSOLUTE(center[Y]) > EPSILON ||
		    ABSOLUTE(center[Z]) > EPSILON) {
		    tab_indent();
		    fprintf(outfile, "translate <%g, %g, %g>\n",
			    center[X], center[Y], center[Z]);
		}

		if (texture_name != NULL) {
		    tab_indent();
		    fprintf(outfile, "texture { %s }", texture_name);
		}
		fprintf(outfile, "\n");

		tab_dec();
		tab_indent();
		fprintf(outfile, "} // torus\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_POLYRAY:
		tab_indent();
		fprintf(outfile, "object { torus %g, %g", iradius, oradius);
		fprintf(outfile, ", <%g, %g, %g>, <%g, %g, %g>",
			center[X], center[Y], center[Z],
			normal[X], normal[Y], normal[Z]);
		if (texture_name != NULL)
		    fprintf(outfile, " %s", texture_name);
		fprintf(outfile, " }\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_RAYSHADE:
		fprintf(outfile, "torus ");
		if (texture_name != NULL)
		    fprintf(outfile, "%s ", texture_name);
		fprintf(outfile, " %g %g %g %g %g %g %g %g\n",
			iradius, oradius,
			center[X], center[Y], center[Z],
			normal[X], normal[Y], normal[Z]);
		break;

	    case OUTPUT_ART:
		tab_indent();
		fprintf(outfile, "torus {\n");
		tab_inc();

		tab_indent();
		fprintf(outfile,
		"radius %g  center(%g, %g, %g)  radius %g  center(%g, %g, %g)\n",
			  iradius, oradius,
			  center[X], center[Y], center[Z],
			  normal[X], normal[Y], normal[Z]);
		tab_dec();
		tab_indent();
		fprintf(outfile, "}\n");
		fprintf(outfile, "\n");
		break;
	}
    } else {
	lib_output_polygon_torus(center, normal, iradius, oradius);
    }
}


#if defined (applec)
/* Mac MPW C "prefers" to keep code segments <32k (small code model) */
#pragma segment "lib3"
#endif


/*-----------------------------------------------------------------*/
/* Generate a box as a set of 4-sided polygons */
static void
lib_output_polygon_box(p1, p2)
    COORD3 p1, p2;
{
    COORD3 box_verts[4];

    /* Sides */
    SET_COORD3(box_verts[0], p1[X], p1[Y], p1[Z]);
    SET_COORD3(box_verts[1], p1[X], p2[Y], p1[Z]);
    SET_COORD3(box_verts[2], p1[X], p2[Y], p2[Z]);
    SET_COORD3(box_verts[3], p1[X], p1[Y], p2[Z]);
    lib_output_polygon(4, box_verts);
    SET_COORD3(box_verts[3], p2[X], p1[Y], p1[Z]);
    SET_COORD3(box_verts[2], p2[X], p2[Y], p1[Z]);
    SET_COORD3(box_verts[1], p2[X], p2[Y], p2[Z]);
    SET_COORD3(box_verts[0], p2[X], p1[Y], p2[Z]);
    lib_output_polygon(4, box_verts);

    /* Front/Back */
    SET_COORD3(box_verts[0], p1[X], p1[Y], p1[Z]);
    SET_COORD3(box_verts[1], p1[X], p2[Y], p1[Z]);
    SET_COORD3(box_verts[2], p2[X], p2[Y], p1[Z]);
    SET_COORD3(box_verts[3], p2[X], p1[Y], p1[Z]);
    lib_output_polygon(4, box_verts);
    SET_COORD3(box_verts[3], p1[X], p1[Y], p2[Z]);
    SET_COORD3(box_verts[2], p1[X], p2[Y], p2[Z]);
    SET_COORD3(box_verts[1], p2[X], p2[Y], p2[Z]);
    SET_COORD3(box_verts[0], p2[X], p1[Y], p2[Z]);
    lib_output_polygon(4, box_verts);

    /* Top/Bottom */
    SET_COORD3(box_verts[0], p1[X], p1[Y], p1[Z]);
    SET_COORD3(box_verts[1], p1[X], p1[Y], p2[Z]);
    SET_COORD3(box_verts[2], p2[X], p1[Y], p2[Z]);
    SET_COORD3(box_verts[3], p2[X], p1[Y], p1[Z]);
    lib_output_polygon(4, box_verts);
    SET_COORD3(box_verts[3], p1[X], p2[Y], p1[Z]);
    SET_COORD3(box_verts[2], p1[X], p2[Y], p2[Z]);
    SET_COORD3(box_verts[1], p2[X], p2[Y], p2[Z]);
    SET_COORD3(box_verts[0], p2[X], p2[Y], p1[Z]);
    lib_output_polygon(4, box_verts);
}

/*-----------------------------------------------------------------*/
/* Output box.  A box is defined by a diagonally opposite corners. */
void
lib_output_box(p1, p2)
    COORD3 p1, p2;
{
    object_ptr new_object;

    if (format == OUTPUT_DELAYED) {
	/* Save all the pertinent information */
	new_object = (object_ptr)malloc(sizeof(struct object_struct));
	if (new_object == NULL)
	    /* Quietly fail */
	    return;
	new_object->object_type  = BOX_OBJ;
	new_object->curve_format = OUTPUT_PATCHES;
	new_object->surf_index   = texture_count;
	COPY_COORD3(new_object->object_data.box.point1, p1);
	COPY_COORD3(new_object->object_data.box.point2, p2);
	new_object->next_object = lib_objects;
	lib_objects = new_object;
    } else {
	switch (format) {
	    case OUTPUT_VIDEO:
	    case OUTPUT_NFF:
	    case OUTPUT_VIVID:
	    case OUTPUT_PLG:
	    case OUTPUT_RAWTRI:
		lib_output_polygon_box(p1, p2);
		break;

	    case OUTPUT_POVRAY_10:
		tab_indent();
		fprintf(outfile, "object { box { <%g %g %g> <%g %g %g> }",
			p1[X], p1[Y], p1[Z], p2[X], p2[Y], p2[Z]);
		if (texture_name != NULL)
		    fprintf(outfile, " texture { %s }", texture_name);
		fprintf(outfile, " }\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_POVRAY_20:
		tab_indent();
		fprintf(outfile, "box { <%g, %g, %g>, <%g, %g, %g>  ",
			p1[X], p1[Y], p1[Z], p2[X], p2[Y], p2[Z]);
		if (texture_name != NULL)
		    fprintf(outfile, " texture { %s }", texture_name);
		fprintf(outfile, " }\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_POLYRAY:
		fprintf(outfile, "object { box <%g, %g, %g>, <%g, %g, %g>",
			p1[X], p1[Y], p1[Z], p2[X], p2[Y], p2[Z]);
		if (texture_name != NULL)
		    fprintf(outfile, " %s", texture_name);
		fprintf(outfile, " }\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_QRT:
		fprintf(outfile, "BEGIN_BBOX\n");
		lib_output_polygon_box(p1, p2);
		fprintf(outfile, "END_BBOX\n");
		break;

	    case OUTPUT_RAYSHADE:
		fprintf(outfile, "box ");
		if (texture_name != NULL)
		    fprintf(outfile, "%s ", texture_name);
		fprintf(outfile, " %g %g %g %g %g %g\n",
			p1[X], p1[Y], p1[Z], p2[X], p2[Y], p2[Z]);
		break;

	    case OUTPUT_ART:
		tab_indent();
		fprintf(outfile, "box {");
		fprintf(outfile, " vertex(%f, %f, %f)\n",
			p1[X], p1[Y], p1[Z]);
		fprintf(outfile, " vertex(%f, %f, %f) }\n",
			p2[X], p2[Y], p2[Z]);
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_RTRACE:
		fprintf(outfile, "2 %d %g %g %g %g %g %g %g\n",
			texture_count, texture_ior,
			(p1[X] + p2[X]) / 2.0,
			(p1[Y] + p2[Y]) / 2.0,
			(p1[Z] + p2[Z]) / 2.0,
			p2[X] - p1[X], p2[Y] - p1[Y], p2[Z] - p1[Z]);
		break;
        }
    }
}


/*-----------------------------------------------------------------*/
/* Given a polygon defined by vertices in verts, determine which of the
   components of the vertex correspond to useful x and y coordinates - with
   these we can pretend the polygon is 2D to do our work on it. */
static void
find_axes(verts)
    COORD3 *verts;
{
    double P1[3], P2[3], x, y, z;

    P1[0] = VERT(1, 0) - VERT(0, 0);
    P1[1] = VERT(1, 1) - VERT(0, 1);
    P1[2] = VERT(1, 2) - VERT(0, 2);

    P2[0] = VERT(2, 0) - VERT(0, 0);
    P2[1] = VERT(2, 1) - VERT(0, 1);
    P2[2] = VERT(2, 2) - VERT(0, 2);

    /* Cross product - don't need to normalize cause we're only interested
      in the size of the components */
    x = fabs(P1[1] * P2[2] - P1[2] * P2[1]);
    y = fabs(P1[2] * P2[0] - P1[0] * P2[2]);
    z = fabs(P1[0] * P2[1] - P1[1] * P2[0]);

    if (x > y && x > z) {
	axis1 = 1;
	axis2 = 2;
    } else if (y > x && y > z) {
	axis1 = 0;
	axis2 = 2;
    } else {
	axis1 = 0;
	axis2 = 1;
    }
}

/*-----------------------------------------------------------------*/
/* Find the left most vertex in the polygon that has vertices m ... n. */
static unsigned int
leftmost_vertex(m, n, verts)
    unsigned int m, n;
    COORD3 *verts;
{
    unsigned int l, i;
    double x;

    /* Assume the first vertex is the farthest to the left */
    l = m;
    x = VERT(m, axis1);

    /* Now see if any of the others are farther to the left */
    for (i=m+1;i<=n;i++) {
       if (VERT(i, axis1) < x) {
	   l = i;
	   x = VERT(i, axis1);
       }
    }
    return l;
}

/*-----------------------------------------------------------------*/
/* Given the leftmost vertex in a polygon, this routine finds another vertex
   can be used to safely split the polygon. */
static unsigned int
split_vertex(l, la, lb, m, n, verts)
    unsigned int l, la, lb, m, n;
    COORD3 *verts;
{
    unsigned int t, k, lpu, lpl;
    double yu, yl;

    yu = MAX(VERT(l, axis2), MAX(VERT(la, axis2), VERT(lb, axis2)));
    yl = MIN(VERT(l, axis2), MIN(VERT(la, axis2), VERT(lb, axis2)));
    if (VERT(lb, axis2) > VERT(la, axis2)) {
	lpu = lb;
	lpl = la;
    } else {
	lpu = la;
	lpl = lb;
    }
    t = (VERT(lb, axis1) > VERT(la, axis1) ? lb : la);
    for (k=m;k<n;k++) {
	if (k != la && k != l && k != lb) {
	    if (VERT(k, axis2) <= yu && VERT(k, axis2) >= yl) {
		if (VERT(k, axis1) < VERT(t, axis1) &&
		    ((VERT(k, axis2) - VERT(l, axis2)) *
		     (VERT(lpu, axis1) - VERT(l, axis1))) <=
		    ((VERT(lpu, axis2) - VERT(l, axis2)) *
		     (VERT(k, axis1) - VERT(l, axis1)))) {
		    if (((VERT(k, axis2) - VERT(l, axis2)) *
			  (VERT(lpl, axis1) - VERT(l, axis1))) >=
			 ((VERT(lpl, axis2) - VERT(l, axis2)) *
			 (VERT(k, axis1) - VERT(l, axis1)))) {
			t = k;
		    }
		}
	    }
	}
    }
    return t;
}

/*-----------------------------------------------------------------*/
/* Test polygon vertices to see if they are linear */
static int
linear_vertices(m, n, verts)
    unsigned int m, n;
    COORD3 *verts;
{
#if defined (applec)
#pragma unused (m,n,verts)
#endif /* applec */
    /* Not doing anything right now */
    return 0;
}

/*-----------------------------------------------------------------*/
/* Shift vertex indices around to make two polygons out of one. */
static void
perform_split(m, m1, n, n1)
    int n, n1, m, m1;
{
    int i, j, k;

    k = n + 3 - m;
    /* Move the new polygon up over the place the current one sits */
    for (j=m1;j<=n1;j++) vbuffer[j+k] = vbuffer[j];

    /* Move top part of remaining polygon */
    for (j=n;j>=n1;j--) vbuffer[j+2] = vbuffer[j];

    /* Move bottom part of remaining polygon */
    k = n1 - m1 + 1;
    for (j=m1;j>=m;j--) vbuffer[j+k] = vbuffer[j];

    /* Copy the new polygon so that it sits before the remaining polygon */
    i = n + 3 - m;
    k = m - m1;
    for (j=m1;j<=n1;j++) vbuffer[j+k] = vbuffer[j+i];
}

/*-----------------------------------------------------------------*/
/* Copy an indirectly referenced triangle into the output triangle buffer */
static void
add_new_triangle(m, verts, norms, out_cnt, out_verts, out_norms)
    unsigned int m, *out_cnt;
    COORD3 *verts, *norms, **out_verts, **out_norms;
{
    if (out_verts != NULL) {
	COPY_COORD3(out_verts[*out_cnt][0], verts[vbuffer[m]]);
	COPY_COORD3(out_verts[*out_cnt][1], verts[vbuffer[m+1]]);
	COPY_COORD3(out_verts[*out_cnt][2], verts[vbuffer[m+2]]);
    }
    if (out_norms != NULL) {
	COPY_COORD3(out_norms[*out_cnt][0], norms[vbuffer[m]]);
	COPY_COORD3(out_norms[*out_cnt][1], norms[vbuffer[m+1]]);
	COPY_COORD3(out_norms[*out_cnt][2], norms[vbuffer[m+2]]);
    }
    *out_cnt += 1;
}

/*-----------------------------------------------------------------*/
static void
split_buffered_polygon(cnt, verts, norms, out_cnt, out_verts, out_norms)
    unsigned int cnt, *out_cnt;
    COORD3 *verts, *norms, **out_verts, **out_norms;
{
    unsigned int i, m, m1, n, n1;
    unsigned int l, la, lb, ls;

    /* No triangles to start with */
    *out_cnt = 0;

    /* Initialize the polygon splitter */
    poly_end[0] = -1;
    poly_end[1] = cnt-1;

    /* Split and push polygons until they turn into triangles */
    for (i=1;i>0;) {
	m = poly_end[i-1] + 1;
	n = poly_end[i];
	if (n - m == 2) {
	    if (!linear_vertices(m, n, verts)) {
		add_new_triangle(m, verts, norms, out_cnt,
			out_verts, out_norms);
	    }
	    i = i - 1;
	} else {
	    l = leftmost_vertex(m, n, verts);
	    la = (l == n ? m : l + 1);
	    lb = (l == m ? n : l - 1);
	    ls = split_vertex(l, la, lb, m, n, verts);
	    if (ls == la || ls == lb) {
		m1 = (la < lb ? la : lb);
		n1 = (la > lb ? la : lb);
	    } else {
		m1 = (l < ls ? l : ls);
		n1 = (l > ls ? l : ls);
	    }
	    perform_split(m, m1, n, n1);
	    poly_end[i++] = m + n1 - m1;
	    poly_end[i] = n + 2;
	}
    }
}

/*-----------------------------------------------------------------*/
/*
 * Split an arbitrary polygon into triangles.
 */
static void
split_polygon(n, vert, norm)
    unsigned int n;
    COORD3 *vert, *norm;
{
    COORD4 tvert[3], v0, v1;
    COORD3 **out_verts, **out_norms;
    unsigned int i, t, out_n;
    object_ptr new_object;

    /* Can't split a NULL vertex list */
    if (vert == NULL) return;
    if (vbuffer == NULL) { /* [esp] Added error */
	lib_initialize();
	/* [are] removed error, go and initialize if it hasn't been done. */
    }

    /* Allocate space to hold the intermediate polygon stacks */
    out_verts = (COORD3 **)malloc((n - 2) * sizeof(COORD3 *));
    if (norm != NULL)
	out_norms = (COORD3 **)malloc((n - 2) * sizeof(COORD3 *));
    else
	out_norms = NULL;
    for (i=0;i<n-2;i++) {
	out_verts[i] = (COORD3 *)malloc(3 * sizeof(COORD3));
	if (norm != NULL)
	    out_norms[i] = (COORD3 *)malloc(3 * sizeof(COORD3));
    }

    /* Start with a strict identity of vertices in verts and vertices in
       the polygon buffer */
    for (i=0;i<n;i++) vbuffer[i] = i;

    /* Make sure we know which axes to look at */
    find_axes(vert);

    out_n = 0;
    split_buffered_polygon(n, vert, norm, &out_n, out_verts, out_norms);

    /* Now output the triangles that we generated */
    for (t=0;t<out_n;t++) {
	PLATFORM_MULTITASK();
	if (format == OUTPUT_DELAYED || format == OUTPUT_PLG) {
	    /* Save all the pertinent information */
	    new_object = (object_ptr)malloc(sizeof(struct object_struct));
	    if (new_object == NULL) return;
	    if (norm == NULL) {
		new_object->object_type  = POLYGON_OBJ;
		new_object->object_data.polygon.tot_vert = 3;
		new_object->object_data.polygon.vert =
			(COORD3 *)malloc(3 * sizeof(COORD3));
		if (new_object->object_data.polygon.vert == NULL) return;
	    } else {
		new_object->object_type  = POLYPATCH_OBJ;
		new_object->object_data.polypatch.tot_vert = 3;
		new_object->object_data.polypatch.vert =
			(COORD3 *)malloc(3 * sizeof(COORD3));
		if (new_object->object_data.polypatch.vert == NULL) return;
		new_object->object_data.polypatch.norm =
			(COORD3 *)malloc(3 * sizeof(COORD3));
		if (new_object->object_data.polypatch.norm == NULL) return;
	    }
	    new_object->curve_format = OUTPUT_PATCHES;
	    new_object->surf_index   = texture_count;
	    for (i=0;i<3;i++) {
		if (norm == NULL)
		    COPY_COORD3(new_object->object_data.polygon.vert[i],
			    out_verts[t][i])
		else {
		    COPY_COORD3(new_object->object_data.polypatch.vert[i],
			    out_verts[t][i])
		    COPY_COORD3(new_object->object_data.polypatch.norm[i],
			    out_norms[t][i])
		}
	    }
	    new_object->next_object = lib_objects;
	    lib_objects = new_object;
	} else {
	    switch (format) {
		case OUTPUT_VIDEO:
		    /* First make sure the display has been opened for
		     * drawing
		     */
		    if (!view_init_flag) {
			lib_create_view_matrix(view.tx, view.from, view.at,
					       view.up, view.resx, view.resy,
					       view.angle, view.aspect);
			display_init(view.resx, view.resy, bk_color);
			view_init_flag = 1;
		    }
		    /* Step through each segment of the polygon, projecting it
		       onto the screen. */
		    for (i=0;i<3;i++) {
			COPY_COORD3(tvert[0], out_verts[t][i]);
			tvert[0][W] = 1.0;
			lib_transform_coord(v0, tvert[0], view.tx);
			COPY_COORD3(tvert[1], out_verts[t][(i+1)%3]);
			tvert[1][W] = 1.0;
			lib_transform_coord(v1, tvert[1], view.tx);
			/* Do the perspective transform on the points */
			v0[X] /= v0[W]; v0[Y] /= v0[W];
			v1[X] /= v1[W]; v1[Y] /= v1[W];
			if (lib_clip_to_box(v0, v1, view_bounds))
			    display_line((int)v0[X], (int)v0[Y],
					(int)v1[X], (int)v1[Y], fg_color);
		    }
		    break;

		case OUTPUT_NFF:
		    if (norm == NULL) {
			fprintf(outfile, "p 3\n");
			for (i=0;i<3;++i) {
			    fprintf(outfile, "%g %g %g\n",
				    out_verts[t][i][X], out_verts[t][i][Y],
				    out_verts[t][i][Z]);
			}
		    } else {
			fprintf(outfile, "pp 3\n");
			for (i=0;i<3;++i) {
			    fprintf(outfile, "%g %g %g %g %g %g\n",
				    out_verts[t][i][X], out_verts[t][i][Y],
				    out_verts[t][i][Z], out_norms[t][i][X],
				    out_norms[t][i][Y], out_norms[t][i][Z]);
			}
		    }
		   break;

		case OUTPUT_POVRAY_10:
		case OUTPUT_POVRAY_20:
		    tab_indent();
		    fprintf(outfile, "object {\n");
		    tab_inc();

		    tab_indent();
		    if (norm == NULL)
			fprintf(outfile, "triangle {\n");
		    else
			fprintf(outfile, "smooth_triangle {\n");
		    tab_inc();

		    for (i=0;i<3;++i) {
			tab_indent();
			if (format == OUTPUT_POVRAY_10) {
			    fprintf(outfile, "<%g %g %g>",
				   out_verts[t][i][X],
				   out_verts[t][i][Y],
				   out_verts[t][i][Z]);
			    if (norm != NULL)
				fprintf(outfile, " <%g %g %g>",
				       out_norms[t][i][X],
				       out_norms[t][i][Y],
				       out_norms[t][i][Z]);
			} else {
			    fprintf(outfile, "<%g, %g, %g>",
				   out_verts[t][i][X],
				   out_verts[t][i][Y],
				   out_verts[t][i][Z]);
			    if (norm != NULL)
			       fprintf(outfile, " <%g, %g, %g>",
				      out_norms[t][i][X],
				      out_norms[t][i][Y],
				      out_norms[t][i][Z]);
			    if (i < 2)
			       fprintf(outfile, ",");
			}
			fprintf(outfile, "\n");
		    } /*for*/

		    tab_dec();
		    tab_indent();
		    fprintf(outfile, "} // tri\n");

		    if (texture_name != NULL) {
		       tab_indent();
		       fprintf(outfile, "texture { %s }\n", texture_name);
		    }

		    tab_dec();
		    tab_indent();
		    fprintf(outfile, "} // object\n");

		    fprintf(outfile, "\n");
		    break;

		case OUTPUT_POLYRAY:
		    if (norm == NULL) {
			tab_indent();
			fprintf(outfile, "object { polygon 3,");
			for (i=0;i<3;i++) {
			    fprintf(outfile, " <%g, %g, %g>",
				    out_verts[t][i][X], out_verts[t][i][Y],
				    out_verts[t][i][Z]);
			    if (i < 2)
				fprintf(outfile, ", ");
			}
		    } else {
			tab_indent();
			fprintf(outfile, "object { patch ");
			for (i=0;i<3;i++) {
			    fprintf(outfile, " <%g, %g, %g>, <%g, %g, %g>",
				    out_verts[t][i][X], out_verts[t][i][Y],
				    out_verts[t][i][Z], out_norms[t][i][X],
				    out_verts[t][i][Y], out_norms[t][i][Z]);
			    if (i < 2)
				fprintf(outfile, ", ");
			}
		    }
		    if (texture_name != NULL)
			fprintf(outfile, " %s", texture_name);
		    fprintf(outfile, " }\n");
		    fprintf(outfile, "\n");
		    break;

		case OUTPUT_VIVID:
		    if (norm == NULL) {
			tab_indent();
			fprintf(outfile, "polygon { points 3 ");
			for (i=0;i<3;i++) {
			    fprintf(outfile, " vertex %g %g %g ",
				    out_verts[t][i][X], out_verts[t][i][Y],
				    out_verts[t][i][Z]);
			}
		    } else {
			fprintf(outfile, "patch {");
			for (i=0;i<3;++i) {
			    fprintf(outfile,
				    " vertex %g %g %g  normal %g %g %g ",
				    out_verts[t][i][X], out_verts[t][i][Y],
				    out_verts[t][i][Z], out_norms[t][i][X],
				    out_norms[t][i][Y], out_norms[t][i][Z]);
			}
		    }
		    fprintf(outfile, " }\n");
		    fprintf(outfile, "\n");
		    break;

		case OUTPUT_QRT:
		    /* Doesn't matter if there are vertex normals,
		     * QRT can't use them.
		     */
		    fprintf(outfile, "TRIANGLE ( ");
		    fprintf(outfile, "loc = (%g, %g, %g), ",
			      out_verts[t][0][X], out_verts[t][0][Y],
			      out_verts[t][0][Z]);
		    fprintf(outfile, "vect1 = (%g, %g, %g), ",
			      out_verts[t][1][X] - out_verts[t][0][X],
			      out_verts[t][1][Y] - out_verts[t][0][Y],
			      out_verts[t][1][Z] - out_verts[t][0][Z]);
		    fprintf(outfile, "vect2 = (%g, %g, %g) ",
			      out_verts[t][2][X] - out_verts[t][0][X],
			      out_verts[t][2][Y] - out_verts[t][0][Y],
			      out_verts[t][2][Z] - out_verts[t][0][Z]);
		    fprintf(outfile, " );\n");
		    break;

		case OUTPUT_RAYSHADE:
		    fprintf(outfile, "triangle ");
		    if (texture_name != NULL)
			fprintf(outfile, "%s ", texture_name);
		    for (i=0;i<3;i++) {
			fprintf(outfile, "%g %g %g ",
				out_verts[t][i][X], out_verts[t][i][Y],
				out_verts[t][i][Z]);
			if (norm != NULL)
			    fprintf(outfile, "%g %g %g ",
				    out_norms[t][i][X], out_norms[t][i][Y],
				    out_norms[t][i][Z]);
		    }
		    fprintf(outfile, "\n");
		    break;

		case OUTPUT_ART:
		    tab_indent();
		    fprintf(outfile, "polygon {\n");
		    tab_inc();

		    tab_indent();
		    for (i=0;i<3;i++) {
			tab_indent();
			fprintf(outfile, "vertex(%f, %f, %f)",
				out_verts[t][i][X], out_verts[t][i][Y],
				out_verts[t][i][Z]);
			if (norm != NULL)
			    fprintf(outfile, ", (%f, %f, %f)\n",
				    out_norms[t][i][X], out_norms[t][i][Y],
				    out_norms[t][i][Z]);
			else
			    fprintf(outfile, "\n");
		    }
		    tab_dec();
		    tab_indent();
		    fprintf(outfile, "}\n");
		    fprintf(outfile, "\n");
		    break;

		case OUTPUT_RTRACE:
		    if (norm == NULL) {
			fprintf(outfile, "5 %d %g 0 0 0 1 1 1 -\n",
				texture_count, texture_ior);
			fprintf(outfile, "3 1 2 3\n\n");
		    } else {
			fprintf(outfile, "6 %d %g 0 0 0 1 1 1 -\n",
				texture_count, texture_ior);
		    }
		    for (i=0;i<3;i++) {
			fprintf(outfile, "%g %g %g",
				out_verts[t][i][X], out_verts[t][i][Y],
				out_verts[t][i][Z]);
			if (norm != NULL) {
			    fprintf(outfile, " %g %g %g",
				    out_norms[t][i][X], out_norms[t][i][Y],
				    out_norms[t][i][Z]);
			}
			fprintf(outfile, "\n");
		    }
		    fprintf(outfile, "\n");
		    break;

		case OUTPUT_RAWTRI:
		    for (i=0;i<3;++i) {
			fprintf(outfile, "%-10.5g %-10.5g %-10.5g  ",
				out_verts[t][i][X], out_verts[t][i][Y],
				out_verts[t][i][Z]);
		    }

#ifdef RAWTRI_WITH_TEXTURES
		    /* raw triangle format extension to do textured raw
		     * triangles */
		    if (texture_name != NULL)
			fprintf(outfile, "%s", texture_name);
		    else
			/* for lack of a better name */
			fprintf(outfile, "texNone");
#endif /* RAWTRI_WITH_TEXTURES */

		    fprintf(outfile, "\n");
		    break;

	    } /* switch */
	} /* else !OUTPUT_DELAYED */
    } /* else for loop */

    /* Clean up intermediate storage */
    for (i=0;i<n-2;i++) {
	free(out_verts[i]);
	if (norm != NULL)
	    free(out_norms[i]);
    }

    free(out_verts);
    if (out_norms != NULL) free(out_norms);
}

/*-----------------------------------------------------------------*/
/*
 * 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]
 *
 */
void
lib_output_polygon(tot_vert, vert)
    int tot_vert;
    COORD3 vert[];
{
    object_ptr new_object;
    int num_vert, i, j;
    COORD3 x;
    COORD4 tvert[3], v0, v1;

    /* First lets do a couple of checks to see if this is a valid polygon */
    for (i=0;i<tot_vert;) {
	/* If there are two adjacent coordinates that degenerate then
	   collapse them down to one */
	SUB3_COORD3(x, vert[i], vert[(i+1)%tot_vert]);
	if (lib_normalize_vector(x) < EPSILON) {
	    for (j=i;j<tot_vert-1;j++)
		memcpy(&vert[j], &vert[j+1], sizeof(COORD4));
	    tot_vert--;
	}
	else {
	    i++;
	}
    }

    if (tot_vert < 3)
	/* No such thing as a poly that only has two sides */
	return;

    if (format == OUTPUT_DELAYED) {
	/* Save all the pertinent information */
	new_object = (object_ptr)malloc(sizeof(struct object_struct));
	new_object->object_data.polygon.vert =
		(COORD3 *)malloc(tot_vert * sizeof(COORD3));
	if (new_object == NULL || new_object->object_data.polygon.vert == NULL)
	    /* Quietly fail */
	    return;
	new_object->object_type  = POLYGON_OBJ;
	new_object->curve_format = OUTPUT_PATCHES;
	new_object->surf_index   = texture_count;
	new_object->object_data.polygon.tot_vert = tot_vert;
	for (i=0;i<tot_vert;i++)
	    COPY_COORD3(new_object->object_data.polygon.vert[i], vert[i])
	new_object->next_object = lib_objects;
	lib_objects = new_object;
    } else {
	switch (format) {
	    case OUTPUT_VIDEO:
		/* First make sure the display has been opened for drawing */
		if (!view_init_flag) {
		    lib_create_view_matrix(view.tx, view.from, view.at,
					   view.up, view.resx, view.resy,
					   view.angle, view.aspect);
		    display_init(view.resx, view.resy, bk_color);
		    view_init_flag = 1;
		}
		/* Step through each segment of the polygon, projecting it
		    onto the screen. */
		for (i=0;i<tot_vert;i++) {
		    COPY_COORD3(tvert[0], vert[i]);
		    tvert[0][W] = 1.0;
		    lib_transform_coord(v0, tvert[0], view.tx);
		    COPY_COORD3(tvert[1], vert[(i+1)%tot_vert]);
		    tvert[1][W] = 1.0;
		    lib_transform_coord(v1, tvert[1], view.tx);
		    /* Do the perspective transform on the points */
		    v0[X] /= v0[W]; v0[Y] /= v0[W];
		    v1[X] /= v1[W]; v1[Y] /= v1[W];
		    if (lib_clip_to_box(v0, v1, view_bounds))
			display_line((int)v0[X], (int)v0[Y],
				(int)v1[X], (int)v1[Y], fg_color);
		}
		break;

	    case OUTPUT_NFF:
		fprintf(outfile, "p %d\n", tot_vert);
		for (num_vert=0;num_vert<tot_vert;++num_vert)
		    fprintf(outfile, "%g %g %g\n",
			    vert[num_vert][X],
			    vert[num_vert][Y],
			    vert[num_vert][Z]);
		break;

	    case OUTPUT_POVRAY_10:
	    case OUTPUT_POVRAY_20:
	    case OUTPUT_QRT:
	    case OUTPUT_PLG:
	    case OUTPUT_RAWTRI:
		/*
		These renderers don't do arbitrary polygons, split the polygon
		into triangles for output
		*/
		split_polygon(tot_vert, vert, NULL);
		break;

	    case OUTPUT_POLYRAY:
		tab_indent();
		fprintf(outfile, "object { polygon %d,", tot_vert);
		for (num_vert = 0; num_vert < tot_vert; num_vert++) {
		    fprintf(outfile, " <%g, %g, %g>",
			     vert[num_vert][X],
			     vert[num_vert][Y],
			     vert[num_vert][Z]);
		    if (num_vert < tot_vert-1)
			fprintf(outfile, ", ");
		}
		if (texture_name != NULL)
		    fprintf(outfile, " %s", texture_name);
		fprintf(outfile, " }\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_VIVID:
		tab_indent();
		fprintf(outfile, "polygon { points %d ", tot_vert);
		for (num_vert = 0; num_vert < tot_vert; num_vert++) {
		    /* Vivid has problems with very long input lines, so in
		     * order to handle polygons with many vertices, we split
		     * the vertices one to a line.
		     */
		    fprintf(outfile, " vertex %g %g %g \n",
			     vert[num_vert][X],
			     vert[num_vert][Y],
			     vert[num_vert][Z]);
		}
		fprintf(outfile, " }\n");
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_RAYSHADE:
		fprintf(outfile, "polygon ");
		if (texture_name != NULL)
		    fprintf(outfile, "%s ", texture_name);
		for (num_vert=0;num_vert<tot_vert;num_vert++) {
		    if (!(num_vert%3)) fprintf(outfile, "\n");
		    fprintf(outfile, "%g %g %g ",
			    vert[num_vert][X],
			    vert[num_vert][Y],
			    vert[num_vert][Z]);
		}
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_RTRACE:
		fprintf(outfile, "5 %d %g 0 0 0 1 1 1 -\n",
			  texture_count, texture_ior);
		fprintf(outfile, "%d ", tot_vert);
		for (num_vert=0;num_vert<tot_vert;num_vert++)
		    fprintf(outfile, "%d ", num_vert+1);
		fprintf(outfile, "\n\n");
		for (num_vert=0;num_vert<tot_vert;num_vert++) {
		    fprintf(outfile, "%g %g %g\n",
			    vert[num_vert][X],
			    vert[num_vert][Y],
			    vert[num_vert][Z]);
		}
		fprintf(outfile, "\n");
		break;

	    case OUTPUT_ART:
		tab_indent();
		fprintf(outfile, "polygon {\n");
		tab_inc();

		for (num_vert=0;num_vert<tot_vert;num_vert++) {
		    tab_indent();
		    fprintf(outfile, "vertex(%f, %f, %f)\n",
			      vert[num_vert][X],
			      vert[num_vert][Y],
			      vert[num_vert][Z]);
		}
		tab_dec();
		tab_indent();
		fprintf(outfile, "}\n");
		fprintf(outfile, "\n");
		break;

	 } /* switch */
    } /* else !OUTPUT_DELAYED */
}


/*-----------------------------------------------------------------*/
/*
 * 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]
 *
 */
void
lib_output_polypatch(tot_vert, vert, norm)
    int tot_vert;
    COORD3 vert[], norm[];
{
    /* As far as I know, there aren't any renderers capable of directly
       generating polygon patches of more than 3 sides.   Therefore we
       will call a routine to split the patch into triangles.
    */
    split_polygon(tot_vert, vert, norm);
}

/*-----------------------------------------------------------------*/
static void
dump_plg_file()
{
    object_ptr temp_obj;
    int i;
    unsigned int fcnt, vcnt;

    fcnt = 0;
    vcnt = 0;
    for (temp_obj = polygon_stack;
	 temp_obj != NULL;
	 temp_obj = temp_obj->next_object) {
	fcnt++;
	vcnt += temp_obj->object_data.polygon.tot_vert;
    }

    fprintf(outfile, "objx %d %d\n", vcnt, fcnt);

    /* Dump all vertices */
    for (temp_obj = polygon_stack;
	 temp_obj != NULL;
	 temp_obj = temp_obj->next_object) {

	PLATFORM_MULTITASK();
	for (i=0;i<temp_obj->object_data.polygon.tot_vert;i++) {
	    fprintf(outfile, "%g %g %g\n",
		   temp_obj->object_data.polygon.vert[i][X],
		   temp_obj->object_data.polygon.vert[i][Y],
		   temp_obj->object_data.polygon.vert[i][Z]);
	}
    }

    /* Dump all faces */
    vcnt = 0;
    for (temp_obj = polygon_stack;
	 temp_obj != NULL;
	 temp_obj = temp_obj->next_object) {

	PLATFORM_MULTITASK();
	fprintf(outfile, "0x11ff %d ", temp_obj->object_data.polygon.tot_vert);
	for (i=0;i<temp_obj->object_data.polygon.tot_vert;i++)
	    fprintf(outfile, "%d ", vcnt + i);
	fprintf(outfile, "\n");
	vcnt += i;
    }
}

/*-----------------------------------------------------------------*/
static void
dump_all_objects()
{
    object_ptr temp_obj;

    if (format == OUTPUT_RTRACE)
	fprintf(outfile, "Objects\n");

    /* Step through all objects dumping them as we go. */
    for (temp_obj = lib_objects;
	 temp_obj != NULL;
	 temp_obj = temp_obj->next_object) {

	PLATFORM_MULTITASK();
	lookup_surface_stats(temp_obj->surf_index,
		&texture_count, &texture_ior);
	switch (temp_obj->object_type) {
	    case BOX_OBJ:
		lib_output_box(temp_obj->object_data.box.point1,
			       temp_obj->object_data.box.point2);
		break;
	    case CONE_OBJ:
		lib_output_cylcone(temp_obj->object_data.cone.base_pt,
				   temp_obj->object_data.cone.apex_pt,
				   temp_obj->curve_format);
		break;
	    case DISC_OBJ:
		lib_output_disc(temp_obj->object_data.disc.center,
				temp_obj->object_data.disc.normal,
				temp_obj->object_data.disc.iradius,
				temp_obj->object_data.disc.oradius,
				temp_obj->curve_format);
		break;
	    case HEIGHT_OBJ:
		lib_output_height(temp_obj->object_data.height.filename,
				  temp_obj->object_data.height.data,
				  temp_obj->object_data.height.height,
				  temp_obj->object_data.height.width,
				  temp_obj->object_data.height.x0,
				  temp_obj->object_data.height.x1,
				  temp_obj->object_data.height.y0,
				  temp_obj->object_data.height.y1,
				  temp_obj->object_data.height.z0,
				  temp_obj->object_data.height.z1);
		break;
	    case POLYGON_OBJ:
		lib_output_polygon(temp_obj->object_data.polygon.tot_vert,
				   temp_obj->object_data.polygon.vert);
		break;
	    case POLYPATCH_OBJ:
		lib_output_polypatch(temp_obj->object_data.polypatch.tot_vert,
				     temp_obj->object_data.polypatch.vert,
				     temp_obj->object_data.polypatch.norm);
		break;
	    case SPHERE_OBJ:
		lib_output_sphere(temp_obj->object_data.sphere.center_pt,
				  temp_obj->curve_format);
		break;
	    case SUPERQ_OBJ:
		lib_output_sq_sphere(temp_obj->object_data.superq.center_pt,
				     temp_obj->object_data.superq.a1,
				     temp_obj->object_data.superq.a2,
				     temp_obj->object_data.superq.a3,
				     temp_obj->object_data.superq.n,
				     temp_obj->object_data.superq.e);
		break;
	    case TORUS_OBJ:
		lib_output_torus(temp_obj->object_data.torus.center,
				 temp_obj->object_data.torus.normal,
				 temp_obj->object_data.torus.iradius,
				 temp_obj->object_data.torus.oradius,
				 temp_obj->curve_format);
		break;
	    default:
		fprintf(outfile, "Bad object type: %d\n",
			  temp_obj->object_type);
		exit(1);
	}
    }

    if (format == OUTPUT_RTRACE)
	fprintf(outfile, "\n");
}

/*-----------------------------------------------------------------*/
static void
reorder_surfaces()
{
    surface_ptr temp_ptr, head_ptr = NULL;

    while (lib_surfaces != NULL) {
	temp_ptr = lib_surfaces;
	lib_surfaces = lib_surfaces->next;
	temp_ptr->next = head_ptr;
	head_ptr = temp_ptr;
    }

    lib_surfaces = head_ptr;
}

/*-----------------------------------------------------------------*/
static void
dump_all_lights()
{
    light_ptr temp_ptr = lib_lights;

    if (format == OUTPUT_RTRACE)
	fprintf(outfile, "Lights\n");

    while (temp_ptr != NULL) {
	lib_output_light(temp_ptr->center_pt);
	temp_ptr = temp_ptr->next;
    }

    if (format == OUTPUT_RTRACE)
	fprintf(outfile, "\n");
}

/*-----------------------------------------------------------------*/
static void
dump_all_surfaces()
{
    surface_ptr temp_ptr = lib_surfaces;

    if (format == OUTPUT_RTRACE)
	fprintf(outfile, "Surfaces\n");

    while (temp_ptr != NULL) {
	lib_output_color(temp_ptr->surf_name, temp_ptr->color, temp_ptr->ka,
			      temp_ptr->kd, temp_ptr->ks, temp_ptr->shine,
			      temp_ptr->ang, temp_ptr->kt, temp_ptr->ior);
	temp_ptr = temp_ptr->next;
    }

    if (format == OUTPUT_RTRACE)
	fprintf(outfile, "\n");
}

/*-----------------------------------------------------------------*/
void
lib_clear_database()
{
    surface_ptr ts1, ts2;
    object_ptr to1, to2;
    light_ptr tl1, tl2;

    outfile = stdout;
    texture_name = NULL;
    texture_count = 0;
    texture_ior = 1.0;
    format = OUTPUT_POVRAY_10;
    u_resolution = OUTPUT_RESOLUTION;
    v_resolution = OUTPUT_RESOLUTION;
    SET_COORD3(bk_color, 0.0, 0.0, 0.0);
    SET_COORD3(fg_color, 0.0, 0.0, 0.0);

    /* Remove all surfaces */
    ts1 = lib_surfaces;
    while (ts1 != NULL) {
	ts2 = ts1;
	ts1 = ts1->next;
	free(ts2->surf_name);
	free(ts2);
    }
    lib_surfaces = NULL;

    /* Remove all objects */
    to1 = lib_objects;
    while (to1 != NULL) {
	to2 = to1;
	to1 = to1->next_object;
	free(to2);
    }
    lib_objects = NULL;

    /* Remove all lights */
    tl1 = lib_lights;
    while (tl1 != NULL) {
	tl2 = tl1;
	tl1 = tl1->next;
	free(tl2);
    }
    lib_lights = NULL;

    /* Reset the view */

    /* Deallocate polygon buffer */
    if (vbuffer != NULL)
	lib_shutdown();
}

/*-----------------------------------------------------------------*/
void
lib_flush_definitions()
{
    unsigned int cnt = 0;

    switch (format) {
       case OUTPUT_RTRACE:
       case OUTPUT_VIDEO:
       case OUTPUT_NFF:
       case OUTPUT_POVRAY_10:
       case OUTPUT_POLYRAY:
       case OUTPUT_VIVID:
       case OUTPUT_QRT:
       case OUTPUT_RAYSHADE:
       case OUTPUT_POVRAY_20:
       case OUTPUT_PLG:
       case OUTPUT_RAWTRI:
	    lib_output_viewpoint(view.from, view.at, view.up, view.angle,
		    view.aspect, view.hither, view.resx, view.resy);

	    lib_output_background_color(bk_color);

	    dump_all_lights();

	    reorder_surfaces();

	    dump_all_surfaces();

	    dump_all_objects();

	    if (format == OUTPUT_RTRACE)
		fprintf(outfile, "Textures\n\n");

	    break;
       case OUTPUT_DELAYED:
	    fprintf(stderr, "Error: Renderer not selected before flushing\n");
	    exit(1);
    }

    if (format == OUTPUT_PLG) {
	/* An extra step is needed to build the polygon file. */
	dump_plg_file();
    }
}

/* End of file */