/*
 * pipes.c - pipe network geometry creation program
 *
 * Mark J. Stock
 * mstock@umich.edu
 *
 * v0.1  1998-08-04  First working version, rewrite from Java
 * v0.2  1998-08-10  Allow curved secions, Radiance output
 * v0.3  2001-12-14  Allow wrapping around edges, better input, output, usage
 * v0.4  2001-12-16  Cylindrical and spherical mappings added
 * v0.5  2005-03-03  Dynamic memory allocation, bug fixes
 * v0.6  2005-03-26  PNG output, better command-line handling, output
 *
 * compile with "cc -o pipes pipes.c -O2 -lm -lpng"
 *
 * Pipes was originally a Java program to create modelling code for
 * a series of non-intersecting, regular-grid pipes. It was slow
 * and I didn't care to support it. Plus, I was a poor programmer
 * back then. I'm not too much better now. Can you tell?
 *
 * To do:
 *   write octree instances for each section of pipe, cartesian only
 *   optionally write text file of center points of pipe path (for post-processing)
 *   write location-specific texture info
 *   use a probability field of values to determine next move, must
 *     define the field somehow...start with a sphere, move to a 
 *     file read-in thingy...brick of bytes, or something
 *   allow branching, track total mass flow to affect minor radius
 *   force pipe to end at valid sides, not just when it can't go
 *     anywhere
 *   read in and write out connection points along the sides
 *   add option for input of random seed
 *
 *   allow wrapping around sides - DONE
 *   write geometry in cylindrical, spherical coordinates - DONE
 *
 * Copyright 1998,2001,2005 Mark J. Stock
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.

 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software Foundation,
 * Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 
 */


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include "png.h"

#define TRUE 1
#define FALSE 0
#define BOOL unsigned char

// these are the maximum allowed dimensions of the grid array
#define MAX_X 5000
#define MAX_Y 5000
#define MAX_Z 5000

// some global output options
int write_rad = FALSE;		// radiance polygons/cylinders
int write_obj3 = FALSE;		// wavefront .obj triangles (uncompressed)
int write_obj2 = FALSE;		// .obj-like segments
int num_v = 0;			// vertex counter

// a 3-tuple, kind of a lame definition, too, huh?
typedef struct vector_record {
   double x;
   double y;
   double z;
} VEC;

// just in case we want to do weird stuff with the coordinate type
// oh, I should say, the "size" value CAN BE NEGATIVE, allow for it.
// #define CELL_COORD_TYPE unsigned char
#define CELL_COORD_TYPE int

// a cell_record is a section of pipe
typedef struct cell_record *cell_ptr;
typedef struct cell_record {
   CELL_COORD_TYPE x;
   CELL_COORD_TYPE y;
   CELL_COORD_TYPE z;
   cell_ptr previous;
} CELL;
CELL size = {0, 0, 0, NULL};
cell_ptr pipe_tail;

// anything we want to save on a side is here
typedef struct side_record {
   int x_low;
   int x_high;
   int y_low;
   int y_high;
   int z_low;
   int z_high;
} SIDE;

/* Validity of each side for pipe start */
SIDE valid_start = {FALSE, FALSE, FALSE, TRUE, FALSE, FALSE};

/* Area of each valid starting side */
SIDE valid_area = {0, 0, 0, 0, 0, 0};

/* Validity of each side for pipe start, use first three only! */
SIDE wrap_side = {FALSE, 0, FALSE, 0, FALSE, 0};

// The cell record
// FALSE means the cell is taken, TRUE means it is free 
// BOOL not_taken[MAX_X][MAX_Y][MAX_Z];
BOOL ***not_taken;

/* the world size of the objects */
// double cell_size = 1.0;	/* cell size */
VEC world = {1.0, 1.0, 1.0};	/* cell size, cartesian world dimensions */

// pipe parameters
double r1 = 0.5;		/* major radius */
double r2;			/* minor radius */


// Alternative coordinates (cartesian is default)

// For cylindrical coords, x becomes radius, y becomes theta,
//   cyl_data[0] is radial offset
//   cyl_data[1] is number of copies to make one revolution
BOOL use_cylindrical = FALSE;
double cyl_data[2] = {1.0, 1.0};

// For spherical coords, x becomes radius, y becomes theta, z becomes phi
//   sph_data[0] is radial offset
//   sph_data[1] is number of copies to make one revolution
//   sph_data[2] is latitude in degrees from equator of lowest z pipe
//   sph_data[3] is latitude in degrees from equator of highest z pipe
BOOL use_spherical = FALSE;
double sph_data[4] = {1.0, 1.0, -80.0, 80.0};

// and some defines...why?
void transform_coord(VEC*);
void find_matrix(CELL*,CELL*,double**);
void find_midpoint(CELL*,CELL*,VEC*);
void find_midpoint2(CELL*,CELL*,double*);
void find_center(CELL*,VEC*);
void find_center2(CELL*,double*);
void write_curve(FILE*,VEC*,VEC*,VEC*,int);
void write_smooth_curve(FILE*,VEC*,VEC*,VEC*,int,int);
void write_smooth_straight(FILE*,VEC*,VEC*,int,int);
void extrude_square_straight(FILE*,int,int,double*,double**,double*,double**);
void extrude_square_end(FILE*,int,int,double*,double**);
void extrude_square_curve(FILE*,int,int,double*,double**,double*,double**);
void write_polygon_4(FILE*,int,int,double*,double*,double*,double*);
void vec_cross(double*,double*i,double*);
int write_pgm_image(char*,int,int,float**);
int write_png_image(char*,int,int,float**);
int** allocate_2d_array_i(int,int);
float** allocate_2d_array_f(int,int);
double** allocate_2d_array_d(int,int);
png_byte** allocate_2d_array_pb(int,int,int);
BOOL*** allocate_3d_array_B(int,int,int);


// let's get it on.
int main(int argc, char* argv[]){

   int i, j, k, count, temp_i;
   char progname[80];
   int seed = 1;
   int num_tries = 0;
   int pipe_count = 0;
   int seg_count = 0;
   int min_length = 10;
   int max_length = 10000;
   int max_pipes = 1;
   int curve_type = -1;	// if 0, use gensurf with smooth polys
                        // if 1, use spheres at corners
                        // if -1, use polygons to make square cross-section
                        // if other, use genworm with number of segments
                        // equal to value of curve_detail, default is 5
   double straight_frac = 0.7;
   char previewoutfile[80];
   int write_pgm = FALSE;
   int write_png = FALSE;
   char geomoutfile[80];
   BOOL temp_b;
   VEC temp_v;
   CELL temp_l;
   cell_ptr current,next_seg;
   FILE *outfp;

   void set_areas();
   void mark_cell_as_occupied(CELL *);
   int choose_next_cell(CELL *,CELL *,double);
   void write_pipe_geom(CELL *,FILE *,int,int);
   void write_data_preview(int,int,char*,int);
   void free_mem(CELL *);

   // the default solution space, voxel units
   size.x = 20;
   size.y = 20;
   size.z = 5;

   r2 = 0.2;

   /* Parse command-line args */
   (void) strcpy(progname,argv[0]);
   if (argc < 3) (void) Usage(progname,0);
   for (i=1; i<argc; i++) {
      if (strncmp(argv[i], "-sxl",4) == 0) {
         valid_start.x_low = TRUE;
      } else if (strncmp(argv[i], "-sxh",4) == 0) {
         valid_start.x_high = TRUE;
      } else if (strncmp(argv[i], "-syl",4) == 0) {
         valid_start.y_low = TRUE;
      } else if (strncmp(argv[i], "-syh",4) == 0) {
         valid_start.y_high = TRUE;
      } else if (strncmp(argv[i], "-szl",4) == 0) {
         valid_start.z_low = TRUE;
      } else if (strncmp(argv[i], "-szh",4) == 0) {
         valid_start.z_high = TRUE;
      } else if (strncmp(argv[i], "-wx", 3) == 0) {
         wrap_side.x_low = TRUE;
      } else if (strncmp(argv[i], "-wy", 3) == 0) {
         wrap_side.y_low = TRUE;
      } else if (strncmp(argv[i], "-wz", 3) == 0) {
         wrap_side.z_low = TRUE;
      } else if (strncmp(argv[i], "-rad", 4) == 0) {
         write_rad = TRUE;
         write_obj3 = FALSE;
         if (argc == i+1) {
            (void) sprintf(geomoutfile,"out.rad");
         } else if (isalnum(argv[i+1][0]) || argv[i+1][0]=='/' || argv[i+1][0]=='.') {
            (void) strcpy(geomoutfile,argv[++i]);
         } else {
            (void) sprintf(geomoutfile,"out.rad");
         }
      } else if (strncmp(argv[i], "-obj", 2) == 0) {
         write_obj3 = TRUE;
         write_rad = FALSE;
         if (argc == i+1) {
            (void) sprintf(geomoutfile,"out.obj");
         } else if (isalnum(argv[i+1][0]) || argv[i+1][0]=='/' || argv[i+1][0]=='.') {
            (void) strcpy(geomoutfile,argv[++i]);
         } else {
            (void) sprintf(geomoutfile,"out.obj");
         }
      } else if (strncmp(argv[i], "-pgm", 3) == 0) {
         write_pgm = TRUE;
         write_png = FALSE;
         if (argc == i+1) {
            (void) sprintf(previewoutfile,"out.pgm");
         } else if (isalnum(argv[i+1][0]) || argv[i+1][0]=='/' || argv[i+1][0]=='.') {
            (void) strcpy(previewoutfile,argv[++i]);
         } else {
            (void) sprintf(previewoutfile,"out.pgm");
         }
      } else if (strncmp(argv[i], "-png", 3) == 0) {
         write_png = TRUE;
         write_pgm = FALSE;
         if (argc == i+1) {
            (void) sprintf(previewoutfile,"out.png");
         } else if (isalnum(argv[i+1][0]) || argv[i+1][0]=='/' || argv[i+1][0]=='.') {
            (void) strcpy(previewoutfile,argv[++i]);
         } else {
            (void) sprintf(previewoutfile,"out.png");
         }
      } else if (strncmp(argv[i], "-cyl", 4) == 0) {
         cyl_data[0] = atof(argv[++i]);
         cyl_data[1] = atof(argv[++i]);
         use_cylindrical = TRUE;
      } else if (strncmp(argv[i], "-sph", 4) == 0) {
         sph_data[0] = atof(argv[++i]);
         sph_data[1] = atof(argv[++i]);
         sph_data[2] = atof(argv[++i]);
         sph_data[3] = atof(argv[++i]);
         use_spherical = TRUE;
      } else if (strncmp(argv[i], "-sq", 3) == 0) {
         curve_type = -1;
      } else if (strncmp(argv[i], "-gs", 3) == 0) {
         curve_type = 0;
      } else if (strncmp(argv[i], "-worm", 2) == 0) {
         curve_type = atoi(argv[++i]); 
         if (curve_type < 2) curve_type = 2;
      } else if (strncmp(argv[i], "-seed", 5) == 0) {
         seed = atoi(argv[++i]);
      } else if (strncmp(argv[i], "-s", 2) == 0) {
         size.x = atoi(argv[++i]); 
         size.y = atoi(argv[++i]); 
         size.z = atoi(argv[++i]);
      } else if (strncmp(argv[i], "-n", 2) == 0) {
         max_pipes = atoi(argv[++i]); 
      } else if (strncmp(argv[i], "-l", 2) == 0) {
         min_length = atoi(argv[++i]); 
         max_length = atoi(argv[++i]);
      } else if (strncmp(argv[i], "-f", 2) == 0) {
         straight_frac = atof(argv[++i]);
      } else if (strncmp(argv[i], "-r", 2) == 0) {
         r2 = atof(argv[++i]);
      } else {
         (void) Usage(progname,0);
      }
   }

   /* Check input and act on user-defined options */

   // initialize random number generator
   srand(seed);

   // make sure sizes are within memory bounds
   if (size.x > MAX_X || size.y > MAX_Y || size.z > MAX_Z) {
      fprintf(stderr,"ERROR: Maximum dimensions are %d %d %d (in x,y,z)\n",MAX_X,MAX_Y,MAX_Z);
      fprintf(stderr,"exiting.\n");
      exit(0);
   }
   if (size.x*size.y*size.z > 500000000) {
      fprintf(stderr,"ERROR: Requested dimensions (%d %d %d) would require\n",size.x,size.y,size.z);
      fprintf(stderr,"   %d MB of RAM.\n",size.x*size.y*size.z*sizeof(BOOL));
      fprintf(stderr,"exiting.\n");
      exit(0);
   }

   // Wavefront .obj files can only be written for square cross-sections
   if (write_obj3 && curve_type != -1) {
      fprintf(stderr,"Warning: Output geometry is Wavefront (.obj), but curve\n");
      fprintf(stderr,"  type was not square, changing to square cross-section.\n");
      curve_type = -1;
   }

   if (!write_obj3 && !write_rad) {
      fprintf(stderr,"Warning: No output geometry to be written. Use -obj\n");
      fprintf(stderr,"  or -rad to write geometry\n");
   }

   // allocate space for the "taken" matrix
   not_taken = allocate_3d_array_B(size.x,size.y,size.z);

   // make sure that valid_start sides are not wrap_sides
   if (wrap_side.x_low && (valid_start.x_low || valid_start.x_high)) {
      fprintf(stderr,"ERROR: Cannot start pipes on wrapped edge (x edge)\n");
      fprintf(stderr,"exiting.\n");
      exit(0);
   }
   if (wrap_side.y_low && (valid_start.y_low || valid_start.y_high)) {
      fprintf(stderr,"ERROR: Cannot start pipes on wrapped edge (y edge)\n");
      fprintf(stderr,"exiting.\n");
      exit(0);
   }
   if (wrap_side.z_low && (valid_start.z_low || valid_start.z_high)) {
      fprintf(stderr,"ERROR: Cannot start pipes on wrapped edge (z edge)\n");
      fprintf(stderr,"exiting.\n");
      exit(0);
   }

   // set the valid_area of the valid_start sides
   (void) set_areas();

   // clean the solution space
   for (i=0; i<size.x; i++) {
      for (j=0; j<size.y; j++) {
         for (k=0; k<size.z; k++) {
            not_taken[i][j][k] = TRUE;
   }  }  }

   /* Open the obj or rad output file */
   if (write_rad || write_obj3) {
      outfp = fopen(geomoutfile,"w");
      if (outfp==NULL) {
         fprintf(stderr,"Could not open output file %s\n",geomoutfile);
         fprintf(stderr,"exiting.\n");
         exit(0);
      }
   }

   /* Loop until there is no space left, or all pipes have been drawn */
   while (pipe_count < max_pipes) {

      /* determine starting point */
      current = (CELL *)malloc(sizeof(CELL));
      /* fprintf(stderr,"\nStarting new pipe\n"); */
      if (find_start_cell(current) == 0) {
         fprintf(stderr,"Could not find valid starting point, stopping after %d pipes.\n",pipe_count);
         free(current);
         break;
      } else {
         // starting point is valid, occupy it and begin
         (void) mark_cell_as_occupied(current);
         seg_count = 1;
      }

      /* determine pipe path: loop while pipe does not end */
      while (seg_count < max_length) {

         /* determine which way the segment goes */
         next_seg = (CELL *)malloc(sizeof(CELL));
         if (choose_next_cell(current,next_seg,straight_frac) == FALSE) {
            free(next_seg);
            break;
         }

         /* add the cell to the linked list for the pipe */
         next_seg->previous = current;

         /* mark the new current cell as used */
         mark_cell_as_occupied(next_seg);

         /* reset current */
         current = next_seg;

         seg_count++;
      }

      /* if it's big enough, write it, start at current and work back */
      if (seg_count >= min_length) {
         fprintf(stderr,"Pipe %d is %d segments long",pipe_count,seg_count);
         if (write_rad || write_obj3) {
            fprintf(stderr,", writing...\n");
            write_pipe_geom(current,outfp,pipe_count,curve_type);
         } else {
            fprintf(stderr,"\n");
         }
         free_mem(current);
         num_tries = 0;
      } else {
         free_mem(current);
         if (num_tries++ > 100) break;
         continue;
      }

      pipe_count++;
   }


   // write preview image
   if (write_pgm || write_png)
      write_data_preview(write_pgm,write_png,previewoutfile,0);

   // take me home.
   exit(0);
}


/*
 * Calculate the areas of each valid starting side
 */
void set_areas() {

   if (valid_start.x_low)
      valid_area.x_low = size.y*size.z;
   if (valid_start.x_high)
      valid_area.x_high = size.y*size.z;
   if (valid_start.y_low)
      valid_area.y_low = size.x*size.z;
   if (valid_start.y_high)
      valid_area.y_high = size.x*size.z;
   if (valid_start.z_low)
      valid_area.z_low = size.y*size.x;
   if (valid_start.z_high)
      valid_area.z_high = size.y*size.x;

   // fprintf(stderr,"Areas are %d %d %d %d %d %d\n",valid_area.x_low,valid_area.x_high,
   //       valid_area.y_low,valid_area.y_high,valid_area.z_low,valid_area.z_high);
}


/*
 * Determine an acceptable starting cell: must be on a valid
 * boundary (per valid_start), and must not be taken
 */
int find_start_cell(CELL *startcell) {

   int area_sum, area_tag;
   int num_tries = 0;

   // being the first cell, it has no previous cell
   startcell->previous = NULL;

   while (num_tries < 100) {

      // first, choose which side to start on
      area_sum = valid_area.x_low + valid_area.x_high + valid_area.y_low + 
         valid_area.y_high + valid_area.z_low + valid_area.z_high;
      area_tag = (int)((double)area_sum * rand()/(RAND_MAX+1.0));
      // fprintf(stderr,"starting cell...%d / %d\n",area_tag,area_sum);
      // fprintf(stderr,"  valid_area.x_[l|h] is %d %d\n",valid_area.x_low,valid_area.x_high);
      // fprintf(stderr,"  valid_area.y_[l|h] is %d %d\n",valid_area.y_low,valid_area.y_high);
      // fprintf(stderr,"  valid_area.z_[l|h] is %d %d\n",valid_area.z_low,valid_area.z_high);

      // then, choose a starting cell based on that side
      if (area_tag < valid_area.x_low) {
         startcell->x = 0;
         startcell->y = area_tag/size.z;
         startcell->z = area_tag%size.z;
      } else if ((area_tag -= valid_area.x_low) < valid_area.x_high) {
         startcell->x = size.x-1;
         startcell->y = area_tag/size.z;
         startcell->z = area_tag%size.z;
      } else if ((area_tag -= valid_area.x_high) < valid_area.y_low) {
         startcell->x = area_tag%size.x;
         startcell->y = 0;
         startcell->z = area_tag/size.x;
      } else if ((area_tag -= valid_area.y_low) < valid_area.y_high) {
         // fprintf(stderr,"area_tag %g, size.x %g, size.y %g\n",area_tag,size.x,size.y);
         startcell->x = area_tag%size.x;
         startcell->y = size.y-1;
         startcell->z = area_tag/size.x;
      } else if ((area_tag -= valid_area.y_high) < valid_area.z_low) {
         startcell->x = area_tag%size.x;
         startcell->y = area_tag/size.x;
         startcell->z = 0;
      } else if ((area_tag -= valid_area.z_low) < valid_area.z_high) {
         startcell->x = area_tag%size.x;
         startcell->y = area_tag/size.x;
         startcell->z = size.z-1;
      }

      // fprintf(stderr,"area_tag is %d, loc is %d %d %d\n",area_tag,startcell->x,startcell->y,startcell->z);

      // return if succesful
      if (not_taken[startcell->x][startcell->y][startcell->z]) return(1);

      // if not, try another one
      num_tries++;
   }

   /* return unsuccesful */
   return(0);
}


/*
 * Mark this cell as taken, no other pipe can go thru it
 */
void mark_cell_as_occupied(CELL *cell) {

   not_taken[cell->x][cell->y][cell->z] = FALSE;
   // fprintf(stderr,"   marked cell %d %d %d as occupied\n",cell->x,cell->y,cell->z);

}


/*
 * This one is pretty complicated: cell is a pointer to the
 * current end segment. new_cell is a pointer to the new
 * segment, frac is the "straight fraction", the probability
 * that the next segment will just go straight.
 */
int choose_next_cell(CELL *cell, CELL *new_cell, double frac) {

   int ct = 0;		/* number of turn cells within bounds and unoccupied, max=4 */
   int choose;		/* which turn cell is next */
   int temp_x, temp_y, temp_z;
   CELL straight, turn[4];
   int go_straight = FALSE;
   int go_turn[4] = {FALSE, FALSE, FALSE, FALSE};

   // fprintf(stderr,"Find next cell\n");

   /* find all valid next cells, start with the straight option */
   if (cell->previous == NULL) {
      /* we're finding the 2nd cell, straight is away from the wall */
      // fprintf(stderr,"   find 2nd valid straight cell\n");
      if (cell->x == 0) {
         straight.x = 1;
         straight.y = cell->y;
         straight.z = cell->z;
      } else if (cell->x == size.x-1) {
         straight.x = size.x-2;
         straight.y = cell->y;
         straight.z = cell->z;
      } else if (cell->y == 0) {
         straight.x = cell->x;
         straight.y = 1;
         straight.z = cell->z;
      } else if (cell->y == size.y-1) {
         straight.x = cell->x;
         straight.y = size.y-2;
         straight.z = cell->z;
      } else if (cell->z == 0) {
         straight.x = cell->x;
         straight.y = cell->y;
         straight.z = 1;
      } else if (cell->z == size.z-1) {
         straight.x = cell->x;
         straight.y = cell->y;
         straight.z = size.z-2;
      }
      // is the straight-ahead cell in-bounds?
      if (is_cell_in_bounds(straight.x,straight.y,straight.z)) {
         // and is that cell occupied or not?
         go_straight = not_taken[straight.x][straight.y][straight.z];
      }

   } else {
      /* we're finding the next cell, straight is away from previous */
      // fprintf(stderr,"   find next valid straight cell, this= %d %d %d, last= %d %d %d\n",
      //    cell->x,cell->y,cell->z,cell->previous->x,cell->previous->y,cell->previous->z);
      temp_x = 2*cell->x - cell->previous->x;
      temp_y = 2*cell->y - cell->previous->y;
      temp_z = 2*cell->z - cell->previous->z;
      go_straight = FALSE;
      if (is_cell_in_bounds(temp_x,temp_y,temp_z)) {
         straight.x = (temp_x+size.x)%size.x;
         straight.y = (temp_y+size.y)%size.y;
         straight.z = (temp_z+size.z)%size.z;
         go_straight = not_taken[straight.x][straight.y][straight.z];
      }
   }
   // fprintf(stderr,"      found it.\n");

   /* now, determine which of the turns are possible */
   if (cell->previous == NULL) {
      /* we're finding the 2nd cell, straight is away from the wall */
      /* do the 4 possible cells around the first cell */
      /* fprintf(stderr,"   find 4 turn cells off first cell\n"); */

      if (wrap_side.z_low || (cell->z > 0 && cell->z < size.z-1)) {
         /* 2 possible cells are +-1 in z-direction */
         turn[ct].x = (cell->x+size.x)%size.x;
         turn[ct].y = (cell->y+size.y)%size.y;
         turn[ct].z = (cell->z + 1+size.z)%size.z;
         if (not_taken[turn[ct].x][turn[ct].y][turn[ct].z]) ct++;
         turn[ct].x = (cell->x+size.x)%size.x;
         turn[ct].y = (cell->y+size.y)%size.y;
         turn[ct].z = (cell->z - 1+size.z)%size.z;
         if (not_taken[turn[ct].x][turn[ct].y][turn[ct].z]) ct++;
      }
      if (wrap_side.y_low || (cell->y > 0 && cell->y < size.y-1)) {
         /* 2 possible cells are +-1 in y-direction */
         turn[ct].x = (cell->x+size.x)%size.x;
         turn[ct].y = (cell->y + 1+size.y)%size.y;
         turn[ct].z = (cell->z+size.z)%size.z;
         if (not_taken[turn[ct].x][turn[ct].y][turn[ct].z]) ct++;
         turn[ct].x = (cell->x+size.x)%size.x;
         turn[ct].y = (cell->y - 1+size.y)%size.y;
         turn[ct].z = (cell->z+size.z)%size.z;
         if (not_taken[turn[ct].x][turn[ct].y][turn[ct].z]) ct++;
      }
      if (wrap_side.x_low || (cell->x > 0 && cell->x < size.x-1)) {
         /* 2 possible cells are +-1 in x-direction */
         turn[ct].x = (cell->x + 1+size.x)%size.x;
         turn[ct].y = (cell->y+size.y)%size.y;
         turn[ct].z = (cell->z+size.z)%size.z;
         if (not_taken[turn[ct].x][turn[ct].y][turn[ct].z]) ct++;
         turn[ct].x = (cell->x - 1+size.x)%size.x;
         turn[ct].y = (cell->y+size.y)%size.y;
         turn[ct].z = (cell->z+size.z)%size.z;
         if (not_taken[turn[ct].x][turn[ct].y][turn[ct].z]) ct++;
      }
      /* fprintf(stderr,"   %d turn cells in bounds and unoccupied\n",ct); */

   } else {

      /* find the 4 arbitrary possible cells */
      // fprintf(stderr,"   find 4 cells off current cell\n");

      if (cell->x - cell->previous->x == 0) {
         /* 2 possible cells are +-1 in x-direction */
         if (is_cell_in_bounds(cell->x + 1,cell->y,cell->z)) {
            turn[ct].x = (cell->x + 1+size.x)%size.x;
            turn[ct].y = (cell->y+size.y)%size.y;
            turn[ct].z = (cell->z+size.z)%size.z;
            if (not_taken[turn[ct].x][turn[ct].y][turn[ct].z]) ct++;
         }
         if (is_cell_in_bounds(cell->x - 1,cell->y,cell->z)) {
            turn[ct].x = (cell->x - 1+size.x)%size.x;
            turn[ct].y = (cell->y+size.y)%size.y;
            turn[ct].z = (cell->z+size.z)%size.z;
            if (not_taken[turn[ct].x][turn[ct].y][turn[ct].z]) ct++;
         }
      }
      if (cell->y - cell->previous->y == 0) {
         /* 2 possible cells are +-1 in y-direction */
         if (is_cell_in_bounds(cell->x,cell->y + 1,cell->z)) {
            turn[ct].x = (cell->x+size.x)%size.x;
            turn[ct].y = (cell->y + 1+size.y)%size.y;
            turn[ct].z = (cell->z+size.z)%size.z;
            if (not_taken[turn[ct].x][turn[ct].y][turn[ct].z]) ct++;
         }
         if (is_cell_in_bounds(cell->x,cell->y - 1,cell->z)) {
            turn[ct].x = (cell->x+size.x)%size.x;
            turn[ct].y = (cell->y - 1+size.y)%size.y;
            turn[ct].z = (cell->z+size.z)%size.z;
            if (not_taken[turn[ct].x][turn[ct].y][turn[ct].z]) ct++;
         }
      }
      if (cell->z - cell->previous->z == 0) {
         /* 2 possible cells are +-1 in z-direction */
         if (is_cell_in_bounds(cell->x,cell->y,cell->z + 1)) {
            turn[ct].x = (cell->x+size.x)%size.x;
            turn[ct].y = (cell->y+size.y)%size.y;
            turn[ct].z = (cell->z + 1+size.z)%size.z;
            if (not_taken[turn[ct].x][turn[ct].y][turn[ct].z]) ct++;
         }
         if (is_cell_in_bounds(cell->x,cell->y,cell->z - 1)) {
            turn[ct].x = (cell->x+size.x)%size.x;
            turn[ct].y = (cell->y+size.y)%size.y;
            turn[ct].z = (cell->z - 1+size.z)%size.z;
            if (not_taken[turn[ct].x][turn[ct].y][turn[ct].z]) ct++;
         }
      }
      // fprintf(stderr,"   %d turn cells in bounds and unoccupied\n",ct);
   }

   /* if there are no available cells, return zero: we can't go anywhere */
   if (go_straight == FALSE && ct == 0) return(0);

   /* select to go straight or turn */
   if ( go_straight && (((double)rand()/(RAND_MAX+1.0) < frac) || (ct == 0)) ) {
      // fprintf(stderr,"   chose straight direction\n");
      /* fprintf(stderr,"s"); */
      new_cell->x = straight.x;      
      new_cell->y = straight.y;      
      new_cell->z = straight.z;      

   } else {
      /* choose one of the valid turn cells */
      if (ct == 1) choose = 0;
      else choose = (int)((double)ct * (double)rand()/(RAND_MAX+1.0));
      /* fprintf(stderr,"."); */
      // fprintf(stderr,"   chose cell %d\n",choose);
      new_cell->x = turn[choose].x;
      new_cell->y = turn[choose].y;
      new_cell->z = turn[choose].z;
   }

   // we found the next cell
   return(TRUE);
}


/*
 * Returns TRUE (1) if the cell indexes point to a cell within
 * the bounds defined by the user
 */
int is_cell_in_bounds(int x, int y, int z) {

   // if (x>-1 && x<size.x && y>-1 && y<size.y && z>-1 && z<size.z) return(TRUE);
   // else return(FALSE);

   if (!wrap_side.x_low)
      if (x<0 || x>size.x-1) return(FALSE);
   if (!wrap_side.y_low)
      if (y<0 || y>size.y-1) return(FALSE);
   if (!wrap_side.z_low)
      if (z<0 || z>size.z-1) return(FALSE);

   return(TRUE);
}


/*
 * Write one pipe, start at its tail and work thru to the
 * first segment, write each cell individually (no long pipes)
 */
void write_pipe_geom(CELL *end,FILE *fp,int pnum,int curve_detail) {

   int seg_ct = 0;
   cell_ptr curr = end;
   cell_ptr next = end->previous;
   cell_ptr last = NULL;
   VEC *p1,*p2,*p3,*p4,*p5,*p6,*p7,*p8;
   double n1[3],n2[3],n3[3];
   double **m1,**m2;
   // int curve_detail = -1;	// if 0, use gensurf with smooth polys
                                // if 1, use spheres at corners
                                // if -1, use polygons to make square cross-section
                                // if other, use genworm with number of segments
                                // equal to value of curve_detail, default is 5
   /*double r1;		/* major radius (of turns) */
   /*double r2;		/* minor radius (of pipe segments) */
   BOOL quit_now = FALSE;
   BOOL must_transform = FALSE;
   BOOL no_next = FALSE;


   // new coordinate transforms can only use limited output detail
   if (use_cylindrical || use_spherical) {
      curve_detail = 1;
      must_transform = TRUE;
   }

   m1 = allocate_2d_array_d(3,3);
   m2 = allocate_2d_array_d(3,3);
   p1 = (VEC*)malloc(sizeof(VEC));
   p2 = (VEC*)malloc(sizeof(VEC));
   p3 = (VEC*)malloc(sizeof(VEC));
   p4 = (VEC*)malloc(sizeof(VEC));
   p5 = (VEC*)malloc(sizeof(VEC));
   p6 = (VEC*)malloc(sizeof(VEC));
   p7 = (VEC*)malloc(sizeof(VEC));
   p8 = (VEC*)malloc(sizeof(VEC));

   fprintf(fp,"\n# starting pipe number %d\n\n",pnum);
   fflush(fp);

   /* there are 12 different permutaions of a turn, 3 of a straight; 
    * the 2 end segments need to run right, flat up against the outer
    * boundary; but, pipes that end in a dead end need to be capped
    * with a sphere, or something special. Plus, we should try to 
    * simplify long straight sections with single cylinders (later) */


   /* while end is still a cell (not NULL) */
   while (curr) {

      /* is this segment one of the ends, or a central piece? */

      if (last == NULL) {
         /* this segment is the tail, last cell created */
         /* fprintf(fp,"# first cell written\n"); */
         (void) find_center(curr,p1);
         (void) find_center2(curr,n1);
         (void) find_midpoint(curr,next,p2);
         (void) find_midpoint2(curr,next,n2);
         (void) find_matrix(next,curr,m2);

         // depending on curve_detail, draw a cylinder-like segment to the next one
         if (curve_detail == 0) {
            /* if using the gensurf feature, need gensurf cylinders! */
            if (must_transform) transform_coord(p1);
            if (must_transform) transform_coord(p2);
            fprintf(fp,"!gensurf c%d p%ds%d",pnum,pnum,seg_ct);
            (void) write_smooth_straight(fp,p1,p2,2,12);
            // draw a larger sphere to denote the ending
            fprintf(fp,"c%d sphere p%ds%d\n0 0 4 %lf %lf %lf  %lf\n",
               pnum,pnum,seg_ct,p1->x,p1->y,p1->z,r2);
         } else if (curve_detail == -1) {
            extrude_square_end(fp,pnum,seg_ct,n2,m2);
         } else {
            if (must_transform) transform_coord(p1);
            if (must_transform) transform_coord(p2);
            fprintf(fp,"c%d cylinder p%ds%dc\n0 0 7 %lf %lf %lf  %lf %lf %lf  %lf\n",
               pnum,pnum,seg_ct,p1->x,p1->y,p1->z,p2->x,p2->y,p2->z,r2);
            // draw a larger sphere to denote the ending
            fprintf(fp,"c%d sphere p%ds%d\n0 0 4 %lf %lf %lf  %lf\n",
               pnum,pnum,seg_ct,p1->x,p1->y,p1->z,r2);
            if (curve_detail == 1) fprintf(fp,"c%d sphere p%ds%ds\n0 0 4 %lf %lf %lf  %lf\n",
               pnum,pnum,seg_ct,p2->x,p2->y,p2->z,r2);
         }

         /* find if the dead end is internal or on an edge - not implemented yet! */
         /* 
         if (curr->z == 0) {
         } else if (curr->z == size.z-1) {
         } else if (curr->y == 0) {
         } else if (curr->y == size.y-1) {
         } else if (curr->x == 0) {
         } else if (curr->x == size.x-1) {
         } else {
            // last cell created is internal
         }
         */

      } else {
         /* this is one of the middle segments, or the first cell, and we need to know
          * what two cells it is between */

         /* if this is the first (created) cell, then we need to determine an
          * artificial "next" cell to help us draw the curve */
         if (next == NULL) {
            quit_now = TRUE;		// because we march through backwards
            next = (CELL *)malloc(sizeof(CELL));
            next->x = 2*curr->x - last->x;
            next->y = 2*curr->y - last->y;
            next->z = 2*curr->z - last->z;
            // if (valid_start.z_low && (curr->z == 0)) next->z = -1;
            // else if (valid_start.z_high && (curr->z == size.z-1)) next->z = size.z;
            // else if (valid_start.y_low && (curr->y == 0)) next->y = -1;
            // else if (valid_start.y_high && (curr->y == size.y-1)) next->y = size.y;
            // else if (valid_start.x_low && (curr->x == 0)) next->x = -1;
            // else if (valid_start.x_high && (curr->x == size.x-1)) next->x = size.x;
            // fprintf(stderr,"curr is %d %d %d\n",curr->x,curr->y,curr->z);
            // fprintf(stderr,"artificial next is %d %d %d\n",next->x,next->y,next->z);
         }

         /* it's either a straight segment within this cell,
          * or it's a bend (then it's tough) */
         if ((next->x == last->x && next->y == last->y) ||
             (next->x == last->x && next->z == last->z) ||
             (next->z == last->z && next->y == last->y)) {

            /* curr is a straight segment */
            // start point and direction matrix
            (void) find_midpoint(last,curr,p1);
            (void) find_midpoint2(last,curr,n1);
            (void) find_matrix(last,curr,m1);
            // end point and direction matrix
            (void) find_midpoint(curr,next,p2);
            (void) find_midpoint2(curr,next,n2);
            (void) find_matrix(curr,next,m2);

            /* fprintf(fp,"# middle cell written - straight\n"); */
            if (curve_detail == 0) {
               /* if using the gensurf feature, need gensurf cylinders! */
               fprintf(fp,"!gensurf c%d p%ds%d",pnum,pnum,seg_ct);
               (void) write_smooth_straight(fp,p1,p2,2,12);
            } else if (curve_detail == -1) {
               if (quit_now) {
                  // it's the other endcap
                  extrude_square_end(fp,pnum,seg_ct,n1,m1);
               } else {
                  // it's a central, straight section
                  extrude_square_straight(fp,pnum,seg_ct,n1,m1,n2,m2);
               }
            } else {
               if (must_transform) transform_coord(p1);
               if (must_transform) transform_coord(p2);
               fprintf(fp,"c%d cylinder p%ds%d\n0 0 7 %lf %lf %lf  %lf %lf %lf  %lf\n",
                  pnum,pnum,seg_ct,p1->x,p1->y,p1->z,p2->x,p2->y,p2->z,r2);
               if (curve_detail == 1) fprintf(fp,"c%d sphere p%ds%ds\n0 0 4 %lf %lf %lf  %lf\n",
                  pnum,pnum,seg_ct,p2->x,p2->y,p2->z,r2);
            }

         } else {

            /* curr is a curved segment */
            (void) find_midpoint(curr,last,p1);
            (void) find_midpoint2(last,curr,n1);
            (void) find_matrix(last,curr,m1);
            (void) find_center(curr,p2);
            (void) find_center2(curr,n3);
            (void) find_midpoint(curr,next,p3);
            (void) find_midpoint2(curr,next,n2);
            (void) find_matrix(curr,next,m2);
            /* fprintf(fp,"# middle cell written - curved\n"); */
            if (curve_detail == 1) {
               /* use a straight segment to denote the curve */
               if (must_transform) transform_coord(p1);
               if (must_transform) transform_coord(p2);
               if (must_transform) transform_coord(p3);
               fprintf(fp,"c%d cylinder p%ds%d0\n0 0 7 %lf %lf %lf  %lf %lf %lf  %lf\n",
                  pnum,pnum,seg_ct,p1->x,p1->y,p1->z,p2->x,p2->y,p2->z,r2);
               fprintf(fp,"c%d sphere p%ds%d1\n0 0 4 %lf %lf %lf  %lf\n",
                  pnum,pnum,seg_ct,p2->x,p2->y,p2->z,r2);
               fprintf(fp,"c%d cylinder p%ds%d2\n0 0 7 %lf %lf %lf  %lf %lf %lf  %lf\n",
                  pnum,pnum,seg_ct,p2->x,p2->y,p2->z,p3->x,p3->y,p3->z,r2);
               fprintf(fp,"c%d sphere p%ds%d3\n0 0 4 %lf %lf %lf  %lf\n",
                  pnum,pnum,seg_ct,p3->x,p3->y,p3->z,r2);
            } else if (curve_detail == 0) {
               /* use a gensurf command with smoothing for the curve */
               if (must_transform) transform_coord(p1);
               if (must_transform) transform_coord(p2);
               if (must_transform) transform_coord(p3);
               fprintf(fp,"!gensurf c%d p%ds%d",pnum,pnum,seg_ct);
               (void) write_smooth_curve(fp,p1,p2,p3,12,12);
            } else if (curve_detail == -1) {
               extrude_square_curve(fp,pnum,seg_ct,n1,m1,n2,m2);
            } else {
               /* use multiple segments/spheres for the curve */
               if (must_transform) transform_coord(p1);
               if (must_transform) transform_coord(p2);
               if (must_transform) transform_coord(p3);
               fprintf(fp,"!genworm c%d p%ds%d",pnum,pnum,seg_ct);
               (void) write_curve(fp,p1,p2,p3,curve_detail);
            }
         }
         if (quit_now) next = NULL;
      }

      fprintf(fp,"# segment %d, cell %d %d %d\n",seg_ct,curr->x,curr->y,curr->z);
      fflush(fp);
      seg_ct++;
      last = curr;
      curr = next;
      if (curr) next = next->previous;
   }
}


/*
 * Write a simple PGM description of the top view of the pipes
 *
 * Still need to write a better version for cylindrical, spherical mappings
 */
void write_data_preview(int write_pgm,int write_png,char* outfile,int pnum) {

   int i,j,k,max_val,raw_val;
   int ii,jj;
   int scaleup = 3;
   int height,width;
   double halfwidth;
   float **a;
   float new_val,already_there_val;
   VEC *pt;
   CELL *acell;
   BOOL write_binary = TRUE;
   FILE *fp;

   // initialize pt and acell
   pt = (VEC*)malloc(sizeof(VEC));
   acell = (CELL*)malloc(sizeof(CELL));

   /*
    * Make the output array
    */

   // decide which angle to view from, different if cyl, sph mappings
   if (use_cylindrical) {
      // force to view from above, xy plane
      // halfwidth = size.x*world.x+cyl_data[0];
      halfwidth = world.x*(size.x+cyl_data[0]);
      height = (int)(2.0*halfwidth);
      width = height;
      // allocate and zero the array space
      a = allocate_2d_array_f(width,height);
      for (i=0; i<width; i++) for (j=0; j<height; j++) a[i][j] = 0.0;
      max_val = size.z-1;
      // march through the not_taken array, to find appropriate color
      for (i=0; i<size.y; i++) {
         acell->y = i;
         for (j=0; j<size.x; j++) {
            acell->x = j;
            raw_val = 0;
            for (k=0; k<max_val; k++)
               if (not_taken[j][i][k] == FALSE)
                  raw_val = k;
            if (raw_val > -1) {
               // transform coordinate to draw onto a[][] array location
               acell->z = raw_val;
               (void) find_center(acell,pt);
               transform_coord(pt);
               // a[(int)(pt->y+halfwidth)][(int)(pt->x+halfwidth)] = raw_val;
               a[(int)(pt->y+halfwidth)][(int)(pt->x+halfwidth)] = (float)raw_val/(float)max_val;
            }
         }
      }

   } else if (use_spherical) {
      // force to view from side, xz plane
      halfwidth = world.x*(size.x+sph_data[0]);
      height = (int)(2.0*halfwidth);
      width = height;
      // allocate and zero the array space
      a = allocate_2d_array_f(width,height);
      for (i=0; i<width; i++) for (j=0; j<height; j++) a[i][j] = 0.0;
      max_val = 255;
      // march through the not_taken array, to find appropriate color
      for (i=0; i<size.y; i++) {
         acell->y = i;
         for (j=0; j<size.x; j++) {
            acell->x = j;
            for (k=0; k<size.z; k++) {
               if (not_taken[j][i][k] == FALSE) {
                  // transform coordinate to draw onto a[][] array location
                  acell->z = k;
                  (void) find_center(acell,pt);
                  transform_coord(pt);
                  // now, transform the y coordinate to a value!
                  new_val = (halfwidth-pt->y)/(float)height;
                  already_there_val = a[(int)(pt->z+halfwidth)][(int)(pt->x+halfwidth)];
                  if (new_val > already_there_val) {
                     a[(int)(pt->z+halfwidth)][(int)(pt->x+halfwidth)] = new_val;
                     if (new_val < 0.) 
                        fprintf(stdout,"y of %g, halfw of %g, height of %d, makes %d\n",pt->y,halfwidth,height,(int)(256*(halfwidth-pt->y)/height));
                  }
               }
            }
         }
      }

   } else if (size.x < size.y && size.x < size.z) {
      // x is the smallest dimension
      width = size.y;
      height = size.z;
      // allocate and zero the array space
      a = allocate_2d_array_f(width,height);
      for (i=0; i<width; i++) for (j=0; j<height; j++) a[i][j] = 0.0;
      max_val = size.x-1;
      for (i=0; i<width; i++) {
         for (j=0; j<height; j++) {
            raw_val = 0;
            for (k=0; k<max_val; k++) if (not_taken[k][i][j] == FALSE) raw_val = k;
            a[i][j] = (float)raw_val/(float)max_val;
         }
      }
   } else if (size.y < size.x && size.y < size.z) {
      // y is the smallest dimension
      width = size.z;
      height = size.x;
      // allocate and zero the array space
      a = allocate_2d_array_f(width,height);
      for (i=0; i<width; i++) for (j=0; j<height; j++) a[i][j] = 0.0;
      max_val = size.y-1;
      for (i=0; i<width; i++) {
         for (j=0; j<height; j++) {
            raw_val = 0;
            for (k=0; k<max_val; k++) if (not_taken[j][k][i] == FALSE) raw_val = k;
            a[i][j] = (float)raw_val/(float)max_val;
         }
      }
   } else {
      // z is the smallest dimension
      width = size.x;
      height = size.y;
      // allocate and zero the array space
      a = allocate_2d_array_f(width,height);
      for (i=0; i<width; i++) for (j=0; j<height; j++) a[i][j] = 0.0;
      max_val = size.z-1;
      for (i=0; i<width; i++) {
         for (j=0; j<height; j++) {
            raw_val = 0;
            for (k=0; k<max_val; k++) if (not_taken[i][j][k] == FALSE) raw_val = k;
            a[i][j] = (float)raw_val/(float)max_val;
         }
      }
   }

   // write the image file
   if (write_pgm) write_pgm_image(outfile,width,height,a);
   if (write_png) write_png_image(outfile,width,height,a);

   // all's well...
   return;
}


/*
 * Write a 2D array to a PGM file
 */
int write_pgm_image(char *filename,int nx,int ny,float **data) {

   int write_binary = TRUE;
   int bit_depth = 8;
   int i,j,printval;
   unsigned char bytes[2];
   FILE *fp;

   // split on bit depth
   if (write_binary) {
      // open file for writing
      fp = fopen(filename,"wb");
      if (fp==NULL) {
         fprintf(stderr,"Could not open output file %s\n",filename);
         fflush(stderr);
         exit(0);
      }
      if (bit_depth == 16) {
         fprintf(fp,"P5\n%d %d\n%d\n",nx,ny,65535);
         for (j=ny-1; j>=0; j--) {
            for (i=0; i<nx; i++) {
               printval = (int)(0.5 + 65536*data[i][j]);
               if (printval<0) printval = 0;
               else if (printval>65535) printval = 65535;
               bytes[0] = (unsigned char)(printval/256);
               bytes[1] = (unsigned char)(printval%256);
               fwrite(&bytes,sizeof(char),2,fp);
               // DOESN"T WORK !!!
            }
         }
      } else {
         fprintf(fp,"P5\n%d %d\n%d\n",nx,ny,255);
         for (j=ny-1; j>=0; j--) {
            for (i=0; i<nx; i++) {
               printval = (int)(0.5 + 256*data[i][j]);
               if (printval<0) printval = 0;
               else if (printval>255) printval = 255;
               bytes[0] = (unsigned char)(printval);
               fwrite(&bytes,sizeof(char),1,fp);
            }
         }
      }
   } else {
      // open file for writing
      fp = fopen(filename,"w");
      if (fp==NULL) {
         fprintf(stderr,"Could not open output file %s\n",filename);
         fflush(stderr);
         exit(0);
      }
      if (bit_depth == 16) {
         fprintf(fp,"P2\n%d %d\n%d\n",nx,ny,65535);
         for (j=ny-1; j>=0; j--) {
            for (i=0; i<nx; i++) {
               printval = (int)(0.5 + 65536*data[i][j]);
               if (printval<0) printval = 0;
               else if (printval>65535) printval = 65535;
               fprintf(fp,"%d\n",printval);
            }
         }
      } else {
         fprintf(fp,"P2\n%d %d\n%d\n",nx,ny,255);
         for (j=ny-1; j>=0; j--) {
            for (i=0; i<nx; i++) {
               printval = (int)(0.5 + 256*data[i][j]);
               if (printval<0) printval = 0;
               else if (printval>255) printval = 255;
               fprintf(fp,"%d\n",printval);
            }
         }
      }
   }

   // close file
   fclose(fp);

   // if all went well, return zero
   return(0);
}


/*
 * Write a 2D array to a PNG file
 */
int write_png_image(char *filename,int nx,int ny,float **data) {

   int bit_depth = 8;
   int i,j,printval;
   FILE *fp;
   png_byte **img;
   png_uint_32 k,height,width;
   png_structp png_ptr;
   png_infop info_ptr;
   png_colorp palette;
   png_voidp user_error_ptr;

   // allocate the space for the special array
   img = allocate_2d_array_pb(nx,ny,bit_depth);

   // convert the floating-point data to png_byte
   if (bit_depth == 16) {
      for (j=ny-1; j>=0; j--) {
         for (i=0; i<nx; i++) {
            printval = (int)(0.5 + 65536*data[i][j]);
            if (printval<0) printval = 0;
            else if (printval>65535) printval = 65535;
            img[ny-1-j][2*i] = (png_byte)(printval/256);
            img[ny-1-j][2*i+1] = (png_byte)(printval%256);
         }
      }
   } else {
      bit_depth = 8;
      for (j=ny-1; j>=0; j--) {
         for (i=0; i<nx; i++) {
            printval = (int)(0.5 + 256*data[i][j]);
            if (printval<0) printval = 0;
            else if (printval>255) printval = 255;
            img[ny-1-j][i] = (png_byte)printval;
         }
      }
   }

   // write the file
   height=ny;
   width=nx;

   // open file for writing
   fp = fopen(filename,"wb");
   if (fp==NULL) {
      fprintf(stderr,"Could not open output file %s\n",filename);
      fflush(stderr);
      exit(0);
   }

   /* Create and initialize the png_struct with the desired error handler
    * functions.  If you want to use the default stderr and longjump method,
    * you can supply NULL for the last three parameters.  We also check that
    * the library version is compatible with the one used at compile time,
    * in case we are using dynamically linked libraries.  REQUIRED.
    */
   png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING,
      NULL, NULL, NULL);
      // user_error_ptr, user_error_fn, user_warning_fn);

   if (png_ptr == NULL) {
      fclose(fp);
      fprintf(stderr,"Could not create png struct\n");
      fflush(stderr);
      exit(0);
      return (-1);
   }

   /* Allocate/initialize the image information data.  REQUIRED */
   info_ptr = png_create_info_struct(png_ptr);
   if (info_ptr == NULL) {
      fclose(fp);
      png_destroy_write_struct(&png_ptr,(png_infopp)NULL);
      return (-1);
   }

   /* Set error handling.  REQUIRED if you aren't supplying your own
    * error handling functions in the png_create_write_struct() call.
    */
   if (setjmp(png_jmpbuf(png_ptr))) {
      /* If we get here, we had a problem reading the file */
      fclose(fp);
      png_destroy_write_struct(&png_ptr, &info_ptr);
      return (-1);
   }

   /* set up the output control if you are using standard C streams */
   png_init_io(png_ptr, fp);

   /* Set the image information here.  Width and height are up to 2^31,
    * bit_depth is one of 1, 2, 4, 8, or 16, but valid values also depend on
    * the color_type selected. color_type is one of PNG_COLOR_TYPE_GRAY,
    * PNG_COLOR_TYPE_GRAY_ALPHA, PNG_COLOR_TYPE_PALETTE, PNG_COLOR_TYPE_RGB,
    * or PNG_COLOR_TYPE_RGB_ALPHA.  interlace is either PNG_INTERLACE_NONE or
    * PNG_INTERLACE_ADAM7, and the compression_type and filter_type MUST
    * currently be PNG_COMPRESSION_TYPE_BASE and PNG_FILTER_TYPE_BASE. REQUIRED
    */
   png_set_IHDR(png_ptr, info_ptr, height, width, bit_depth,
      PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE,
      PNG_FILTER_TYPE_BASE);

   /* Optional gamma chunk is strongly suggested if you have any guess
    * as to the correct gamma of the image. */
   png_set_gAMA(png_ptr, info_ptr, 2.2);

   /* Write the file header information.  REQUIRED */
   png_write_info(png_ptr, info_ptr);

   /* One of the following output methods is REQUIRED */
   // png_write_image(png_ptr, row_pointers);
   png_write_image(png_ptr, img);

   /* It is REQUIRED to call this to finish writing the rest of the file */
   png_write_end(png_ptr, info_ptr);

   /* clean up after the write, and free any memory allocated */
   png_destroy_write_struct(&png_ptr, &info_ptr);

   // close file
   fclose(fp);

   // if all went well, return zero
   return(0);
}


/*
 * This function writes basic usage information to stderr,
 * and then quits. Too bad.
 */
int Usage(char progname[80],int status) {

   /* Usage for pipes */
   static char **cpp, *help_message[] =
   {
       "where options include:",
       "  -help          writes this help information",
       "  -rad [name]    write triangles to Radiance file, turns off -obj",
       "  -obj [name]    write triangles to wavefront file, turns off -rad",
       "  -pgm [name]    write PGM preview image, turns off -png",
       "  -png [name]    write PNG preview image, turns off -pgm",
       "  -s x y z       integer size of available space, default is 20 20 5",
       "  -n val         number of pipes, def= 1",
       "  -l min max     length of pipes, minimum and maximum, def= 10 10000",
       "  -f val         probability that a segment will go straight if it can;",
       "                 real value, def= 0.7, range 0-1",
       "  -r val         radius of extrusion; range 0:0.5, default 0.2",
       "  -sq            use square cross-section, turns off -gs and -worm",
       "  -gs            use gensurf extrusions, turns off -sq and -worm",
       "  -worm num      use genworm extrusions with num segments",
       "  -[wx|wy|wz]    volume is periodic (wraps around) in [x|y|z] dimension",
       "  -s[x|y|z][l|h] start pipes on x,y,z low,high face; ex: -syh -szl",
       "  -cyl orad nc   use cylindrical projection, x is radial direction, y is theta,",
       "                 z is still z; orad is the radial offset, in cells; nc is the",
       "                 number of arc copies required to fill the whole circle;",
       "                 example: -cyl 10 4 means that the closest any pipe segments",
       "                 will get to the center is 10 units, and the resulting pipe",
       "                 volume will cover a quadrant of a circle (90 degrees)",
       "  -sph orad nc phi1 phi2",
       "                 use a spherical projection, again, x is radial direction,",
       "                 y is theta (longitude), but z is phi (latitude); phi1 is the",
       "                 latitude in degrees from the equator of the lowest z pipe,",
       "                 phi2 is the same for the highest z pipe;",
       "                 example: -sph 2 4 0 89 will fill one octant of a sphere",
       " ",
       "Options may be abbreviated to an unambiguous length.",
       "Filenames may be omitted, output filenames will then be out.xxx",
       NULL
   };

   fprintf(stdout, "usage:\n  %s [-options ...]\n\n", progname);
   for (cpp = help_message; *cpp; cpp++)
      fprintf(stdout, "%s\n", *cpp);
   exit(status);
   return(0);
}


/*
 * Free the memory associated with the last pipe
 */
void free_mem(CELL *end) {

   CELL *prev;

   /* while end is still a cell */
   while (end) {
      prev = end->previous;
      free(end);
      end = prev;
   }
}


/*
 * Transform the coordinate to either cylindrical or spherical coordinates!
 */
void transform_coord(VEC *point) {

   double rad,theta,phi;

   if (use_cylindrical) {
      rad = point->x + cyl_data[0];
      theta = 2.0*M_PI/cyl_data[1]*point->y/size.y/world.y;
      point->x = rad * cos(theta);
      point->y = rad * sin(theta);
   } else if (use_spherical) {
      rad = point->x + sph_data[0];
      theta = 2.0*M_PI/sph_data[1]*point->y/size.y/world.y;
      phi = M_PI/180.0 * (point->z/size.z/world.z*(sph_data[3]-sph_data[2])+sph_data[2]);
      point->x = rad * cos(theta) * cos(phi);
      point->y = rad * sin(theta) * cos(phi);
      point->z = rad * sin(phi);
   }

}


/*
 * Return the midpoint of the common face of two cells
 * If the cells are separated across a wrapped face, return the value
 * on the "start" side of the face.
 */
void find_midpoint(CELL *start, CELL *end, VEC *point) {

   // double csd2 = cell_size/2.0;
   if (wrap_side.x_low && (abs(start->x-end->x) > 1)) {
      if (start->x == 0) point->x = 0.0;
      else point->x = (double)size.x * world.x;
      // fprintf(stdout,"wrapping in x, start %d, end %d, point %g\n",start->x,end->x,point->x);
   } else {
      point->x = 0.5 * world.x * (1 + (int)start->x + (int)end->x);
   }
   if (wrap_side.y_low && (abs(start->y-end->y) > 1)) {
      if (start->y == 0) point->y = 0.0;
      else point->y = (double)size.y * world.y;
   } else {
      point->y = 0.5 * world.y * (1 + (int)start->y + (int)end->y);
   }
   if (wrap_side.z_low && (abs(start->z-end->z) > 1)) {
      if (start->z == 0) point->z = 0.0;
      else point->z = (double)size.z * world.z;
   } else {
      point->z = 0.5 * world.z * (1 + (int)start->z + (int)end->z);
   }
}


/*
 * Return the midpoint of the common face of two cells
 * If the cells are separated across a wrapped face, return the value
 * on the "start" side of the face.
 */
void find_midpoint2(CELL *start, CELL *end, double *point) {

   // double csd2 = cell_size/2.0;
   if (wrap_side.x_low && (abs(start->x-end->x) > 1)) {
      if (start->x == 0) point[0] = 0.0;
      else point[0] = (double)size.x * world.x;
      // fprintf(stdout,"wrapping in x, start %d, end %d, point %g\n",start->x,end->x,point[0]);
   } else {
      point[0] = 0.5 * world.x * (1 + (int)start->x + (int)end->x);
   }
   if (wrap_side.y_low && (abs(start->y-end->y) > 1)) {
      if (start->y == 0) point[1] = 0.0;
      else point[1] = (double)size.y * world.y;
   } else {
      point[1] = 0.5 * world.y * (1 + (int)start->y + (int)end->y);
   }
   if (wrap_side.z_low && (abs(start->z-end->z) > 1)) {
      if (start->z == 0) point[2] = 0.0;
      else point[2] = (double)size.z * world.z;
   } else {
      point[2] = 0.5 * world.z * (1 + (int)start->z + (int)end->z);
   }
}


/*
 * Return the cosine matrix of the point at the center of the common face
 * of two cells
 * If the cells are separated across a wrapped face, return the value
 * on the "start" side of the face.
 * The local z-direction (m[2]) points along the direction of extrusion
 * The local x-direction (m[0]) points along a dir perp to z
 * The local y-direction (m[1]) points along a dir perp to both x and z
 */
void find_matrix(CELL *start, CELL *end, double **m) {

   int i,j;

   for (i=0;i<3;i++) for (j=0;j<3;j++) m[i][j] = 0.;

   if (abs(start->x-end->x) == 0) {
      if (abs(start->y-end->y) == 0) {
         // segment runs along z-axis
         m[2][2] = (double)(end->z - start->z);
         if (fabs(m[2][2]) > 1.1) m[2][2] = -1.*m[2][2];
         m[1][1] = -1.*m[2][2];
         vec_cross(m[1],m[2],m[0]);
      } else {
         // segment runs along y-axis
         m[2][1] = (double)(end->y - start->y);
         if (fabs(m[2][1]) > 1.1) m[2][1] = -1.*m[2][1];
         m[1][0] = -1.*m[2][1];
         vec_cross(m[1],m[2],m[0]);
      }
   } else {
         // segment runs along x-axis
         m[2][0] = (double)(end->x - start->x);
         if (fabs(m[2][0]) > 1.1) m[2][0] = -1.*m[2][0];
         m[1][2] = -1.*m[2][0];
         vec_cross(m[1],m[2],m[0]);
   }

   // fprintf(stderr,"matrix is\n");
   // fprintf(stderr,"  %g %g %g\n",m[0][0],m[0][1],m[0][2]);
   // fprintf(stderr,"  %g %g %g\n",m[1][0],m[1][1],m[1][2]);
   // fprintf(stderr,"  %g %g %g\n",m[2][0],m[2][1],m[2][2]);

   return;
}

void vec_cross(double *x,double *y,double *z) {
   z[0] = x[1]*y[2]-x[2]*y[1];
   z[1] = x[2]*y[0]-x[0]*y[2];
   z[2] = x[0]*y[1]-x[1]*y[0];
   return;
}


/*
 * Return the center of the cell, cartesian only
 */
void find_center(CELL *cell, VEC *point) {
   point->x = 0.5*world.x * (1 + 2*(int)cell->x);
   point->y = 0.5*world.y * (1 + 2*(int)cell->y);
   point->z = 0.5*world.z * (1 + 2*(int)cell->z);
}


/*
 * Return the center of the cell, cartesian only
 */
void find_center2(CELL *cell, double *point) {
   point[0] = 0.5*world.x * (1 + 2*(int)cell->x);
   point[1] = 0.5*world.y * (1 + 2*(int)cell->y);
   point[2] = 0.5*world.z * (1 + 2*(int)cell->z);
}


/*
 * Draw a curved section of pipe in the given cell using genworm
 * using st(art), ce(nter), and en(d) points
 */
void write_curve(FILE *fp, VEC *st, VEC *ce, VEC *en, int segs) {

   double po2 = M_PI/2.0;

   /* put them together */
   fprintf(fp," '%lf + %lf*(1-cos(t*%lf)) + %lf*sin(t*%lf)'",st->x,en->x-ce->x,po2,ce->x-st->x,po2);
   fprintf(fp," '%lf + %lf*(1-cos(t*%lf)) + %lf*sin(t*%lf)'",st->y,en->y-ce->y,po2,ce->y-st->y,po2);
   fprintf(fp," '%lf + %lf*(1-cos(t*%lf)) + %lf*sin(t*%lf)'",st->z,en->z-ce->z,po2,ce->z-st->z,po2);
   fprintf(fp," '%lf' %d\n",r2,segs);	/* r2 is global */
}


/*
 * Draw a curved section of pipe in the given cell using gensurf formulae
 * using st(art), ce(nter), and en(d) points
 */
void write_smooth_curve(FILE *fp, VEC *st, VEC *ce, VEC *en, int sres, int tres) {

   double po2 = M_PI/2.0;
   double tp = M_PI*2.0;
   double length;
   VEC *e2;

   e2 = (VEC *)malloc(sizeof(VEC));

   /* create vector normal to the plane of the curve */
   e2->x = (ce->y-st->y)*(en->z-ce->z) - (ce->z-st->z)*(en->y-ce->y);
   e2->y = (ce->z-st->z)*(en->x-ce->x) - (ce->x-st->x)*(en->z-ce->z);
   e2->z = (ce->x-st->x)*(en->y-ce->y) - (ce->y-st->y)*(en->x-ce->x);
   length = sqrt(pow(e2->x,2) + pow(e2->y,2) + pow(e2->z,2));
   e2->x = r2 * e2->x / length;
   e2->y = r2 * e2->y / length;
   e2->z = r2 * e2->z / length;

   /* s is the parametric curve, t is the circumference */
   fprintf(fp," '%lf+%lf*(1-cos(s*%lf))+%lf*sin(s*%lf) +",st->x,en->x-ce->x,po2,ce->x-st->x,po2);
   fprintf(fp," (%lf*cos(s*%lf)+%lf*sin(s*%lf))*cos(t*%lf)*%lf +",ce->x-en->x,po2,ce->x-st->x,po2,tp,r2/r1);
   fprintf(fp," %lf*sin(t*%lf)'",e2->x,tp);

   fprintf(fp," '%lf+%lf*(1-cos(s*%lf))+%lf*sin(s*%lf) +",st->y,en->y-ce->y,po2,ce->y-st->y,po2);
   fprintf(fp," (%lf*cos(s*%lf)+%lf*sin(s*%lf))*cos(t*%lf)*%lf +",ce->y-en->y,po2,ce->y-st->y,po2,tp,r2/r1);
   fprintf(fp," %lf*sin(t*%lf)'",e2->y,tp);

   fprintf(fp," '%lf+%lf*(1-cos(s*%lf))+%lf*sin(s*%lf) +",st->z,en->z-ce->z,po2,ce->z-st->z,po2);
   fprintf(fp," (%lf*cos(s*%lf)+%lf*sin(s*%lf))*cos(t*%lf)*%lf +",ce->z-en->z,po2,ce->z-st->z,po2,tp,r2/r1);
   fprintf(fp," %lf*sin(t*%lf)'",e2->z,tp);

   fprintf(fp," %d %d -s\n",sres,tres);

}


/*
 * Draw a straight section of pipe in the given cell using gensurf formulae
 * this is neccessary because gensurf actually defines small polygons for
 * the curved sections, and using a cylinder primitive for these straight
 * sections reveals cracks in the seam between straight and curved sections.
 * Hence, we must make straight sections out of polygons, too. 
 */
void write_smooth_straight(FILE *fp, VEC *st, VEC *en, int sres, int tres) {

   double tp = M_PI*2.0;
   VEC *e1;
   VEC *e2;

   e1 = (VEC *)malloc(sizeof(VEC));
   e2 = (VEC *)malloc(sizeof(VEC));

   /* create vector normal to the plane of the curve */
   if (fabs(en->x-st->x) > 0.001) {
      e1->x = 0.0;
      e1->y = r2;
      e1->z = 0.0;
      e2->x = 0.0;
      e2->y = 0.0;
      e2->z = r2;
   } else if (fabs(en->y-st->y) > 0.001) {
      e1->x = r2;
      e1->y = 0.0;
      e1->z = 0.0;
      e2->x = 0.0;
      e2->y = 0.0;
      e2->z = r2;
   } else {
      e1->x = 0.0;
      e1->y = r2;
      e1->z = 0.0;
      e2->x = r2;
      e2->y = 0.0;
      e2->z = 0.0;
   }

   /* s is the parametric curve, t is the circumference */
   fprintf(fp," '%lf+%lf*s +",st->x,en->x-st->x);
   fprintf(fp," %lf*cos(t*%lf) + %lf*sin(t*%lf)'",e1->x,tp,e2->x,tp);

   fprintf(fp," '%lf+%lf*s +",st->y,en->y-st->y);
   fprintf(fp," %lf*cos(t*%lf) + %lf*sin(t*%lf)'",e1->y,tp,e2->y,tp);

   fprintf(fp," '%lf+%lf*s +",st->z,en->z-st->z);
   fprintf(fp," %lf*cos(t*%lf) + %lf*sin(t*%lf)'",e1->z,tp,e2->z,tp);

   fprintf(fp," %d %d -s\n",sres,tres);

}


/*
 * write polygons for a staright segment of square-section
 */
void extrude_square_straight(FILE *fp,int pnum,int seg_ct,double *p1,double **m1,double *p2,double **m2) {
   int i,j;
   double n[8][3];

   // define the 8 endnodes
   for (i=0;i<3;i++) n[0][i] = p1[i] - r2*m1[0][i] - r2*m1[1][i];
   for (i=0;i<3;i++) n[1][i] = p1[i] - r2*m1[0][i] + r2*m1[1][i];
   for (i=0;i<3;i++) n[2][i] = p1[i] + r2*m1[0][i] + r2*m1[1][i];
   for (i=0;i<3;i++) n[3][i] = p1[i] + r2*m1[0][i] - r2*m1[1][i];
   for (i=0;i<3;i++) n[4][i] = p2[i] - r2*m2[0][i] - r2*m2[1][i];
   for (i=0;i<3;i++) n[5][i] = p2[i] - r2*m2[0][i] + r2*m2[1][i];
   for (i=0;i<3;i++) n[6][i] = p2[i] + r2*m2[0][i] + r2*m2[1][i];
   for (i=0;i<3;i++) n[7][i] = p2[i] + r2*m2[0][i] - r2*m2[1][i];

   // fprintf(stderr,"8 endnodes are\n");
   // fprintf(stderr,"  %g %g %g\n",n[0][0],n[0][1],n[0][2]);
   // fprintf(stderr,"  %g %g %g\n",n[1][0],n[1][1],n[1][2]);
   // fprintf(stderr,"  %g %g %g\n",n[2][0],n[2][1],n[2][2]);
   // fprintf(stderr,"  %g %g %g\n",n[3][0],n[3][1],n[3][2]);
   // fprintf(stderr,"  %g %g %g\n",n[4][0],n[4][1],n[4][2]);
   // fprintf(stderr,"  %g %g %g\n",n[5][0],n[5][1],n[5][2]);
   // fprintf(stderr,"  %g %g %g\n",n[6][0],n[6][1],n[6][2]);
   // fprintf(stderr,"  %g %g %g\n",n[7][0],n[7][1],n[7][2]);
   // exit(0);

   // write the requisite polygons
   write_polygon_4(fp,pnum,seg_ct,n[3],n[2],n[6],n[7]);
   write_polygon_4(fp,pnum,seg_ct,n[2],n[1],n[5],n[6]);
   write_polygon_4(fp,pnum,seg_ct,n[1],n[0],n[4],n[5]);
   write_polygon_4(fp,pnum,seg_ct,n[0],n[3],n[7],n[4]);

   return;
}


/*
 * write polygons for a staright segment of square-section
 */
void extrude_square_end(FILE *fp,int pnum,int seg_ct,double *p1,double **m1) {

   int i;
   double n[8][3];

   // define the 8 endnodes
   for (i=0;i<3;i++) n[0][i] = p1[i] - r2*m1[0][i] - r2*m1[1][i];
   for (i=0;i<3;i++) n[1][i] = p1[i] - r2*m1[0][i] + r2*m1[1][i];
   for (i=0;i<3;i++) n[2][i] = p1[i] + r2*m1[0][i] + r2*m1[1][i];
   for (i=0;i<3;i++) n[3][i] = p1[i] + r2*m1[0][i] - r2*m1[1][i];
   for (i=0;i<3;i++) n[4][i] = p1[i] - r2*m1[0][i] - r2*m1[1][i] + (0.5+r2)*m1[2][i];
   for (i=0;i<3;i++) n[5][i] = p1[i] - r2*m1[0][i] + r2*m1[1][i] + (0.5+r2)*m1[2][i];
   for (i=0;i<3;i++) n[6][i] = p1[i] + r2*m1[0][i] + r2*m1[1][i] + (0.5+r2)*m1[2][i];
   for (i=0;i<3;i++) n[7][i] = p1[i] + r2*m1[0][i] - r2*m1[1][i] + (0.5+r2)*m1[2][i];

   // write the requisite polygons
   write_polygon_4(fp,pnum,seg_ct,n[3],n[2],n[6],n[7]);
   write_polygon_4(fp,pnum,seg_ct,n[2],n[1],n[5],n[6]);
   write_polygon_4(fp,pnum,seg_ct,n[1],n[0],n[4],n[5]);
   write_polygon_4(fp,pnum,seg_ct,n[0],n[3],n[7],n[4]);

   // write the endcap
   write_polygon_4(fp,pnum,seg_ct,n[7],n[6],n[5],n[4]);

   return;
}


/*
 * write polygons for a "curved" segment of square-section
 */
void extrude_square_curve(FILE *fp,int pnum,int seg_ct,double *p1,double **m1,double *p2,double **m2) {

   int i;
   double p3[3],p4[3];
   double n[8][3];
   double m3[3][3];

   // write the intro
   for (i=0;i<3;i++) p3[i] = p1[i] + (0.5-r2)*m1[2][i];
   extrude_square_straight(fp,pnum,seg_ct,p1,m1,p3,m1);

   // write the outtro
   for (i=0;i<3;i++) p4[i] = p2[i] - (0.5-r2)*m2[2][i];
   extrude_square_straight(fp,pnum,seg_ct,p4,m2,p2,m2);

   // write the requisite polygons for the central cube
   // opposite the entry
   for (i=0;i<3;i++) p3[i] = p1[i] + (0.5+r2)*m1[2][i];
   for (i=0;i<3;i++) n[0][i] = p3[i] - r2*m1[0][i] - r2*m1[1][i];
   for (i=0;i<3;i++) n[1][i] = p3[i] - r2*m1[0][i] + r2*m1[1][i];
   for (i=0;i<3;i++) n[2][i] = p3[i] + r2*m1[0][i] + r2*m1[1][i];
   for (i=0;i<3;i++) n[3][i] = p3[i] + r2*m1[0][i] - r2*m1[1][i];
   write_polygon_4(fp,pnum,seg_ct,n[3],n[2],n[1],n[0]);

   // opposite the exit
   for (i=0;i<3;i++) p3[i] = p2[i] - (0.5+r2)*m2[2][i];
   for (i=0;i<3;i++) n[0][i] = p3[i] - r2*m2[0][i] - r2*m2[1][i];
   for (i=0;i<3;i++) n[1][i] = p3[i] - r2*m2[0][i] + r2*m2[1][i];
   for (i=0;i<3;i++) n[2][i] = p3[i] + r2*m2[0][i] + r2*m2[1][i];
   for (i=0;i<3;i++) n[3][i] = p3[i] + r2*m2[0][i] - r2*m2[1][i];
   write_polygon_4(fp,pnum,seg_ct,n[0],n[1],n[2],n[3]);

   /*
   fprintf(stderr,"matrix m1 is\n");
   fprintf(stderr,"  %g %g %g\n",m1[0][0],m1[0][1],m1[0][2]);
   fprintf(stderr,"  %g %g %g\n",m1[1][0],m1[1][1],m1[1][2]);
   fprintf(stderr,"  %g %g %g\n",m1[2][0],m1[2][1],m1[2][2]);

   fprintf(stderr,"matrix m2 is\n");
   fprintf(stderr,"  %g %g %g\n",m2[0][0],m2[0][1],m2[0][2]);
   fprintf(stderr,"  %g %g %g\n",m2[1][0],m2[1][1],m2[1][2]);
   fprintf(stderr,"  %g %g %g\n",m2[2][0],m2[2][1],m2[2][2]);
   */

   // define the vector perpendicular to the inlet and outlet
   // find the matrix for this plane
   vec_cross(m1[2],m2[2],m3[2]);
   // fprintf(stderr,"matrix m3 is\n");
   // fprintf(stderr,"  %g %g %g\n",m3[2][0],m3[2][1],m3[2][2]);
   for (i=0;i<3;i++) p4[i] = p1[i] + 0.5*m1[2][i] + r2*m3[2][i];
   for (i=0;i<3;i++) n[0][i] = p4[i] - r2*m1[2][i] - r2*m2[2][i];
   for (i=0;i<3;i++) n[1][i] = p4[i] - r2*m1[2][i] + r2*m2[2][i];
   for (i=0;i<3;i++) n[2][i] = p4[i] + r2*m1[2][i] + r2*m2[2][i];
   for (i=0;i<3;i++) n[3][i] = p4[i] + r2*m1[2][i] - r2*m2[2][i];
   write_polygon_4(fp,pnum,seg_ct,n[3],n[2],n[1],n[0]);

   // define the vector perpendicular to the inlet and outlet
   // find the matrix for this plane
   vec_cross(m2[2],m1[2],m3[2]);
   // fprintf(stderr,"matrix m3 is\n");
   // fprintf(stderr,"  %g %g %g\n",m3[2][0],m3[2][1],m3[2][2]);
   for (i=0;i<3;i++) p4[i] = p1[i] + 0.5*m1[2][i] + r2*m3[2][i];
   for (i=0;i<3;i++) n[0][i] = p4[i] - r2*m1[2][i] - r2*m2[2][i];
   for (i=0;i<3;i++) n[1][i] = p4[i] - r2*m1[2][i] + r2*m2[2][i];
   for (i=0;i<3;i++) n[2][i] = p4[i] + r2*m1[2][i] + r2*m2[2][i];
   for (i=0;i<3;i++) n[3][i] = p4[i] + r2*m1[2][i] - r2*m2[2][i];
   write_polygon_4(fp,pnum,seg_ct,n[0],n[1],n[2],n[3]);

   // fprintf(stderr,"\n");
   return;
}


/*
 * write a 4-node polygon
 */
void write_polygon_4(FILE *fp,int pnum,int seg_ct,double *p1,double *p2,double *p3,double *p4) {

   // if radiance polygon
   if (write_rad) {
   fprintf(fp,"c%d polygon p%ds%d\n0 0 12\n",pnum,pnum,seg_ct);
   fprintf(fp,"%lf %lf %lf\n",p1[0],p1[1],p1[2]);
   fprintf(fp,"%lf %lf %lf\n",p2[0],p2[1],p2[2]);
   fprintf(fp,"%lf %lf %lf\n",p3[0],p3[1],p3[2]);
   fprintf(fp,"%lf %lf %lf\n",p4[0],p4[1],p4[2]);
   }

   // if obj-triangles
   else if (write_obj3) {
   fprintf(fp,"v %lf %lf %lf\n",p1[0],p1[1],p1[2]);
   fprintf(fp,"v %lf %lf %lf\n",p2[0],p2[1],p2[2]);
   fprintf(fp,"v %lf %lf %lf\n",p3[0],p3[1],p3[2]);
   fprintf(fp,"v %lf %lf %lf\n",p4[0],p4[1],p4[2]);
   fprintf(fp,"f %d %d %d\n",num_v+1,num_v+2,num_v+3);
   fprintf(fp,"f %d %d %d\n",num_v+1,num_v+3,num_v+4);
   num_v += 4;
   }

   return;
}


/*
 * allocate memory for a two-dimensional array of ints
 */
int** allocate_2d_array_i(int nx, int ny) {

   int i,j;
   int **array = (int **)malloc(nx * sizeof(int *));

   array[0] = (int *)malloc(nx * ny * sizeof(int));

   for (i=1; i<nx; i++)
      array[i] = array[0] + i * ny;

   return(array);
}


/*
 * allocate memory for a two-dimensional array of floats
 */
float** allocate_2d_array_f(int nx, int ny) {

   int i,j;
   float **array = (float **)malloc(nx * sizeof(float *));

   array[0] = (float *)malloc(nx * ny * sizeof(float));

   for (i=1; i<nx; i++)
      array[i] = array[0] + i * ny;

   return(array);
}


/*
 * allocate memory for a two-dimensional array of doubles
 */
double** allocate_2d_array_d(int nx, int ny) {

   int i,j;
   double **array = (double **)malloc(nx * sizeof(double *));

   array[0] = (double *)malloc(nx * ny * sizeof(double));

   for (i=1; i<nx; i++)
      array[i] = array[0] + i * ny;

   return(array);
}


/*
 * allocate memory for a two-dimensional array of png_byte
 */
png_byte** allocate_2d_array_pb(int nx, int ny, int depth) {

   int i,j,bytesperpixel;
   png_byte **array;

   if (depth <= 8) bytesperpixel = 1;
   else bytesperpixel = 2;
   array = (png_byte **)malloc(ny * sizeof(png_byte *));
   array[0] = (png_byte *)malloc(bytesperpixel * nx * ny * sizeof(png_byte));

   for (i=1; i<ny; i++)
      array[i] = array[0] + i * bytesperpixel * nx;

   return(array);
}


/*
 * allocate memory for a three-dimensional array of BOOLs
 */
BOOL*** allocate_3d_array_B(int nx, int ny, int nz) {

   int i,j;
   BOOL ***array = (BOOL ***)malloc(nx * sizeof(BOOL **));

   array[0] = (BOOL **)malloc(nx * ny * sizeof(BOOL *));
   array[0][0] = (BOOL *)malloc(nx * ny * nz * sizeof(BOOL));

   for (i=1; i<nx; i++)
      array[i] = array[0] + i * ny;

   for (i=0; i<nx; i++) {
      if (i!=0)
         array[i][0] = array[0][0] + i * ny * nz;
      for (j=1; j<ny; j++)
         array[i][j] = array[i][0] + j * nz;
   }

   return(array);
}