#define XMIN   0.0
#define XMAX  20.0
#define YMIN   0.0
#define YMAX  20.0
#define ZMIN -15.0
#define ZMAX   0.0
#define STARTBOXSIZE 5.00
#define PERIOD 10.0

#define OCTNODEPERBLOCK 4092

/* The next line is used to outsmart some very stupid compilers.  If your    */
/*   compiler is smarter, feel free to replace the "int" with "void".        */
/*   Not that it matters.                                                    */

#define VOID int

#include <stdio.h>

#ifndef _STDLIB_H_
extern void *malloc();
extern void free();
extern void exit();
#endif

/* Labels that signify whether a record consists primarily of pointers or of */
/*   floating-point words.  Used to make decisions about data alignment.     */

enum wordtype {POINTER, FLOATINGPOINT};

enum octnodetype {LEAF, PARENT};

struct octnode {
  struct octnode *neighbor[4][2];
  enum octnodetype mytype;
  int treedepth;
  int numerator[3];
  int denominator[3];
};

/* A type used to allocate memory.  firstblock is the first block of items.  */
/*   nowblock is the block from which items are currently being allocated.   */
/*   nextitem points to the next slab of free memory for an item.            */
/*   deaditemstack is the head of a linked list (stack) of deallocated items */
/*   that can be recycled.  unallocateditems is the number of items that     */
/*   remain to be allocated from nowblock.                                   */
/*                                                                           */
/* Traversal is the process of walking through the entire list of items, and */
/*   is separate from allocation.  Note that a traversal will visit items on */
/*   the "deaditemstack" stack as well as live items.  pathblock points to   */
/*   the block currently being traversed.  pathitem points to the next item  */
/*   to be traversed.  pathitemsleft is the number of items that remain to   */
/*   be traversed in pathblock.                                              */
/*                                                                           */
/* itemwordtype is set to POINTER or FLOATINGPOINT, and is used to suggest   */
/*   what sort of word the record is primarily made up of.  alignwithnull is */
/*   a 0/1 value that determines whether new records should be aligned with  */
/*   NULL, modulo the record length.  itembytes and itemwords are the length */
/*   of a record in bytes (after rounding up) and words.  itemsperblock is   */
/*   the number of items allocated at once in a single block.  items is the  */
/*   number of currently allocated items.  maxitems is the maximum number of */
/*   items that have been allocated at once; it is the current number of     */
/*   items plus the number of records kept on deaditemstack.                 */

struct memorypool {
  VOID **firstblock, **nowblock;
  VOID *nextitem;
  VOID *deaditemstack;
  VOID **pathblock;
  VOID *pathitem;
  enum wordtype itemwordtype;
  int alignwithnull;
  int itembytes, itemwords;
  int itemsperblock;
  int items, maxitems;
  int unallocateditems;
  int pathitemsleft;
};

void initlookup_h();
/*
{
  return;
}
*/

double lookup_h();
/*
double x;
double y;
double z;
{
  return 0.005 * (x * x + y * y + z * z) + 0.5;
}
*/

/********* Memory management routines begin here                     *********/
/**                                                                         **/
/**                                                                         **/

/* A useful forward declaration.                                             */

void poolrestart();

/*****************************************************************************/
/*                                                                           */
/*  poolinit()   Initialize a pool of memory for allocation of items.        */
/*                                                                           */
/*  This routine initializes the machinery for allocating items.  A `pool'   */
/*  is created whose records have size at least `bytecount'.  Items will be  */
/*  allocated in `itemcount'-item blocks.  Each item is assumed to be a      */
/*  collection of words, and either pointers or floating-point values are    */
/*  assumed to be the "primary" word type.  (The "primary" word type is used */
/*  to determine alignment of items.)  If the align flag is set, then items  */
/*  will be aligned with NULL in memory.                                     */
/*                                                                           */
/*  The memory allocation routines are the kludgiest part of the whole       */
/*  program because when the align flag is set, I force every item to be     */
/*  aligned with NULL modulo `itemwords'.  This is to make it possible to    */
/*  store both a pointer to a triangle or shell edge, and an orientation as  */
/*  well in a single pointer.  (Such pointers will point *into* a triangle   */
/*  or shell edge, and the alignment property is used to determine where the */
/*  item starts.)                                                            */
/*                                                                           */
/*  This may cause all sorts of problems for a C compiler that represents    */
/*  NULL with an address larger than the addresses being allocated from the  */
/*  heap.                                                                    */
/*                                                                           */
/*  "There are terrible temptations that it required strength, strength and  */
/*  courage, to yield to."  - Oscar Wilde                                    */
/*                                                                           */
/*  Don't change this routine unless you understand it.                      */
/*                                                                           */
/*****************************************************************************/

void poolinit(pool, bytecount, itemcount, wtype, align)
struct memorypool *pool;
int bytecount;
int itemcount;
enum wordtype wtype;
int align;
{
  int wordsize;
  int pointers;

  /* Initialize values in the pool. */
  pool->itemwordtype = wtype;
  pool->alignwithnull = align;
  wordsize = (pool->itemwordtype == POINTER) ? sizeof(VOID *) : sizeof(double);
  /* Make sure there's room for a pointer for the dead item stack. */
  if (bytecount < sizeof(VOID *)) {
    bytecount = sizeof(VOID *);
  }
  if (sizeof(VOID *) > wordsize) {
    /* Round up the number of bytes to an integral number of pointers. */
    pointers = (bytecount + sizeof(VOID *) - 1) / sizeof(VOID *);
    pool->itemwords = pointers * (sizeof(VOID *) / wordsize);
  } else {
    /* Round up the number of bytes to an integral number of words. */
    pool->itemwords = (bytecount + wordsize - 1) / wordsize;
  }
  pool->itembytes = pool->itemwords * wordsize;
  pool->itemsperblock = itemcount;

  /* Allocate a block of items.  Space for `itemsperblock' items and one    */
  /*   pointer (to point to the next block) are allocated.  If alignment is */
  /*   necessary, an extra item is allocated to facilitate alignment.       */
  if (pool->alignwithnull) {
    pool->firstblock = (VOID **) malloc((pool->itemsperblock + 1)
                                        * pool->itembytes
                                        + (sizeof(VOID *) >= sizeof(double)
                                           ? sizeof(VOID *) : sizeof(double)));
  } else {
    pool->firstblock = (VOID **) malloc(pool->itemsperblock
                                        * pool->itembytes
                                        + (sizeof(VOID *) >= sizeof(double)
                                           ? sizeof(VOID *) : sizeof(double)));
  }
  if (pool->firstblock == (VOID **) NULL) {
    printf("Error:  Out of memory.\n");
    exit(1);
  }
  /* Set the next block pointer to NULL. */
  *(pool->firstblock) = (VOID *) NULL;
  poolrestart(pool);
}

/*****************************************************************************/
/*                                                                           */
/*  poolrestart()   Deallocate all items in a pool.                          */
/*                                                                           */
/*  The pool is returned to its starting state, except that no memory is     */
/*  freed to the operating system.  Rather, the previously allocated blocks  */
/*  are ready to be reused.                                                  */
/*                                                                           */
/*****************************************************************************/

void poolrestart(pool)
struct memorypool *pool;
{
  VOID **pointeritem;
  double *realitem;

  pool->items = 0;
  pool->maxitems = 0;

  /* Set the currently active block. */
  pool->nowblock = pool->firstblock;
  /* Find the first item in the pool.  Increment by the size of */
  /*   (VOID *) or double, whichever is larger.                 */
  pool->nextitem = sizeof(VOID *) >= sizeof(double)
                 ? (VOID *) (pool->nowblock + 1)
                 : (VOID *) (((double *) pool->nowblock) + 1);
  if (pool->alignwithnull) {
    if (pool->itemwordtype == POINTER) {
      /* Align the item with NULL modulo `itemwords' pointer-sized words. */
      pointeritem = (VOID **) pool->nextitem;
      pointeritem += pool->itemwords
                   - ((pointeritem - (VOID **) NULL) % pool->itemwords);
      pool->nextitem = (VOID *) pointeritem;
    } else {
      /* Align the item with NULL modulo `itemwords' double-sized words. */
      realitem = (double *) pool->nextitem;
      realitem += pool->itemwords
                - ((realitem - (double *) NULL) % pool->itemwords);
      pool->nextitem = (VOID *) realitem;
    }
  }
  /* There are lots of unallocated items left in this block. */
  pool->unallocateditems = pool->itemsperblock;
  /* The stack of deallocated items is empty. */
  pool->deaditemstack = (VOID *) NULL;
}

/*****************************************************************************/
/*                                                                           */
/*  pooldeinit()   Free to the operating system all memory taken by a pool.  */
/*                                                                           */
/*****************************************************************************/

void pooldeinit(pool)
struct memorypool *pool;
{
  while (pool->firstblock != (VOID **) NULL) {
    pool->nowblock = (VOID **) *(pool->firstblock);
    free(pool->firstblock);
    pool->firstblock = pool->nowblock;
  }
}

/*****************************************************************************/
/*                                                                           */
/*  poolalloc()   Allocate space for an item.                                */
/*                                                                           */
/*****************************************************************************/

VOID *poolalloc(pool)
struct memorypool *pool;
{
  VOID *newitem;
  VOID **newblock;
  VOID **pointeritem;
  double *realitem;

  /* First check the linked list of dead items.  If the list is not   */
  /*   empty, allocate an item from the list rather than a fresh one. */
  if (pool->deaditemstack != (VOID *) NULL) {
    newitem = pool->deaditemstack;               /* Take first item in list. */
    pool->deaditemstack = * (VOID **) pool->deaditemstack;
  } else {
    /* Check if there are any free items left in the current block. */
    if (pool->unallocateditems == 0) {
      /* Check if another block must be allocated. */
      if (*(pool->nowblock) == (VOID *) NULL) {
        /* Allocate a new block of items, pointed to by the previous block. */
        if (pool->alignwithnull) {
          newblock = (VOID **) malloc((pool->itemsperblock + 1)
                                      * pool->itembytes
                                      + (sizeof(VOID *) >= sizeof(double)
                                         ? sizeof(VOID *) : sizeof(double)));
        } else {
          newblock = (VOID **) malloc(pool->itemsperblock
                                      * pool->itembytes
                                      + (sizeof(VOID *) >= sizeof(double)
                                         ? sizeof(VOID *) : sizeof(double)));
        }
        if (newblock == (VOID **) NULL) {
          printf("Error:  Out of memory.\n");
          exit(1);
        }
        *(pool->nowblock) = (VOID *) newblock;
        /* The next block pointer is NULL. */
        *newblock = (VOID *) NULL;
      }
      /* Move to the new block. */
      pool->nowblock = (VOID **) *(pool->nowblock);
      /* Find the first item in the block.  Increment by the size of */
      /*   (VOID *) or double, whichever is larger.                  */
      pool->nextitem = sizeof(VOID *) >= sizeof(double)
                     ? (VOID *) (pool->nowblock + 1)
                     : (VOID *) (((double *) pool->nowblock) + 1);
      if (pool->alignwithnull) {
        /* Align first item of block with NULL modulo `itemwords'. */
        if (pool->itemwordtype == POINTER) {
          pointeritem = (VOID **) pool->nextitem;
          pointeritem += pool->itemwords
                       - ((pointeritem - (VOID **) NULL) % pool->itemwords);
          pool->nextitem = (VOID *) pointeritem;
        } else {
          realitem = (double *) pool->nextitem;
          realitem += pool->itemwords
                    - ((realitem - (double *) NULL) % pool->itemwords);
          pool->nextitem = (VOID *) realitem;
        }
      }
      /* There are lots of unallocated items left in this block. */
      pool->unallocateditems = pool->itemsperblock;
    }
    /* Allocate a new item. */
    newitem = pool->nextitem;
    /* Advance `nextitem' pointer to next free item in block. */
    if (pool->itemwordtype == POINTER) {
      pool->nextitem = (VOID *) ((VOID **) pool->nextitem + pool->itemwords);
    } else {
      pool->nextitem = (VOID *) ((double *) pool->nextitem + pool->itemwords);
    }
    pool->unallocateditems--;
    pool->maxitems++;
  }
  pool->items++;
  return newitem;
}

/*****************************************************************************/
/*                                                                           */
/*  pooldealloc()   Deallocate space for an item.                            */
/*                                                                           */
/*  The deallocated space is stored in a queue for later reuse.  This may    */
/*  cause unaligned accesses on machines whose pointer sizes are larger than */
/*  their double sizes.  This typically causes an operating system trap and  */
/*  slightly slower performance, but it happens infrequently enough that     */
/*  I don't mind, especially as the problem usually arises only if single    */
/*  precision floating-point arithmetic is used.  One could eliminate the    */
/*  unaligned accesses by padding `point's so their length is always a       */
/*  multiple of the pointer size, but I don't think it's worth the wasted    */
/*  memory.                                                                  */
/*                                                                           */
/*  Unfortunately, some operating systems have the temerity to print a       */
/*  warning message to the terminal for every unaligned access.  There's     */
/*  usually a command that will disable such bad behavior.                   */
/*                                                                           */
/*****************************************************************************/

void pooldealloc(pool, dyingitem)
struct memorypool *pool;
VOID *dyingitem;
{
  /* Push freshly killed item onto stack. */
  *((VOID **) dyingitem) = pool->deaditemstack;
  pool->deaditemstack = dyingitem;
  pool->items--;
}

/*****************************************************************************/
/*                                                                           */
/*  traversalinit()   Prepare to traverse the entire list of items.          */
/*                                                                           */
/*  This routine is used in conjunction with traverse().                     */
/*                                                                           */
/*****************************************************************************/

void traversalinit(pool)
struct memorypool *pool;
{
  VOID **pointeritem;
  double *realitem;

  /* Begin the traversal in the first block. */
  pool->pathblock = pool->firstblock;
  /* Find the first item in the block.  Increment by the size of */
  /*   (VOID *) or double, whichever is larger.                  */
  pool->pathitem = sizeof(VOID *) >= sizeof(double)
                 ? (VOID *) (pool->pathblock + 1)
                 : (VOID *) (((double *) pool->pathblock) + 1);
  if (pool->alignwithnull) {
    /* Find an aligned item. */
    if (pool->itemwordtype == POINTER) {
      pointeritem = (VOID **) pool->pathitem;
      pointeritem += pool->itemwords
                   - ((pointeritem - (VOID **) NULL) % pool->itemwords);
      pool->pathitem = (VOID *) pointeritem;
    } else {
      realitem = (double *) pool->pathitem;
      realitem += pool->itemwords
                - ((realitem - (double *) NULL) % pool->itemwords);
      pool->pathitem = (VOID *) realitem;
    }
  }
  /* Set the number of items left in the current block. */
  pool->pathitemsleft = pool->itemsperblock;
}

/*****************************************************************************/
/*                                                                           */
/*  traverse()   Find the next item in the list.                             */
/*                                                                           */
/*  This routine is used in conjunction with traversalinit().  Be forewarned */
/*  that this routine successively returns all items in the list, including  */
/*  deallocated ones on the deaditemqueue.  It's up to you to figure out     */
/*  which ones are actually dead.  Why?  I don't want to allocate extra      */
/*  space just to demarcate dead items.  It can usually be done more         */
/*  space-efficiently by a routine that knows something about the structure  */
/*  of the item.                                                             */
/*                                                                           */
/*****************************************************************************/

VOID *traverse(pool)
struct memorypool *pool;
{
  VOID *newitem;
  VOID **pointeritem;
  double *realitem;

  /* Stop upon exhausting the list of items. */
  if (pool->pathitem == pool->nextitem) {
    return (VOID *) NULL;
  }
  /* Check whether any untraversed items remain in the current block. */
  if (pool->pathitemsleft == 0) {
    /* Find the next block. */
    pool->pathblock = (VOID **) *(pool->pathblock);
    /* Find the first item in the block.  Increment by the size of */
    /*   (VOID *) or double, whichever is larger.                  */
    pool->pathitem = sizeof(VOID *) >= sizeof(double)
                   ? (VOID *) (pool->pathblock + 1)
                   : (VOID *) (((double *) pool->pathblock) + 1);
    if (pool->alignwithnull) {
      /* Find an aligned item. */
      if (pool->itemwordtype == POINTER) {
        pointeritem = (VOID **) pool->pathitem;
        pointeritem += pool->itemwords
                     - ((pointeritem - (VOID **) NULL) % pool->itemwords);
        pool->pathitem = (VOID *) pointeritem;
      } else {
        realitem = (double *) pool->pathitem;
        realitem += pool->itemwords
                  - ((realitem - (double *) NULL) % pool->itemwords);
        pool->pathitem = (VOID *) realitem;
      }
    }
    /* Set the number of items left in the current block. */
    pool->pathitemsleft = pool->itemsperblock;
  }
  newitem = pool->pathitem;
  /* Find the next item in the block. */
  if (pool->itemwordtype == POINTER) {
    pool->pathitem = (VOID *) ((VOID **) pool->pathitem + pool->itemwords);
  } else {
    pool->pathitem = (VOID *) ((double *) pool->pathitem + pool->itemwords);
  }
  pool->pathitemsleft--;
  return newitem;
}

/**                                                                         **/
/**                                                                         **/
/********* Memory management routines end here                       *********/

void initialgrid(octnodes)
struct memorypool *octnodes;
{
  struct octnode *newnode;
  struct octnode *xadjacent, *yadjacent, *zadjacent;
  struct octnode *xadjacent_y, *yadjacent_y;
  struct octnode *xadjacent_x;
  int xnum, ynum, znum;
  int x, y, z;

  xnum = (int) ((XMAX - XMIN) / STARTBOXSIZE + 0.5);
  if (xnum < 1) {
    xnum = 1;
  }
  ynum = (int) ((YMAX - YMIN) / STARTBOXSIZE + 0.5);
  if (ynum < 1) {
    ynum = 1;
  }
  znum = (int) ((ZMAX - ZMIN) / STARTBOXSIZE + 0.5);
  if (znum < 1) {
    znum = 1;
  }

  xadjacent_x = (struct octnode *) NULL;
  for (x = 0; x < xnum; x++) {
    xadjacent_y = xadjacent_x;
    yadjacent_y = (struct octnode *) NULL;
    for (y = 0; y < ynum; y++) {
      xadjacent = xadjacent_y;
      yadjacent = yadjacent_y;
      zadjacent = (struct octnode *) NULL;
      if (x > 0) {
        xadjacent_y = xadjacent_y->neighbor[1][1];
      }
      for (z = 0; z < znum; z++) {
        newnode = (struct octnode *) poolalloc(octnodes);
        newnode->mytype = LEAF;
        newnode->treedepth = 0;
        newnode->numerator[0] = x;
        newnode->denominator[0] = xnum;
        newnode->numerator[1] = y;
        newnode->denominator[1] = ynum;
        newnode->numerator[2] = z;
        newnode->denominator[2] = znum;

        if (x == 0) {
          newnode->neighbor[0][0] = (struct octnode *) NULL;
        } else {
          newnode->neighbor[0][0] = xadjacent;
          xadjacent->neighbor[0][1] = newnode;
          xadjacent = xadjacent->neighbor[2][1];
        }
        newnode->neighbor[0][1] = (struct octnode *) NULL;

        if (y == 0) {
          newnode->neighbor[1][0] = (struct octnode *) NULL;
        } else {
          newnode->neighbor[1][0] = yadjacent;
          yadjacent->neighbor[1][1] = newnode;
          yadjacent = yadjacent->neighbor[2][1];
        }
        newnode->neighbor[1][1] = (struct octnode *) NULL;

        if (z == 0) {
          newnode->neighbor[2][0] = (struct octnode *) NULL;
          yadjacent_y = newnode;
          if (y == 0) {
            xadjacent_x = newnode;
          }
        } else {
          newnode->neighbor[2][0] = zadjacent;
          zadjacent->neighbor[2][1] = newnode;
        }
        zadjacent = newnode;
        newnode->neighbor[2][1] = (struct octnode *) NULL;
      }
    }
  }
}

void splitnode(octnodes, cutnode)
struct memorypool *octnodes;
struct octnode *cutnode;
{
  struct octnode *child[8];
  struct octnode *friend, *friendchild;
  enum octnodetype friendtype;
  int mydepth, frienddepth;
  int i, j, k, ii;

  mydepth = cutnode->treedepth;
  for (i = 0; i < 8; i++) {
    child[i] = (struct octnode *) poolalloc(octnodes);
    child[i]->mytype = LEAF;
    child[i]->treedepth = mydepth + 1;
  }
  for (i = 0; i < 8; i++) {
    for (j = 0; j < 3; j++) {
      child[i]->denominator[j] = cutnode->denominator[j] * 2;
      child[i]->numerator[j] = cutnode->numerator[j] * 2;
      if ((i & (1 << j)) != 0) {
        child[i]->numerator[j]++;
      }
      child[i]->neighbor[j][(i & (1 << j)) == 0] = child[i ^ (1 << j)];
    }
  }
  for (j = 0; j < 3; j++) {
    for (k = 0; k < 2; k++) {
      friend = cutnode->neighbor[j][k];
      if (friend == (struct octnode *) NULL) {
        friendtype = LEAF;
      } else {
        frienddepth = friend->treedepth;
        friendtype = friend->mytype;
        if (frienddepth < mydepth) {
          splitnode(octnodes, friend);
          friend = cutnode->neighbor[j][k];
        }
      }
      for (i = 0; i < 8; i++) {
        if ((i & (1 << j)) == (k << j)) {
          if (friendtype == LEAF) {
            child[i]->neighbor[j][k] = friend;
          } else {
            ii = i ^ (1 << j);
            friendchild = friend->neighbor[ii >> 1][ii & 1];
            child[i]->neighbor[j][k] = friendchild;
            friendchild->neighbor[j][1 - k] = child[i];
          }
        }
      }
    }
  }
  cutnode->mytype = PARENT;
  for (i = 0; i < 8; i++) {
    cutnode->neighbor[i >> 1][i & 1] = child[i];
  }
}

void quadrefine(octnodes)
struct memorypool *octnodes;
{
  struct octnode *testnode;
  double x0, y0, z0;
  double x1, y1, z1;
  double h;
  double hmaybe;
  double xwidth, ywidth, zwidth;
  int changes;

  do {
    changes = 0;
    traversalinit(octnodes);
    testnode = (struct octnode *) traverse(octnodes);
    while (testnode != (struct octnode *) NULL) {
      if (testnode->mytype == LEAF) {
        x0 = XMIN + ((double) testnode->numerator[0]
                     / (double) testnode->denominator[0]) * (XMAX - XMIN);
        y0 = YMIN + ((double) testnode->numerator[1]
                     / (double) testnode->denominator[1]) * (YMAX - YMIN);
        z0 = ZMIN + ((double) testnode->numerator[2]
                     / (double) testnode->denominator[2]) * (ZMAX - ZMIN);
        x1 = XMIN + (((double) testnode->numerator[0] + 1.0)
                     / (double) testnode->denominator[0]) * (XMAX - XMIN);
        y1 = YMIN + (((double) testnode->numerator[1] + 1.0)
                     / (double) testnode->denominator[1]) * (YMAX - YMIN);
        z1 = ZMIN + (((double) testnode->numerator[2] + 1.0)
                     / (double) testnode->denominator[2]) * (ZMAX - ZMIN);
        h = lookup_h(x0, y0, z0, PERIOD);
        hmaybe = lookup_h(x1, y0, z0, PERIOD);
        h = (h < hmaybe) ? h : hmaybe;
        hmaybe = lookup_h(x0, y1, z0, PERIOD);
        h = (h < hmaybe) ? h : hmaybe;
        hmaybe = lookup_h(x1, y1, z0, PERIOD);
        h = (h < hmaybe) ? h : hmaybe;
        hmaybe = lookup_h(x0, y0, z1, PERIOD);
        h = (h < hmaybe) ? h : hmaybe;
        hmaybe = lookup_h(x1, y0, z1, PERIOD);
        h = (h < hmaybe) ? h : hmaybe;
        hmaybe = lookup_h(x0, y1, z1, PERIOD);
        h = (h < hmaybe) ? h : hmaybe;
        hmaybe = lookup_h(x1, y1, z1, PERIOD);
        h = (h < hmaybe) ? h : hmaybe;
        xwidth = (XMAX - XMIN) / (double) testnode->denominator[0];
        ywidth = (YMAX - YMIN) / (double) testnode->denominator[1];
        zwidth = (ZMAX - ZMIN) / (double) testnode->denominator[2];
        if ((h < xwidth) || (h < ywidth) || (h < zwidth)) {
          splitnode(octnodes, testnode);
          changes = 1;
        }
      }
      testnode = (struct octnode *) traverse(octnodes);
    }
    printf("%d boxes so far.\n", octnodes->items);
  } while (changes);
}

int countquadnodes(octnodes)
struct memorypool *octnodes;
{
  struct octnode *outnode;
  enum octnodetype xtype, ytype, ztype;
  enum octnodetype xytype, yztype, zxtype;
  int nodenumber;

  nodenumber = 0;
  traversalinit(octnodes);
  outnode = (struct octnode *) traverse(octnodes);
  while (outnode != (struct octnode *) NULL) {
    if (outnode->mytype == LEAF) {

      nodenumber++;

      if (outnode->neighbor[0][0] == (struct octnode *) NULL) {
        xtype = LEAF;
        xytype = LEAF;
      } else {
        xtype = outnode->neighbor[0][0]->mytype;
        if ((outnode->neighbor[0][0]->neighbor[1][0]
             == (struct octnode *) NULL)
            || (outnode->neighbor[0][0]->treedepth < outnode->treedepth)) {
          xytype = LEAF;
        } else {
          xytype = outnode->neighbor[0][0]->neighbor[1][0]->mytype;
        }
      }
      if (outnode->neighbor[1][0] == (struct octnode *) NULL) {
        ytype = LEAF;
        yztype = LEAF;
      } else {
        ytype = outnode->neighbor[1][0]->mytype;
        if ((outnode->neighbor[1][0]->neighbor[2][0]
             == (struct octnode *) NULL)
            || (outnode->neighbor[1][0]->treedepth < outnode->treedepth)) {
          yztype = LEAF;
        } else {
          yztype = outnode->neighbor[1][0]->neighbor[2][0]->mytype;
        }
      }
      if (outnode->neighbor[2][0] == (struct octnode *) NULL) {
        ztype = LEAF;
        zxtype = LEAF;
      } else {
        ztype = outnode->neighbor[2][0]->mytype;
        if ((outnode->neighbor[2][0]->neighbor[0][0]
             == (struct octnode *) NULL)
            || (outnode->neighbor[2][0]->treedepth < outnode->treedepth)) {
          zxtype = LEAF;
        } else {
          zxtype = outnode->neighbor[2][0]->neighbor[0][0]->mytype;
        }
      }

      if (xtype == PARENT) {
        nodenumber++;
      }
      if (ytype == PARENT) {
        nodenumber++;
      }
      if (ztype == PARENT) {
        nodenumber++;
      }
      if ((xtype == PARENT) || (ytype == PARENT) || (xytype == PARENT)) {
        nodenumber++;
      }
      if ((ytype == PARENT) || (ztype == PARENT) || (yztype == PARENT)) {
        nodenumber++;
      }
      if ((ztype == PARENT) || (xtype == PARENT) || (zxtype == PARENT)) {
        nodenumber++;
      }

      if (outnode->neighbor[0][1] == (struct octnode *) NULL) {
        nodenumber++;
        if (outnode->neighbor[1][1] == (struct octnode *) NULL) {
          nodenumber++;
          if (outnode->neighbor[2][1] == (struct octnode *) NULL) {
            nodenumber++;
          }
        }
        if (outnode->neighbor[2][1] == (struct octnode *) NULL) {
          nodenumber++;
        }
        if (ytype == PARENT) {
          nodenumber++;
        }
        if (ztype == PARENT) {
          nodenumber++;
        }
      }
      if (outnode->neighbor[1][1] == (struct octnode *) NULL) {
        nodenumber++;
        if (outnode->neighbor[2][1] == (struct octnode *) NULL) {
          nodenumber++;
        }
        if (xtype == PARENT) {
          nodenumber++;
        }
        if (ztype == PARENT) {
          nodenumber++;
        }
      }
      if (outnode->neighbor[2][1] == (struct octnode *) NULL) {
        nodenumber++;
        if (xtype == PARENT) {
          nodenumber++;
        }
        if (ytype == PARENT) {
          nodenumber++;
        }
      }
    }

    outnode = (struct octnode *) traverse(octnodes);
  }

  return nodenumber;
}

void writequadnodes(octnodes, nodefilename)
struct memorypool *octnodes;
char *nodefilename;
{
  struct octnode *outnode;
  FILE *nodefile;
  enum octnodetype xtype, ytype, ztype;
  enum octnodetype xytype, yztype, zxtype;
  double smallx, smally, smallz;
  double middlex, middley, middlez;
  double bigx, bigy, bigz;
  int nodenumber;

  nodefile = fopen(nodefilename, "w");
  if (nodefile == (FILE *) NULL) {
    printf("  Error:  Cannot create file %s.\n", nodefilename);
    exit(1);
  }
  fprintf(nodefile, "%d  3  0  0\n", countquadnodes(octnodes));
  nodenumber = 0;
  traversalinit(octnodes);
  outnode = (struct octnode *) traverse(octnodes);
  while (outnode != (struct octnode *) NULL) {
    if (outnode->mytype == LEAF) {

      smallx = ((double) outnode->numerator[0]
                / (double) outnode->denominator[0]) * (XMAX - XMIN) + XMIN;
      smally = ((double) outnode->numerator[1]
                / (double) outnode->denominator[1]) * (YMAX - YMIN) + YMIN;
      smallz = ((double) outnode->numerator[2]
                / (double) outnode->denominator[2]) * (ZMAX - ZMIN) + ZMIN;
      middlex = ((double) (2 * outnode->numerator[0] + 1)
                 / (double) (2 * outnode->denominator[0]))
              * (XMAX - XMIN) + XMIN;
      middley = ((double) (2 * outnode->numerator[1] + 1)
                 / (double) (2 * outnode->denominator[1]))
              * (YMAX - YMIN) + YMIN;
      middlez = ((double) (2 * outnode->numerator[2] + 1)
                 / (double) (2 * outnode->denominator[2]))
              * (ZMAX - ZMIN) + ZMIN;
      bigx = ((double) (outnode->numerator[0] + 1)
              / (double) outnode->denominator[0]) * (XMAX - XMIN) + XMIN;
      bigy = ((double) (outnode->numerator[1] + 1)
              / (double) outnode->denominator[1]) * (YMAX - YMIN) + YMIN;
      bigz = ((double) (outnode->numerator[2] + 1)
              / (double) outnode->denominator[2]) * (ZMAX - ZMIN) + ZMIN;

      fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
              smallx, smally, smallz);

      if (outnode->neighbor[0][0] == (struct octnode *) NULL) {
        xtype = LEAF;
        xytype = LEAF;
      } else {
        xtype = outnode->neighbor[0][0]->mytype;
        if ((outnode->neighbor[0][0]->neighbor[1][0]
             == (struct octnode *) NULL)
            || (outnode->neighbor[0][0]->treedepth < outnode->treedepth)) {
          xytype = LEAF;
        } else {
          xytype = outnode->neighbor[0][0]->neighbor[1][0]->mytype;
        }
      }
      if (outnode->neighbor[1][0] == (struct octnode *) NULL) {
        ytype = LEAF;
        yztype = LEAF;
      } else {
        ytype = outnode->neighbor[1][0]->mytype;
        if ((outnode->neighbor[1][0]->neighbor[2][0]
             == (struct octnode *) NULL)
            || (outnode->neighbor[1][0]->treedepth < outnode->treedepth)) {
          yztype = LEAF;
        } else {
          yztype = outnode->neighbor[1][0]->neighbor[2][0]->mytype;
        }
      }
      if (outnode->neighbor[2][0] == (struct octnode *) NULL) {
        ztype = LEAF;
        zxtype = LEAF;
      } else {
        ztype = outnode->neighbor[2][0]->mytype;
        if ((outnode->neighbor[2][0]->neighbor[0][0]
             == (struct octnode *) NULL)
            || (outnode->neighbor[2][0]->treedepth < outnode->treedepth)) {
          zxtype = LEAF;
        } else {
          zxtype = outnode->neighbor[2][0]->neighbor[0][0]->mytype;
        }
      }

      if (xtype == PARENT) {
        fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                smallx, middley, middlez);
      }
      if (ytype == PARENT) {
        fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                middlex, smally, middlez);
      }
      if (ztype == PARENT) {
        fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                middlex, middley, smallz);
      }
      if ((xtype == PARENT) || (ytype == PARENT) || (xytype == PARENT)) {
        fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                smallx, smally, middlez);
      }
      if ((ytype == PARENT) || (ztype == PARENT) || (yztype == PARENT)) {
        fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                middlex, smally, smallz);
      }
      if ((ztype == PARENT) || (xtype == PARENT) || (zxtype == PARENT)) {
        fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                smallx, middley, smallz);
      }

      if (outnode->neighbor[0][1] == (struct octnode *) NULL) {
        fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                bigx, smally, smallz);
        if (outnode->neighbor[1][1] == (struct octnode *) NULL) {
          fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                  bigx, bigy, smallz);
          if (outnode->neighbor[2][1] == (struct octnode *) NULL) {
            fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                    bigx, bigy, bigz);
          }
        }
        if (outnode->neighbor[2][1] == (struct octnode *) NULL) {
          fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                  bigx, smally, bigz);
        }
        if (ytype == PARENT) {
          fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                  bigx, smally, middlez);
        }
        if (ztype == PARENT) {
          fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                  bigx, middley, smallz);
        }
      }
      if (outnode->neighbor[1][1] == (struct octnode *) NULL) {
        fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                smallx, bigy, smallz);
        if (outnode->neighbor[2][1] == (struct octnode *) NULL) {
          fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                  smallx, bigy, bigz);
        }
        if (xtype == PARENT) {
          fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                  smallx, bigy, middlez);
        }
        if (ztype == PARENT) {
          fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                  middlex, bigy, smallz);
        }
      }
      if (outnode->neighbor[2][1] == (struct octnode *) NULL) {
        fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                smallx, smally, bigz);
        if (xtype == PARENT) {
          fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                  smallx, middley, bigz);
        }
        if (ytype == PARENT) {
          fprintf(nodefile, "%4d    %.20g  %.20g  %.20g\n", ++nodenumber,
                  middlex, smally, bigz);
        }
      }
    }

    outnode = (struct octnode *) traverse(octnodes);
  }

  fclose(nodefile);
}

int main()
{
  struct memorypool octnodes;

/*  initlookup_h("/afs/cs/project/quake/3/jxu/lookuph/cut1_final.dat"); */
  poolinit(&octnodes, sizeof(struct octnode), OCTNODEPERBLOCK, POINTER, 0);
  initialgrid(&octnodes);
  printf("Refining the initial grid of %d boxes.\n", octnodes.items);
  quadrefine(&octnodes);
  writequadnodes(&octnodes, "try.0.node");

  return 0;
}