# $Id: slice3.py,v 1.8 2003/09/30 23:43:47 travo Exp $
# Copyright (c) 1996, 1997, The Regents of the University of California.
# All rights reserved.  See Legal.htm for full text and disclaimer.

#  SLICE3.PY
# find 2D slices of a 3D hexahedral mesh

#  $Id: slice3.py,v 1.8 2003/09/30 23:43:47 travo Exp $
#

# Change (ZCM 12/4/96) Apparently _draw3_list, which is global in
# pl3d.py, can be fetched from there, but once this has been done,
# assignments to it over there are not reflected in the copy here.
# This has been fixed by creating an access function.

from Numeric import *
from shapetest import *
from types import *
from pl3d import *
from arrayfns import *
from gistC import *

 #
 # Caveats:
 # (A) Performance is reasonably good, but may still be a factor of
 #     several slower than what could be achieved in compiled code.
 # (B) Only a simple in-memory mesh model is implemented here.
 #     However, hooks are supplied for more interesting possibilities,
 #     such as a large binary file resident mesh data base.
 # (C) There is a conceptual difficulty with _walk3 for the case
 #     of a quad face all four of whose edges are cut by the slicing
 #     plane.  This can only happen when two opposite corners are
 #     above and the other two below the slicing plane.  There are
 #     three possible ways to connect the four intersection points in
 #     two pairs: (1) // (2) \\ and (3) X  There is a severe problem
 #     with (1) and (2) in that a consistent decision must be made
 #     when connecting the points on the two cells which share the
 #     face - that is, each face must carry information on which way
 #     it is triangulated.  For a regular 3D mesh, it is relatively
 #     easy to come up with a consistent scheme for triangulating faces,
 #     but for a general unstructured mesh, each face itself must carry
 #     this information.  This presents a huge challenge for data flow,
 #     which I don't believe is worthwhile.  Because the X choice is
 #     unique, and I don't see why we shouldn't use it here.
 #     For contouring routines, we reject the X choice on esthetic
 #     grounds, and perhaps that will prove to be the case here as
 #     well - but I believe we should try the simple way out first.
 #     In this case, we are going to be filling these polygons with
 #     a color representing a function value in the cell.  Since the
 #     adjacent cells should have nearly the same values, the X-traced
 #     polygons will have nearly the same color, and I doubt there will
 #     be an esthetic problem.  Anyway, the slice3 implemented
 #     below produces the unique X (bowtied) polygons, rather than
 #     attempting to choose between // or \\ (non-bowtied) alternatives.
 #     Besides, in the case of contours, the trivial alternating
 #     triangulation scheme is just as bad esthetically as every
 #     zone triangulated the same way!

def plane3 (normal, point) :

   """
   plane3(normal, point)
         or plane3([nx,ny,nz], [px,py,pz])

     returns [nx,ny,nz,pp] for the specified plane.
   """

   # the normal doesn't really need to be normalized, but this
   # has the desirable side effect of blowing up if normal==0
   newnorm = zeros (4, Float)
   newnorm [0:3] = normal / sqrt (sum (normal*normal))
   newnorm [3] = sum (multiply (normal, point))
   return newnorm

_Mesh3Error = "Mesh3Error"

def mesh3 (x, y = None, z = None, ** kw) :

   """
    mesh3(x,y,z)
         or mesh3(x,y,z, funcs = [f1,f2,...])
         or mesh3(xyz, funcs = [f1,f2,...])
         or mesh3(nxnynz, dxdydz, x0y0z0, funcs = [f1,f2,...])

     make mesh3 argument for slice3, xyz3, getv3, etc., functions.
     X, Y, and Z are each 3D coordinate arrays.  The optional F1, F2,
     etc. are 3D arrays of function values (e.g. density, temperature)
     which have one less value along each dimension than the coordinate
     arrays.  The "index" of each zone in the returned mesh3 is
     the index in these cell-centered Fi arrays, so every index from
     one through the total number of cells indicates one real cell.
     The Fi arrays can also have the same dimensions as X, Y, or Z
     in order to represent point-centered quantities.

     If X has four dimensions and the length of the first is 3, then
     it is interpreted as XYZ (which is the quantity actually stored
     in the returned cell list).

     If X is a vector of 3 integers, it is interpreted as [nx,ny,nz]
     of a uniform 3D mesh, and the second and third arguments are
     [dx,dy,dz] and [x0,y0,z0] respectively.  (DXDYDZ represent the
     size of the entire mesh, not the size of one cell, and NXNYNZ are
     the number of cells, not the number of points.)

     Added by ZCM 1/13/97: if x, y, and z are one-dimensional of
     the same length and if the keyword verts exists and yields
     an NCELLS by 8 integer array, then we have an unstructured
     rectangular mesh, and the subscripts of cell i's vertices
     are verts[i, 0:8].

     Added by ZCM 10/10/97: if x, y, and z are one-dimensional
     of the same length or not, and verts does not exist, then
     we have a structured reectangular mesh with unequally spaced
     nodes.

     Other sorts of meshes are possible -- a mesh which lives
     in a binary file is an obvious example -- which would need
     different workers for xyz3, getv3, getc3, and iterator3
     iterator3_rect may be more general than the other three;
     as long as the cell dimensions are the car of the list
     which is the 2nd car of m3, it will work.
   """

   dims = shape (x)
   if len (dims) == 1 and y != None and len (x) == len (y) \
      and z != None and len(x) == len (z) and kw.has_key ("verts") :
      virtuals = [xyz3_irreg, getv3_irreg,
                  getc3_irreg, iterator3_irreg]
      dims = kw ["verts"]
      if type (dims) != ListType :
         m3 = [virtuals, [dims, array ( [x, y, z])], []]
      else : # Irregular mesh with more than one cell type
         sizes = ()
         for nv in dims :
            sizes = sizes + (shape (nv) [0],) # no. cells of this type
         totals = [sizes [0]]
         for i in range (1, len (sizes)) :
            totals.append (totals [i - 1] + sizes [i]) #total cells so far
         m3 = [virtuals, [dims, array ( [x, y, z]), sizes, totals], []]
      if kw.has_key ("funcs") :
         funcs = kw ["funcs"]
      else :
         funcs = []
      i = 0
      for f in funcs:
         if len (f) != len (x) and len (f) != shape (dims) [0] :
            # if vertex-centered, f must be same size as x.
            # if zone centered, its length must match number of cells.
            raise _Mesh3Error, "F" + `i` + " is not a viable 3D cell value"
         m3 [2] = m3 [2] + [f]
         i = i + 1
      return m3

   virtuals = [xyz3_rect, getv3_rect, getc3_rect, iterator3_rect]
   if len (dims) == 4 and dims [0] == 3 and min (dims) >= 2 :
      xyz = x
      dims = dims [1:4]
   elif len (dims) == 1 and len (x) == 3 and type (x [0]) == IntType \
      and y != None and z != None and len (y) == len (z) == 3 :
      xyz = array ([y, z])
      dims = (1 + x [0], 1 + x [1], 1 + x [2])
      virtuals [0] = xyz3_unif
   elif len (dims) == 1 and y != None and z != None and len (y.shape) == 1 \
      and len (z.shape) == 1 and x.typecode () == y.typecode () == \
      z.typecode () == Float : 
      # regular mesh with unequally spaced points
      dims = array ( [len (x), len (y), len (z)], Int)
      xyz = [x, y, z] # has to be a list since could be different lengths
      virtuals [0] = xyz3_unif
   else :
      if len (dims) != 3 or min (dims) < 2 or \
         y == None or len (shape (y)) != 3 or shape (y) != dims or \
         z == None or len (shape (z)) != 3 or shape (z) != dims:
         raise _Mesh3Error, "X,Y,Z are not viable 3D coordinate mesh arrays"
      xyz = array ( [x, y, z])
   dim_cell = (dims [0] - 1, dims [1] - 1, dims [2] - 1)
   m3 = [virtuals, [dim_cell, xyz], []]
   if kw.has_key ("funcs") :
      funcs = kw ["funcs"]
   else :
      funcs = []
   i = 0
   for f in funcs:
      if len (f.shape) == 3 and \
         ( (f.shape [0] == dims [0] and f.shape [1] == dims [1] and
            f.shape [2] == dims [2]) or (f.shape [0] == dim_cell [0] and
            f.shape [1] == dim_cell [1] and f.shape [2] == dim_cell [2])) :
         m3 [2] = m3 [2] + [f]
         i = i + 1
      else :
         raise _Mesh3Error, "F" + `i` + " is not a viable 3D cell value"

   return m3

 # Ways that a list of polygons can be extracted:
 # Basic idea:
 #   (1) At each *vertex* of the cell list, a function value is defined.
 #       This is the "slicing function", perhaps the equation of a plane,
 #       perhaps some other vertex-centered function.
 #   (2) The slice3 routine returns a list of cells for which the
 #       function value changes sign -- that is, cells for which some
 #       vertices have positive function values, and some negative.
 #       The function values and vertex coordinates are also returned.
 #   (3) The slice3 routine computes the points along the *edges*
 #       of each cell where the function value is zero (assuming linear
 #       variation along each edge).  These points will be vertices of
 #       the polygons.  The routine also sorts the vertices into cyclic
 #       order.
 #   (4) A "color function" can be used to assign a "color" or other
 #       value to each polygon.  If this function depends only on the
 #       coordinates of the polygon vertices (e.g.- 3D lighting), then
 #       the calculation can be done elsewhere.  There are two other
 #       possibilities:  The color function might be a cell-centered
 #       quantity, or a vertex-centered quantity (like the slicing
 #       function) on the mesh.  In these cases, slice3 already
 #       has done much of the work, since it "knows" cell indices,
 #       edge interpolation coefficients, and the like.
 #
 # There are two particularly important cases:
 # (1) Slicing function is a plane, coloring function is either a
 #     vertex or cell centered mesh function.  Coloring function
 #     might also be a *function* of one or more of the predefined
 #     mesh functions.  If you're eventually going to sweep the whole
 #     mesh, you want to precalculate it, otherwise on-the-fly might
 #     be better.
 # (2) Slicing function is a vertex centered mesh function,
 #     coloring function is 3D shading (deferred).
 #
 # fslice(m3, vertex_list)
 # vertex_list_iterator(m3, vertex_list, mesh3)
 # fcolor(m3, vertex_list, fslice_1, fslice_2)
 #   the coloring function may need the value of fslice at the vertices
 #   in order to compute the color values by interpolation
 # two "edge functions": one to detect edges where sign of fslice changes,
 #   second to interpolate for fcolor
 #   second to interpolate for fcolor
 #
 # slice3(m3, fslice, &nverts, &xyzverts, <fcolor>)

_Slice3Error = "Slice3Error"


def slice3 (m3, fslice, nverts, xyzverts, * args, ** kw) :

   """
   slice3 (m3, fslice, nverts, xyzverts)
         or color_values= slice3(m3, fslice, nverts, xyzverts, fcolor)
         or color_values= slice3(m3, fslice, nverts, xyzverts, fcolor, 1)

     slice the 3D mesh M3 using the slicing function FSLICE, returning
     the list [NVERTS, XYZVERTS, color].  Note that it is impossible to
     pass arguments as addresses, as yorick does in this routine.
     NVERTS is the number of vertices in each polygon of the slice, and
     XYZVERTS is the 3-by-sum(NVERTS) list of polygon vertices.  If the
     FCOLOR argument is present, the values of that coloring function on
     the polygons are returned as the value of the slice3 function
     (numberof(color_values) == numberof(NVERTS) == number of polygons).

     If the slice function FSLICE is a function, it should be of the
     form:
        func fslice(m3, chunk)
     returning a list of function values on the specified chunk of the
     mesh m3.  The format of chunk depends on the type of m3 mesh, so
     you should use only the other mesh functions xyz3 and getv3 which
     take m3 and chunk as arguments.  The return value of fslice should
     have the same dimensions as the return value of getv3; the return
     value of xyz3 has an additional first dimension of length 3.

     If FSLICE is a list of 4 numbers, it is taken as a slicing plane
     with the equation FSLICE(+:1:3)*xyz(+)-FSLICE(4), as returned by
     plane3.

     If FSLICE is a single integer, the slice will be an isosurface for
     the FSLICEth variable associated with the mesh M3.  In this case,
     the keyword value= must also be present, representing the value
     of that variable on the isosurface.

     If FCOLOR is nil, slice3 returns nil.  If you want to color the
     polygons in a manner that depends only on their vertex coordinates
     (e.g.- by a 3D shading calculation), use this mode.

     If FCOLOR is a function, it should be of the form:
        func fcolor(m3, cells, l, u, fsl, fsu, ihist)
     returning a list of function values on the specified cells of the
     mesh m3.  The cells argument will be the list of cell indices in
     m3 at which values are to be returned.  l, u, fsl, fsu, and ihist
     are interpolation coefficients which can be used to interpolate
     from vertex centered values to the required cell centered values,
     ignoring the cells argument.  See getc3 source code.
     The return values should always have dimsof(cells).

     If FCOLOR is a single integer, the slice will be an isosurface for
     the FCOLORth variable associated with the mesh M3.

     If the optional argument after FCOLOR is non-nil and non-zero,
     then the FCOLOR function is called with only two arguments:
        func fcolor(m3, cells)

     The keyword argument NODE, if present and nonzero, is a signal
     to return node-centered values rather than cell-centered
     values. (ZCM 4/16/97)

   """

   global _poly_permutations

   iso_index = None
   if type (fslice) != FunctionType :
      if not kw.has_key ("value") and not is_scalar (fslice) and \
         len (shape (fslice)) == 1 and len (fslice) == 4 :
         normal = fslice [0:3]
         projection = fslice [3]
         fslice = _plane_slicer
      elif is_scalar (fslice) and type (fslice) == IntType :
         if not kw.has_key ("value") :
            raise _Slice3Error, \
               "value= keyword required when FSLICE is mesh variable"
         _value = kw ["value"]
         iso_index = fslice
         fslice = _isosurface_slicer
      else :
         raise _Slice3Error, \
            "illegal form of FSLICE argument, try help,slice3"

   if kw.has_key ("node") :
      node = kw ["node"]
   else :
      node = 0

   # will need cell list if fcolor function to be computed
   need_clist = len (args) > 0
   if len (args) > 1 :
      nointerp = args [1]
   else :
      nointerp = None

   if need_clist :
      fcolor = args [0]
      if fcolor == None :
         need_clist = 0
   else :
      fcolor = None
   
   # test the different possibilities for fcolor
   if need_clist and type (fcolor) != FunctionType :
      if not is_scalar (fcolor) or type (fcolor) != IntType :
         raise _Slice3Error, \
            "illegal form of FCOLOR argument, try help,slice3"

   # chunk up the m3 mesh and evaluate the slicing function to
   # find those cells cut by fslice==0
   # chunking avoids potentially disastrously large temporaries
   got_xyz = 0
   ntotal = 0
   # The following are used only for an irregular mesh, to
   # give the sizes of each portion of the mesh.
   ntotal8 = 0
   ntotal6 = 0
   ntotal5 = 0
   ntotal4 = 0
   # The following are used only for an irregular mesh, to
   # give the indices of the different types of chunk in the
   # results list.
   i8 = []
   i6 = []
   i5 = []
   i4 = []
   itot = [i4, i5, i6, i8]
   nchunk = 0
   results = []
   chunk = iterator3 (m3)
   cell_offsets = [0, 0, 0, 0]
   while chunk != None :

      # get the values of the slicing function at the vertices of
      # this chunk
      if fslice == _isosurface_slicer :
         fs = fslice (m3, chunk, iso_index, _value)
         # an isosurface slicer brings back a list [vals, None]
         # where vals is simply an array of the values of the
         # iso_index'th function on the vertices of the specified
         # chunk, or a triple, consisting of the array of
         # values, an array of relative cell numbers in the
         # chunk, and an offset to add to the preceding to
         # get absolute cell numbers.
      elif fslice == _plane_slicer :
         fs = fslice (m3, chunk, normal, projection)
         # In the case of a plane slice, fs is a list [vals, _xyz3]
         # (or [ [vals, clist, cell_offset], _xyz3] in the irregular case)
         # where _xyz3 is the array of vertices of the chunk. _xyz3
         # is ncells by 3 by something (in the irregular case),
         # ncells by 3 by 2 by 2 by 2 in the regular case,
         # and 3 by ni by nj by nk otherwise. vals will be
         # the values of the projections of the corresponding
         # vertex on the normal to the plane, positive if in
         # front, and negative if in back.
      else :
         fs = fslice (m3, chunk)
      if node == 1 and fcolor != None and fcolor != FunctionType :
         # need vertex-centered data
         col = getv3 (fcolor, m3, chunk)
         if type (col) == ListType :
            col = col [0]
      else :
         col = None
      # ZCM 2/24/97 Elimination of _xyz3 as a global necessitates the following:
      # (_xyz3 comes back as the last element of the list fs)
      _xyz3 = fs [1]
      fs = fs [0]
      irregular = type (fs) == ListType
      if irregular :
         cell_offset = fs [2]

      # will need cell list if fslice did not compute xyz
      got_xyz = _xyz3 != None
      need_clist = need_clist or not got_xyz

      # If the m3 mesh is totally unstructured, the chunk should be
      # arranged so that fslice returns an ncells-by-2-by-2-by-2
      # (or ncells-by-3-by-2 or ncells-by-5 or ncells-by-4) array
      # of vertex values of the slicing function. Note that a
      # chunk of an irregular mesh always consists of just one
      # kind of cell.
      # On the other hand, if the mesh vertices are arranged in a
      # rectangular grid (or a few patches of rectangular grids), the
      # chunk should be the far less redundant rectangular patch.
      if (irregular) :
         # fs is a 2-sequence, of which the first element is an ncells-by-
         # 2-by-2-by-2 (by-3-by-2, by-5, or by-4) array, and the second
         # is the array of corresponding cell numbers.
         # here is the fastest way to generate the required cell list
         dims = shape (fs [0])
         dim1 = dims [0]
         slice3_precision = 0.0
         if len (dims) == 4 : # hex case
            # Note that the sum below will be between 1 and 7
            # precisely if f changes sign in the cell.
            critical_cells = bitwise_and (add.reduce \
               (reshape (ravel (transpose (less (fs [0], slice3_precision))), \
               (8, dim1))), 7)
            if (sum (critical_cells) != 0) :
               clist = take (fs [1], nonzero (critical_cells))
               ntotal8 = ntotal8 + len (clist)
            else :
               clist = None
            i8.append (len (results))
            cell_offsets [3] = cell_offset
         elif len (dims) == 3 : # prism case
            # Note that the sum below will be between 1 and 5
            # precisely if f changes sign in the cell.
            critical_cells = add.reduce \
               (reshape (ravel (transpose (less (fs [0], slice3_precision))), \
               (6, dim1)))
            critical_cells = logical_and (greater (critical_cells, 0),
                                         less (critical_cells, 6))
            if (sum (critical_cells) != 0) :
               clist = take (fs [1], nonzero (critical_cells))
               ntotal6 = ntotal6 + len (clist)
            else :
               clist = None
            i6.append (len (results))
            cell_offsets [2] = cell_offset
         elif dims [1] == 5 : # pyramid case
            # Note that the sum below will be between 1 and 4
            # precisely if f changes sign in the cell.
            critical_cells = add.reduce \
               (reshape (ravel (transpose (less (fs [0], slice3_precision))), \
               (5, dim1)))
            critical_cells = logical_and (greater (critical_cells, 0),
                                         less (critical_cells, 5))
            if (sum (critical_cells) != 0) :
               clist = take (fs [1], nonzero (critical_cells))
               ntotal5 = ntotal5 + len (clist)
            else :
               clist = None
            i5.append (len (results))
            cell_offsets [1] = cell_offset
         else : # tet case
            # Note that the sum below will be between 1 and 3
            # precisely if f changes sign in the cell.
            critical_cells = bitwise_and (add.reduce \
               (reshape (ravel (transpose (less (fs [0], slice3_precision))), \
               (4, dim1))), 3)
            if (sum (critical_cells) != 0) :
               clist = take (fs [1], nonzero (critical_cells))
               ntotal4 = ntotal4 + len (clist)
            else :
               clist = None
            i4.append (len (results))
            cell_offsets [0] = cell_offset
      else :
         dims = shape (fs)
         # fs is an ni-by-nj-by-nk array
         # result of the zcen is 0, 1/8, 2/8, ..., 7/8, or 1
#        slice3_precision = max (ravel (abs (fs))) * (-1.e-12)
         slice3_precision = 0
         clist1 = ravel (zcen_ (zcen_ (zcen_
            (array (less (fs, slice3_precision), Float), 0), 1), 2))
         clist1 = logical_and (less (clist1, .9), greater (clist1, .1))
         if sum (clist1) > 0 :
            clist = nonzero (clist1)
            ntotal = ntotal + len (clist)
         else :
            clist = None
         i8.append (len (results)) # Treat regular case as hex

      if clist != None :
         #  we need to save:
         # (1) the absolute cell indices of the cells in clist
         # (2) the corresponding ncells-by-2-by-2-by-2 (by-3-by-2,
         #     by-5, or by-4) list of fslice
         #     values at the vertices of these cells
         if (irregular) :
            # extract the portions of the data indexed by clist
            fs = take (fs [0], clist)
            if got_xyz :
               _xyz3 = take (_xyz3, clist)
            if col :
               col = take (col, clist)
         else :
            # extract the to_corners portions of the data indexed by clist
            indices = to_corners3 (clist, dims [1], dims [2])
            no_cells = shape (indices) [0]
            indices = ravel (indices)
            fs = reshape (take (ravel (fs), indices),\
               (no_cells, 2, 2, 2))
            if got_xyz :
               new_xyz3 = zeros ( (no_cells, 3, 2, 2, 2), Float )
               new_xyz3 [:, 0, ...] = reshape (take (ravel (_xyz3 [0, ...]),\
                  indices), (no_cells, 2, 2, 2))
               new_xyz3 [:, 1, ...] = reshape (take (ravel (_xyz3 [1, ...]),\
                  indices), (no_cells, 2, 2, 2))
               new_xyz3 [:, 2, ...] = reshape (take (ravel (_xyz3 [2, ...]),\
                  indices), (no_cells, 2, 2, 2))
               _xyz3 = new_xyz3
               del new_xyz3
            if col != None :
               col = reshape (take (ravel (col), indices), (no_cells, 2, 2, 2))
               # NB: col represents node colors, and is only used
               # if those are requested.
         # here, the iterator converts to absolute cell indices without
         # incrementing the chunk
         if (need_clist) :
            clist = iterator3 (m3, chunk, clist)
         else :
            clist = None
         nchunk = nchunk + 1
         need_vert_col = col != None
         results.append ( [clist, fs, _xyz3, col])
      else :
         results.append ( [None, None, None, None])
      chunk = iterator3 (m3, chunk)
      # endwhile chunk != None

   # collect the results of the chunking loop
   if not ntotal and not (ntotal8 + ntotal6 + ntotal5 + ntotal4) :
      return None
   if ntotal : # (regular mesh, but can be handled same as hex)
      ntotal8 = ntotal
      i8 = range (len (results))
      itot [3] = i8
   ntot = [ntotal4, ntotal5, ntotal6, ntotal8]
   new_results = []
   for i in range (len (ntot)) :
      # This loop processes each kind of cell independently,
      # the results to be combined at the end.
      if ntot [i] == 0 : # No cells of type i
         continue
      if need_clist :
         clist = zeros (ntot [i], Int)
         fs = zeros ( (ntot [i], _no_verts [i]), Float )
         if got_xyz :
            xyz = zeros ( (ntot [i], 3, _no_verts [i]), Float )
         else :
            xyz = None
      if need_vert_col :
         col = zeros ( (ntot [i], _no_verts [i]), Float )
      else :
         col = None
      k = 0

     # collect the results of the chunking loop
      for j in range (len (itot [i])) :
         l = k
         k = k + len (results [itot [i] [j]] [0])
         if need_clist :
            clist [l:k] = results [itot [i] [j]] [0]
         fs [l:k] = reshape (results [itot [i] [j]] [1], (k - l, _no_verts [i]))
         if xyz != None :
            xyz [l:k] = reshape (results [itot [i] [j]] [2],
               (k - l, 3, _no_verts [i]))
         if col != None :
            col [l:k] = reshape (results [itot [i] [j]] [3],
               (k - l, _no_verts [i]))
      if not got_xyz :
         # zcm 2/4/97 go to absolute cell list again
         if i > 0 and len (m3 [1]) > 2 :
            adder = m3 [1] [3] [i - 1]
         else :
            adder = 0
         xyz = reshape (xyz3 (m3, clist + adder), (ntot [i], 3, _no_verts [i]))
      # produce the lists of edge intersection points
      # -- generate (nsliced)x12 (9, 8, 6) array of edge mask values
      # (mask non-zero if edge is cut by plane)
      below = less (fs, 0.0)
      # I put the following into C for speed
      mask = find_mask (below, _node_edges [i])
      list = nonzero (mask)
      edges = array (list, copy = 1)
      cells = edges / _no_edges [i]
      edges = edges % _no_edges [i]
      # construct edge endpoint indices in fs, xyz arrays
      # the numbers are the endpoint indices corresponding to
      # the order of the _no_edges [i] edges in the mask array
      lower = take (_lower_vert [i], edges) + _no_verts [i] * cells
      upper = take (_upper_vert [i], edges) + _no_verts [i] * cells
      fsl = take (ravel (fs), lower)
      fsu = take (ravel (fs), upper)
      # following denominator guaranteed non-zero
      denom = fsu - fsl
      fsu = fsu / denom
      fsl = fsl / denom
      new_xyz = zeros ( (len (lower), 3), Float )
      new_xyz [:, 0] = reshape ( (take (ravel (xyz [:, 0]), lower) * fsu - \
         take (ravel (xyz [:, 0]), upper) * fsl), (len (lower),))
      new_xyz [:, 1] = reshape ( (take (ravel (xyz [:, 1]), lower) * fsu - \
         take (ravel (xyz [:, 1]), upper) * fsl), (len (lower),))
      new_xyz [:, 2] = reshape ( (take (ravel (xyz [:, 2]), lower) * fsu - \
         take (ravel (xyz [:, 2]), upper) * fsl), (len (lower),))
      xyz = new_xyz
      del new_xyz
      if col != None :
         # Extract subset of the data the same way
         col = take (ravel (col), lower) * fsu - \
            take (ravel (col), upper) * fsl
      # The xyz array is now the output xyzverts array,
      # but for the order of the points within each cell.

      # give each sliced cell a "pattern index" between 0 and 255
      # (non-inclusive) representing the pattern of its 8 corners
      # above and below the slicing plane
      p2 = left_shift (ones (_no_verts [i], Int) , array (
         [0, 1, 2, 3, 4, 5, 6, 7], Int) [0: _no_verts [i]])
      pattern = transpose (sum (transpose (multiply (below, p2))))

      # broadcast the cell's pattern onto each of its sliced edges
      pattern = take (pattern, list / _no_edges [i])
      # Let ne represent the number of edges of this type of cell,
      # and nv the number of vertices.
      # To each pattern, there corresponds a permutation of the
      # twelve edges so that they occur in the order in which the
      # edges are to be connected.  Let each such permuation be
      # stored as a list of integers from 0 to ne - 1 such that
      # sorting the integers into increasing order rearranges the edges at
      # the corresponding indices into the correct order.  (The position
      # of unsliced edges in the list is arbitrary as long as the sliced
      # edges are in the proper order relative to each other.)
      # Let these permutations be stored in a ne-by-2**nv - 2 array
      # _poly_permutations (see next comment for explanation of 4 * ne):
      pattern = take (ravel (transpose (_poly_permutations [i])), 
         _no_edges [i] * (pattern - 1) + edges) + 4 * _no_edges [i] * cells
      order = argsort (pattern)
      xyz1 = zeros ( (len (order), 3), Float )
      xyz1 [:,0] = take (ravel (xyz [:,0]), order)
      xyz1 [:,1] = take (ravel (xyz [:,1]), order)
      xyz1 [:,2] = take (ravel (xyz [:,2]), order)
      xyz = xyz1
      if col != None :
         col = take (col, order)
      edges = take (edges, order)
      pattern = take (pattern, order)
      # cells(order) is same as cells by construction */

      # There remains only the question of splitting the points in
      # a single cell into multiple disjoint polygons.
      # To do this, we need one more precomputed array: poly_splits
      # should be another ne-by-2**nv - 2 array with values between 0 and 3
      # 0 for each edge on the first part, 1 for each edge on the
      # second part, and so on up to 3 for each edge on the fourth
      # part.  The value on unsliced edges can be anything, say 0.
      # With a little cleverness poly_splits can be combined with
      # _poly_permutations, by putting _poly_permutations =
      # _poly_permutations(as described above) + _no_edges [i]*poly_splits
      # (this doesn't change the ordering of _poly_permutations).
      # I assume this has been done here:
      pattern = pattern / _no_edges [i]
      # now pattern jumps by 4 between cells, smaller jumps within cells
      # get the list of places where a new value begins, and form a
      # new pattern with values that increment by 1 between each plateau
      pattern = dif_ (pattern, 0)
      nz = nonzero (pattern)
      list = zeros (len (nz) + 1, Int)
      list [1:] = nz + 1
      newpat = zeros (len (pattern) + 1, Int)
      newpat [0] = 1
      newpat [1:] = cumsum (not_equal (pattern, 0)) + 1
      pattern = newpat
      nverts = histogram (pattern) [1:]
      xyzverts = xyz

      # finally, deal with any fcolor function
      if fcolor == None :
         new_results.append ( [nverts, xyzverts, None])
         continue

      # if some polys have been split, need to split clist as well
      if len (list) > len (clist) :
         clist = take (clist, take (cells, list))
      if col == None :
         if nointerp == None :
            if type (fcolor) == FunctionType :
               col = fcolor (m3, clist + cell_offsets [i], lower, upper, fsl,
                  fsu, pattern - 1)
            else :
               col = getc3 (fcolor, m3, clist + cell_offsets [i], lower, upper,
                  fsl, fsu, pattern - 1)
         else :
            if type (fcolor) == FunctionType :
               col = fcolor (m3, clist + cell_offsets [i])
            else :
               col = getc3 (fcolor, m3, clist + cell_offsets [i])
      new_results.append ( [nverts, xyzverts, col])
   # New loop to consolidate the return values
   nv_n = 0
   xyzv_n = 0
   col_n = 0
   for i in range (len (new_results)) :
      nv_n = nv_n + len (new_results [i] [0])
      xyzv_n = xyzv_n + shape (new_results [i] [1]) [0]
      if new_results [i] [2] != None :
         col_n = col_n + len (new_results [i] [2])
   nverts = zeros (nv_n, Int)
   xyzverts = zeros ( (xyzv_n, 3), Float )
   if col_n != 0 :
      col = zeros (col_n, Float )
   else :
      col = None
   nv_n1 = 0
   xyzv_n1 = 0
   col_n1 = 0
   for i in range (len (new_results)) :
      nv_n2 = len (new_results [i] [0])
      xyzv_n2 = shape (new_results [i] [1]) [0]
      nverts [nv_n1:nv_n1 + nv_n2] = new_results [i] [0]
      xyzverts [xyzv_n1:xyzv_n1 + xyzv_n2] = new_results [i] [1]
      if new_results [i] [2] != None :
         col_n2 = len (new_results [i] [2])
         col [col_n1:col_n1 + col_n2] = new_results [i] [2]
         col_n1 = col_n1 + col_n2
      nv_n1 = nv_n1 + nv_n2
      xyzv_n1 = xyzv_n1 + xyzv_n2
   return [nverts, xyzverts, col]

 # The iterator3 function combines three distinct operations:
 # (1) If only the M3 argument is given, return the initial
 #     chunk of the mesh.  The chunk will be no more than
 #     _chunk3_limit cells of the mesh.
 # (2) If only M3 and CHUNK are given, return the next CHUNK,
 #     or [] if there are no more chunks.
 # (3) If M3, CHUNK, and CLIST are all specified, return the
 #     absolute cell index list corresponding to the index list
 #     CLIST of the cells in the CHUNK.
 #     Do not increment the chunk in this case.
 #
 # The form of the CHUNK argument and return value for cases (1)
 # and (2) is not specified, but it must be recognized by the
 # xyz3 and getv3 functions which go along with this iterator3.
 # (For case (3), CLIST and the return value are both ordinary
 #  index lists.)

_Slice3MeshError = "Slice3MeshError"

def slice3mesh (xyz, * args, ** kw) :

   """
   slice3mesh returns a triple [nverts, xyzverts, color]
    nverts is no_cells long and the ith entry tells how many
       vertices the ith cell has.
    xyzverts is sum (nverts) by 3 and gives the vertex
       coordinates of the cells in order.
    color, if present, is len (nverts) long and contains
       a color value for each cell in the mesh.

   There are a number of ways to call slice3mesh:

      slice3mesh (z, color = None)

   z is a two dimensional array of cell function values, assumed
      to be on a uniform mesh nx by ny cells (assuming z is nx by ny)
      nx being the number of cells in the x direction, ny the number
      in the y direction.
   color, if specified, is either an nx by ny array
      of cell-centered values by which the surface is to
      be colored, or an nx +1 by ny + 1 array of vertex-
      centered values, which will be averaged over each
      cell to give cell-centered values.

      slice3mesh (nxny, dxdy, x0y0, z, color = None)

   In this case, slice3mesh accepts the specification for
   a regular 2d mesh: nxny is the number of cells in the
   x direction and the y direction; x0y0 are the initial
   values of x and y; and dxdy are the increments in the
   two directions. z is the height of a surface above
   the xy plane and must be dimensioned nx + 1 by ny + 1.
   color, if specified, is as above.

      slice3mesh (x, y, z, color = None)

   z is as above, an nx by ny array of function values
   on a mesh of the same dimensions. There are two choices
   for x and y: they can both be one-dimensional, dimensioned
   nx and ny respectively, in which case they represent a
   mesh whose edges are parallel to the axes; or else they
   can both be nx by ny, in which case they represent a
   general quadrilateral mesh.
   color, if specified, is as above.
   """

   two_d = 0
   if kw.has_key ("smooth") :
      smooth = kw ["smooth"]
   else :
      smooth = 0
   if len (args) == 0 :
      # Only the z argument is present
      if len (shape (xyz)) != 2 :
         raise _Slice3MeshError, \
            "z must be two dimensional."
      else :
         z = xyz
         ncx = shape (xyz) [0]
         ncy = shape (xyz) [1]
         x = arange (ncx, typecode = Float )
         y = arange (ncy, typecode = Float )
   elif len (args) == 3 :
      # must be the (nxny, dxdy, x0y0, z...) form
      ncx = xyz [0] + 1
      ncy = xyz [1] + 1
      x = arange (ncx, typecode = Float ) * args [0] [0] + args [1] [0]
      y = arange (ncy, typecode = Float ) * args [0] [1] + args [1] [1]
      z = args [2]
      if (ncx, ncy) != shape (z) :
         raise _Slice3MeshError, \
            "The shape of z must match the shape of x and y."
   elif len (args) == 2 :
      # must be the x, y, z format
      x = xyz
      y = args [0]
      z = args [1]
      dims = shape (x)
      if len (dims) == 2 :
         two_d = 1
         if dims != shape (y) or dims != shape (z) :
            raise _Slice3MeshError, \
               "The shapes of x, y, and z must match."
         ncx = dims [0]
         ncy = dims [1]
      elif len (dims) == 1 :
         ncx = dims [0]
         ncy = len (y)
         if (ncx, ncy) != shape (z) :
            raise _Slice3MeshError, \
               "The shape of z must match the shape of x and y."
      else :
         raise _Slice3MeshError, \
            "Unable to decipher arguments to slice3mesh."
   else :
      raise _Slice3MeshError, \
         "Unable to decipher arguments to slice3mesh."

   nverts = ones ( (ncx - 1) *  (ncy - 1), Int) * 4

   ncxx = arange (ncx - 1, typecode = Int) * (ncy)
   ncyy = arange (ncy - 1, typecode = Int)

   if kw.has_key ("color") :
      color = kw ["color"]
   else :
      color = None
   if color != None :
#     col = array (len (nverts), Float )
      if shape (color) == (ncx - 1, ncy - 1) :
         col = color
      elif shape (color) == (ncx, ncy) and smooth == 0 :
         col = ravel (color)
         # Lower left, upper left, upper right, lower right
         col = 0.25 * (take (col, ravel (add.outer ( ncxx, ncyy))) +
            take (col, ravel (add.outer ( ncxx, ncyy + 1))) +
            take (col, ravel (add.outer ( ncxx + ncy, ncyy + 1))) +
            take (col, ravel (add.outer ( ncxx + ncy, ncyy))))
      elif shape (color) == (ncx, ncy) and smooth != 0 :
         # Node-centered colors are wanted (smooth plots)
         col = ravel (color)
         col = ravel (transpose (array ( [
            take (col, ravel (add.outer ( ncxx, ncyy))),
            take (col, ravel (add.outer ( ncxx, ncyy + 1))),
            take (col, ravel (add.outer ( ncxx + ncy, ncyy + 1))),
            take (col, ravel (add.outer ( ncxx + ncy, ncyy)))])))
      else :
         raise _Slice3MeshError, \
            "color must be cell-centered or vertex centered."
   else :
      col = None
   xyzverts = zeros ( (4 * (ncx -1) * (ncy -1), 3), Float )

   if not two_d :
      x1 = multiply.outer (ones (ncy - 1, Float), x [0:ncx - 1])
      x2 = multiply.outer (ones (ncy - 1, Float), x [1:ncx])
      xyzverts [:, 0] = ravel (transpose (array ([x1, x1, x2, x2])))
      del x1, x2
      y1 = multiply.outer (y [0:ncy - 1], ones (ncx - 1))
      y2 = multiply.outer (y [1:ncy], ones (ncx - 1))
      xyzverts [:, 1] = ravel (transpose (array ([y1, y2, y2, y1])))
      del y1, y2
   else :
      newx = ravel (x)
      xyzverts [:, 0] = ravel (transpose (array ( [
         take (newx, ravel (add.outer ( ncxx, ncyy))),
         take (newx, ravel (add.outer ( ncxx, ncyy + 1))),
         take (newx, ravel (add.outer ( ncxx + ncy, ncyy + 1))),
         take (newx, ravel (add.outer ( ncxx + ncy, ncyy)))])))
      newy = ravel (y)
      xyzverts [:, 1] = ravel (transpose (array ( [
         take (newy, ravel (add.outer ( ncxx, ncyy))),
         take (newy, ravel (add.outer ( ncxx, ncyy + 1))),
         take (newy, ravel (add.outer ( ncxx + ncy, ncyy + 1))),
         take (newy, ravel (add.outer ( ncxx + ncy, ncyy)))])))
   newz = ravel (z)
   xyzverts [:, 2] = ravel (transpose (array ( [
      take (newz, ravel (add.outer ( ncxx, ncyy))),
      take (newz, ravel (add.outer ( ncxx, ncyy + 1))),
      take (newz, ravel (add.outer ( ncxx + ncy, ncyy + 1))),
      take (newz, ravel (add.outer ( ncxx + ncy, ncyy)))])))
      
   return [nverts, xyzverts, col]
   
def iterator3 (m3 , chunk = None, clist = None) :

   """
   iterator3 (m3)
   iterator3 (m3, chunk, clist)
   iterator3_rect (m3)
   iterator3_rect (m3, chunk, clist)
   iterator3_irreg (m3)
   iterator3_irreg (m3, chunk, clist)

   The iterator3 functions combine three distinct operations:
   (1) If only the M3 argument is given, return the initial
       chunk of the mesh.  The chunk will be no more than
       chunk3_limit cells of the mesh.
   (2) If only M3 and CHUNK are given, return the next CHUNK,
       or None if there are no more chunks.
   (3) If M3, CHUNK, and CLIST are all specified, return the
       absolute cell index list corresponding to the index list
       CLIST of the cells in the CHUNK.
       Do not increment the chunk in this case.

   The form of the CHUNK argument and return value for cases (1)
   and (2) is not specified, but it must be recognized by the
   xyz3 and getv3 functions which go along with this iterator3.
   (For case (3), CLIST and the return value are both ordinary
   index lists.)
   In the irregular case, it is guaranteed that the returned chunk
   consists of only one type of cell (tetrahedra, hexahedra,
   pyramids, or prisms).
   """

   return m3 [0] [3] (m3, chunk, clist)

# biggest temporary is 3 doubles times this,
# perhaps 4 or 5 doubles times this is most at one time
_chunk3_limit = 10000

def iterator3_rect (m3, chunk, clist) :

   """
   iterator3 (m3)
   iterator3 (m3, chunk, clist)
   iterator3_rect (m3)
   iterator3_rect (m3, chunk, clist)
   iterator3_irreg (m3)
   iterator3_irreg (m3, chunk, clist)

   The iterator3 functions combine three distinct operations:
   (1) If only the M3 argument is given, return the initial
       chunk of the mesh.  The chunk will be no more than
       chunk3_limit cells of the mesh.
   (2) If only M3 and CHUNK are given, return the next CHUNK,
       or None if there are no more chunks.
   (3) If M3, CHUNK, and CLIST are all specified, return the
       absolute cell index list corresponding to the index list
       CLIST of the cells in the CHUNK.
       Do not increment the chunk in this case.

   The form of the CHUNK argument and return value for cases (1)
   and (2) is not specified, but it must be recognized by the
   xyz3 and getv3 functions which go along with this iterator3.
   (For case (3), CLIST and the return value are both ordinary
   index lists.)
   In the irregular case, it is guaranteed that the returned chunk
   consists of only one type of cell (tetrahedra, hexahedra,
   pyramids, or prisms).
   """

#  Note: if you look at the yorick version of this routine, you
#  will see that the significance of the subscripts is reversed.
#  This is because we do things in row-major order.

   global _chunk3_limit
   
   if chunk == None :
      dims = m3 [1] [0]      # [ni,nj,nk] cell dimensions
      [ni, nj, nk] = [dims [0], dims [1], dims [2]]
      njnk = nj * nk
      if _chunk3_limit <= nk :
         # stuck with 1D chunks
         ck = (nk - 1) / _chunk3_limit + 1
         cj = ci = 0
      elif _chunk3_limit <= njnk :
         # 2D chunks
         ci = ck = 0
         cj = (njnk - 1) / _chunk3_limit + 1
      else :
         # 3D chunks
         cj = ck = 0
         ci = (njnk * ni - 1) / _chunk3_limit + 1
      chunk = array ( [[ci == 0, cj == 0, ck == 0],
                       [not ci, nj * (ci != 0) + (ck != 0),
                        nk * ( (cj + ci) != 0)],
                       [ci, cj, ck], [ni, nj, nk]])
   else :
      ni = chunk [3,0]
      nj = chunk [3,1]
      nk = chunk [3,2]
      njnk = nj * nk
      offsets = array ( [njnk, nj, 1], Int)
      if clist != None :
         # add offset for this chunk to clist and return
         return sum (offsets * ( chunk [0] - 1)) + clist

   # increment to next chunk
   xi = chunk [1, 0]
   xj = chunk [1, 1]
   xk = chunk [1, 2]

   np = chunk [2, 2]
   if (np) :
      # 1D chunks
      if xk == nk :
         if xj == nj :
            if xi == ni : return None
            xi = xi + 1
            xj = 1;
         else :
            xj = xj + 1
         xk = 0
      ck = xk + 1
      step = ck / np 
      frst = ck % np     # first frst steps are step+1
      if (xk < (step + 1) * frst) : step = step + 1
      xk = xk + step
      chunk [0] = array ( [xi, xj, ck])
      chunk [1] = array ( [xi, xj, xk])
   else :
      np = chunk [2, 1]
      if (np) :
         if (xj == nj) :
            if (xi == ni) : return None
            xi = xi + 1
            xj = 0
         cj = xj + 1
         step = nj / np
         frst = nj % np    # first frst steps are step+1
         if (xj < (step + 1) * frst) : step = step + 1
         xj = xj + step
         chunk [0, 0:2] = array ( [xi, cj])
         chunk [1, 0:2] = array ( [xi, xj])
      else :
         if xi == ni : return None
         ci = xi + 1
         np = chunk [2, 0]
         step = ni / np
         frst = ni % np    # first frst steps are step+1
         if (xi < (step + 1) * frst) : step = step + 1
         xi = xi + step
         chunk [0, 0] = ci
         chunk [1, 0] = xi
   return chunk

def iterator3_irreg (m3, chunk, clist) :

   """
   iterator3 (m3)
   iterator3 (m3, chunk, clist)
   iterator3_rect (m3)
   iterator3_rect (m3, chunk, clist)
   iterator3_irreg (m3)
   iterator3_irreg (m3, chunk, clist)

   The iterator3 functions combine three distinct operations:
   (1) If only the M3 argument is given, return the initial
       chunk of the mesh.  The chunk will be no more than
       chunk3_limit cells of the mesh.
   (2) If only M3 and CHUNK are given, return the next CHUNK,
       or None if there are no more chunks.
   (3) If M3, CHUNK, and CLIST are all specified, return the
       absolute cell index list corresponding to the index list
       CLIST of the cells in the CHUNK.
       Do not increment the chunk in this case.

   The form of the CHUNK argument and return value for cases (1)
   and (2) is not specified, but it must be recognized by the
   xyz3 and getv3 functions which go along with this iterator3.
   (For case (3), CLIST and the return value are both ordinary
   index lists.)
   In the irregular case, it is guaranteed that the returned chunk
   consists of only one type of cell (tetrahedra, hexahedra,
   pyramids, or prisms).
 
   iterator3_irreg Does the same thing as iterator3_rect only for an 
   irregular rectangular mesh. It simply splits a large mesh into smaller
   parts. Whether this is necessary I am not sure.
   Certainly it makes it easier in the irregular case to handle
   the four different types of cells separately.
   if clist is present, in the irregular case it is already
   the list of absolute cell indices, so it is simply returned.
   This and other routines to do with irregular meshes return a
   chunk which is a 2-list. The first item delimits the chunk;
   the second gives a list of corresponding cell numbers.
   """

   global _chunk3_limit

   if clist != None:
      return clist

   dims = m3 [1] [0]     # ncells by _no_verts array of subscripts
                         # (or a list of from one to four of same)

   if type (dims) != ListType :
      if chunk == None:     # get the first chunk
         return [ [0, min (shape (dims) [0], _chunk3_limit)],
                  arange (0, min (shape (dims) [0], _chunk3_limit),
                  typecode = Int)]
      else :                # iterate to next chunk
         start = chunk [0] [1]
         if start >= shape(dims) [0] :
            return None
         else :
            return [ [start, min (shape (dims) [0], start + _chunk3_limit)],
                     arange (start, min (shape (dims) [0],
                                               start + _chunk3_limit),
                     typecode = Int)]
   else :
      totals = m3 [1] [3] # cumulative totals of numbers of cells
      if chunk == None :
         return [ [0, min (totals [0], _chunk3_limit)],
                  arange (0, min (totals [0], _chunk3_limit),
                  typecode = Int)]
      else :                # iterate to next chunk
         start = chunk [0] [1]
         if start >= totals [-1] :
            return None
         else :
            for i in range (len (totals)) :
               if start < totals [i] :
                  break
            return [ [start, min (totals [i], start + _chunk3_limit)],
                     arange (start,
                        min (totals [i], start + _chunk3_limit), 
                        typecode = Int)]


def getv3 (i, m3, chunk) :

   """
   getv3(i, m3, chunk)

     return vertex values of the Ith function attached to 3D mesh M3
     for cells in the specified CHUNK.  The CHUNK may be a list of
     cell indices, in which case getv3 returns a 2x2x2x(dimsof(CHUNK))
     list of vertex coordinates.  CHUNK may also be a mesh-specific data
     structure used in the slice3 routine, in which case getv3 may
     return a (ni)x(nj)x(nk) array of vertex values.  For meshes which
     are logically rectangular or consist of several rectangular
     patches, this is up to 8 times less data, with a concomitant
     performance advantage.  Use getv3 when writing slicing functions
     for slice3.
   """

   return m3 [0] [1] (i, m3, chunk)

_Getv3Error = "Getv3Error"

def getv3_rect (i, m3, chunk) :

   """
   getv3_rect(i, m3, chunk) does the job for a regular rectangular
     mesh.
   """

   fi = m3 [2]
   i = i - 1
   if i < 0 or is_scalar (fi) or i >= len (fi) :
      raise _Getv3Error, "no such mesh function as F" + `i`
   dims = m3 [1] [0]
   if dims == shape (fi [i]) :
      raise _Getv3Error, "mesh function F" + `i` + " is not vertex-centered"
   if len (shape (chunk)) != 1 :
      c = chunk
      # The difference here is that our arrays are 0-based, while
      # yorick's are 1-based; and the last element in a range is not
      # included in the result array.
      return fi [i] [c [0, 0] - 1:1 + c [1, 0], c [0, 1] - 1:1 + c [1, 1] ,
                     c [0, 2] - 1:1 + c [1, 2]]
   else :
      # Need to create an array of fi values the same size and shape
      # as what to_corners3 returns.
      # To avoid exceedingly arcane calculations attempting to
      # go backwards to a cell list, this branch returns the list
      # [<function values>, chunk]
      # Then it is trivial for slice3 to find a list of cell
      # numbers in which fi changes sign.
      indices = to_corners3 (chunk, dims [0] + 1, dims [1] + 1)
      no_cells = shape (indices) [0]
      indices = ravel (indices)
      retval = reshape (take (ravel (fi [i]), indices), (no_cells, 2, 2, 2))

      return [retval, chunk]

def getv3_irreg (i, m3, chunk) :

   """
     for an irregular mesh, returns a 3-list whose elements are:
     (1) the function values for the ith function on the vertices of the
     given chunk. (The function values must have the same dimension
     as the coordinates; there is no attempt to convert zone-centered
     values to vertex-centered values.)
     (2) an array of relative cell numbers within the list of cells
     of this type.
     (3) a number that can be added to these relative numbers to gives
     the absolute cell numbers for correct access to their coordinates
     and function values.
   """
   
   fi = m3 [2]
   i = i - 1
   if i < 0 or is_scalar (fi) or i >= len (fi) :
      raise _Getv3Error, "no such function as F" + `i`
   # len (fi [i]) and the second dimension of m3 [1] [1] (xyz) should
   # be the same, i. e., there is a value associated with each coordinate.
   if len (fi [i]) != len (m3 [1] [1] [0]) :
      raise _Getv3Error, "mesh function F" + `i` + " is not vertex-centered."

   verts = m3 [1] [0]
   oldstart = chunk [0] [0]
   oldfin = chunk [0] [1]
   no_cells = oldfin - oldstart

   if type (verts) != ListType : # Only one kind of cell in mesh
      indices = ravel (verts [oldstart:oldfin])
   else : # A list of possibly more than one kind
      sizes = m3 [1] [2]
      totals = m3 [1] [3]
      for j in range (len (totals)) :
         if oldfin <= totals [j] :
            break
      verts = verts [j]
      if j > 0 :
         start = oldstart - totals [j - 1]
         fin = oldfin - totals [j - 1]
      else :
         start = oldstart 
         fin = oldfin
      indices = ravel (verts [start:fin])

   tc = shape (verts) [1]
   # ZCM 2/4/97 the array of cell numbers must be relative
   if tc == 8 : # hex cells
      return [ reshape (take (fi [i], indices), (no_cells, 2, 2, 2)),
              arange (0, no_cells, typecode = Int), oldstart]
   elif tc == 6 : # pyramids
      return [ reshape (take (fi [i], indices), (no_cells, 3, 2)),
              arange (0, no_cells, typecode = Int), oldstart]
   else : # tetrahedron or pyramid
      return [ reshape (take (fi [i], indices), (no_cells, tc)),
              arange (0, no_cells, typecode = Int), oldstart]

_Getc3Error = "Getc3Error"

def getc3 (i, m3, chunk, *args) :

   """
   getc3(i, m3, chunk)
         or getc3(i, m3, clist, l, u, fsl, fsu, cells)

     return cell values of the Ith function attached to 3D mesh M3
     for cells in the specified CHUNK.  The CHUNK may be a list of
     cell indices, in which case getc3 returns a (dimsof(CHUNK))
     list of vertex coordinates.  CHUNK may also be a mesh-specific data
     structure used in the slice3 routine, in which case getc3 may
     return a (ni)x(nj)x(nk) array of vertex values.  There is no
     savings in the amount of data for such a CHUNK, but the gather
     operation is cheaper than a general list of cell indices.
     Use getc3 when writing colorng functions for slice3.

     If CHUNK is a CLIST, the additional arguments L, U, FSL, and FSU
     are vertex index lists which override the CLIST if the Ith attached
     function is defined on mesh vertices.  L and U are index lists into
     the (dimsof(CLIST))x2x2x2 vertex value array, say vva, and FSL
     and FSU are corresponding interpolation coefficients; the zone
     centered value is computed as a weighted average of involving these
     coefficients.  The CELLS argument is required by histogram to do
     the averaging.  See the source code for details.
     By default, this conversion (if necessary) is done by averaging
     the eight vertex-centered values.
    """

   if len (args) == 0 :
      l = None
      u = None
      fsl = None
      fsu = None
      cells = None
   elif len (args) == 5 :
      l = args [0]
      u = args [1]
      fsl = args [2]
      fsu = args [3]
      cells = args [4]
   else :
      raise _Getc3Error, "getc3 requires either three or eight parameters."

   return m3 [0] [2] (i, m3, chunk, l, u, fsl, fsu, cells)

def getc3_rect (i, m3, chunk, l, u, fsl, fsu, cells) :

   """
   getc3_rect (i, m3, chunk, l, u, fsl, fsu, cells) does the job
     for a regular rectangular mesh.
   """

   fi = m3 [2]
   m3 = m3 [1]
   if ( i < 1 or i > len (fi)) :
      raise _Getc3Error, "no such mesh function as F" + `i - 1`
   dims = m3 [0]
   if shape (fi [i - 1]) == dims :
      # it is a cell-centered quantity
      if len (shape (chunk)) != 1 :
         c = chunk
         # The difference here is that our arrays are 0-based, while
         # yorick's are 1-based; and the last element in a range is not
         # included in the result array.
         return fi [i - 1] [c [0, 0] - 1:1 + c [1, 0],
                            c [0, 1] - 1:1 + c [1, 1] ,
                            c [0, 2] - 1:1 + c [1, 2]]
      else :
         [k, l. m] = dims
         return reshape (take (ravel (fi [i - 1]), chunk),
            (len (chunk), k, l, m))
   else :
      # it is vertex-centered, so we take averages to get cell quantity
      if len (shape (chunk)) != 1 :
         c = chunk
         # The difference here is that our arrays are 0-based, while
         # yorick's are 1-based; and the last element in a range is not
         # included in the result array.
         return zcen_ (zcen_( zcen_ (
               (fi [i - 1] [c [0, 0] - 1:1 + c [1, 0],
                            c [0, 1] - 1:1 + c [1, 1] ,
                            c [0, 2] - 1:1 + c [1, 2]]), 0), 1), 2)
      else :
         indices = to_corners3 (chunk, dims [1] + 1,  dims [2] + 1)
         no_cells = shape (indices) [0]
         indices = ravel (indices)
         corners = take (ravel (fi [i - 1]), indices)
         if l == None :
            return 0.125 * sum (transpose (reshape (corners, (no_cells, 8))))
         else :
            # interpolate corner values to get edge values
            corners = (take (corners, l) * fsu -
               take (corners, u) * fsl) / (fsu -fsl)
            # average edge values (vertex values of polys) on each poly
            return histogram (cells, corners) / histogram (cells)

def getc3_irreg (i, m3, chunk, l, u, fsl, fsu, cells) :

   """
      Same thing as getc3_rect, i. e., returns the same type of
      data structure, but from an irregular mesh.
      m3 [1] is a 2-list; m3[1] [0] is an array whose ith element
         is an array of coordinate indices for the ith cell,
         or a list of up to four such arrays.
         m3 [1] [1] is the 3 by nverts array of coordinates.
      m3 [2] is a list of arrays of vertex-centered or cell-centered
         data.
      chunk may be a list, in which case chunk [0] is a 2-sequence
      representing a range of cell indices; or it may be a one-dimensional
      array, in which case it is a nonconsecutive set of cell indices.
      It is guaranteed that all cells indexed by the chunk are the
      same type.
   """

   fi = m3 [2]
   if i < 1 or i > len (fi) :
      raise _Getc3Error, "no such mesh function as F" + `i - 1`
   verts = m3 [1] [0]
   if type (verts) == ListType :
      sizes = m3 [1] [2]
      totals = m3 [1] [3]
   if type (verts) == ListType and totals [-1] == len (fi [i - 1]) or \
      type (verts) != ListType and shape (verts) [0] == len (fi [i - 1]) :
      # cell-centered case
      if type (chunk) == ListType :
         return fi [i - 1] [chunk [0] [0]:chunk [0] [1]]
      elif type (chunk) == ArrayType and len (shape (chunk)) == 1 :
         return take (fi [i - 1], chunk)
      else :
         raise _Getc3Error, "chunk argument is incomprehensible."

   if len (fi [i - 1]) != shape (m3 [1] [1]) [1] :
      raise _Getc3Error, "F" + `i - 1` + " has the wrong size to be " \
         "either zone-centered or node-centered."
   # vertex-centered case
   # First we need to pick up the vertex subscripts, which are
   # also the fi [i - 1] subscripts.
   if type (verts) != ListType :
      if type (chunk) == ListType :
         indices = verts [chunk [0] [0]:chunk [0] [1]]
      elif type (chunk) == ArrayType and len (shape (chunk)) == 1 :
         indices = take (verts, chunk)
      else :
         raise _Getc3Error, "chunk argument is incomprehensible."
   else :
      # We have a list of vertex subscripts, each for a different
      # type of cell; need to extract the correct list:
      if type (chunk) == ListType :
         start = chunk [0] [0]
         fin = chunk [0] [1]
         for j in range (len (totals)) :
            if fin <= totals [j] :
               break
         verts = verts [j]
         if j > 0 :
            start = start - totals [j - 1]
            fin = fin - totals [j - 1]
         indices = verts [start:fin]
      elif type (chunk) == ArrayType and len (shape (chunk)) == 1 :
         for j in range (len (totals)) :
            if chunk [-1] <= totals [j] :
               break
         verts = verts [j]
         ch = chunk
         if j > 0 :
            ch = chunk - totals [j - 1]
         indices = take (verts, ch)
      else :
         raise _Getc3Error, "chunk argument is incomprehensible."

   shp = shape (indices)
   no_cells = shp [0]
   indices = ravel (indices)
   corners = take (fi [i - 1], indices)
   if l == None :
      return (1. / shp [1]) * transpose ((sum (transpose (reshape (corners,
         (no_cells, shp [1]))) [0:shp [1]])))
   else :
      # interpolate corner values to get edge values
      corners = (take (corners, l) * fsu -
         take (corners, u) * fsl) / (fsu -fsl)
      # average edge values (vertex values of polys) on each poly
      return histogram (cells, corners) / histogram (cells)

_no_verts = array ( [4, 5, 6, 8])
_no_edges = array ( [6, 8, 9, 12])

# Lower and upper vertex subscripts for each edge
_lower_vert4 = array ( [0, 0, 0, 1, 2, 3], Int)
_lower_vert5 = array ( [0, 0, 0, 0, 1, 2, 3, 4], Int)
_lower_vert6 = array ( [0, 1, 0, 1, 2, 3, 0, 2, 4], Int)
_lower_vert8 = array ( [0, 1, 2, 3, 0, 1, 4, 5, 0, 2, 4, 6], Int)
_lower_vert = [_lower_vert4, _lower_vert5, _lower_vert6, _lower_vert8]
_upper_vert4 = array ( [1, 2, 3, 2, 3, 1], Int)
_upper_vert5 = array ( [1, 2, 3, 4, 2, 3, 4, 1], Int)
_upper_vert6 = array ( [4, 5, 2, 3, 4, 5, 1, 3, 5], Int)
_upper_vert8 = array ( [4, 5, 6, 7, 2, 3, 6, 7, 1, 3, 5, 7], Int)
_upper_vert = [_upper_vert4, _upper_vert5, _upper_vert6, _upper_vert8]

_node_edges8_s = array ( [ [1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0],
                        [0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0],
                        [0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0],
                        [0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0],
                        [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0],
                        [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
                        [0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1],
                        [0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1]], Int)
_node_edges8 = array ( [ [0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1],
                        [0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1],
                        [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
                        [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0],
                        [0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0],
                        [0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0],
                        [0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0],
                        [1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0]], Int)
_node_edges6_s = array ( [ [1, 0, 1, 0, 0, 0, 1, 0, 0],
                        [0, 1, 0, 1, 0, 0, 1, 0, 0],
                        [0, 0, 1, 0, 1, 0, 0, 1, 0],
                        [0, 0, 0, 1, 0, 1, 0, 1, 0],
                        [1, 0, 0, 0, 1, 0, 0, 0, 1],
                        [0, 1, 0, 0, 0, 1, 0, 0, 1]], Int)
_node_edges6 = array ( [ [0, 1, 0, 0, 0, 1, 0, 0, 1],
                        [1, 0, 0, 0, 1, 0, 0, 0, 1],
                        [0, 0, 0, 1, 0, 1, 0, 1, 0],
                        [0, 0, 1, 0, 1, 0, 0, 1, 0],
                        [0, 1, 0, 1, 0, 0, 1, 0, 0],
                        [1, 0, 1, 0, 0, 0, 1, 0, 0]], Int)
_node_edges4_s = array ( [ [1, 1, 1, 0, 0, 0],
                        [1, 0, 0, 1, 0, 1],
                        [0, 1, 0, 1, 1, 0],
                        [0, 0, 1, 0, 1, 1]], Int)
_node_edges4 = array ( [ [0, 0, 1, 0, 1, 1],
                        [0, 1, 0, 1, 1, 0],
                        [1, 0, 0, 1, 0, 1],
                        [1, 1, 1, 0, 0, 0]], Int)
_node_edges5_s = array ( [ [1, 1, 1, 1, 0, 0, 0, 0],
                        [1, 0, 0, 0, 1, 0, 0, 1],
                        [0, 1, 0, 0, 1, 1, 0, 0],
                        [0, 0, 1, 0, 0, 1, 1, 0],
                        [0, 0, 0, 1, 0, 0, 1, 1]], Int)
_node_edges5 = array ( [ [0, 0, 0, 1, 0, 0, 1, 1],
                        [0, 0, 1, 0, 0, 1, 1, 0],
                        [0, 1, 0, 0, 1, 1, 0, 0],
                        [1, 0, 0, 0, 1, 0, 0, 1],
                        [1, 1, 1, 1, 0, 0, 0, 0]], Int)

_node_edges = [_node_edges4_s, _node_edges5_s, _node_edges6_s, _node_edges8_s]
_node_edges3 = [_node_edges4, _node_edges5, _node_edges6, _node_edges8]

def _construct3 (itype) :
   global _node_edges
   global _no_verts
   global _no_edges
   i = arange (1, 2**_no_verts [itype] - 1, typecode = Int)
   if itype == 0 :
      below = transpose (not_equal (array ( [bitwise_and (i, 8),
                                             bitwise_and (i, 4),
                                             bitwise_and (i, 2),
                                             bitwise_and (i, 1)]), 0))
   elif itype == 1 :
      below = transpose (not_equal (array ( [bitwise_and (i, 16),
                                             bitwise_and (i, 8),
                                             bitwise_and (i, 4),
                                             bitwise_and (i, 2),
                                             bitwise_and (i, 1)]), 0))
   elif itype == 2 :
      below = transpose (not_equal (array ( [bitwise_and (i, 32),
                                             bitwise_and (i, 16),
                                             bitwise_and (i, 8),
                                             bitwise_and (i, 4),
                                             bitwise_and (i, 2),
                                             bitwise_and (i, 1)]), 0))
   elif itype == 3 :
      below = transpose (not_equal (array ( [bitwise_and (i, 128),
                                             bitwise_and (i, 64),
                                             bitwise_and (i, 32),
                                             bitwise_and (i, 16),
                                             bitwise_and (i, 8),
                                             bitwise_and (i, 4),
                                             bitwise_and (i, 2),
                                             bitwise_and (i, 1)]), 0))
   # For some reason the node edges for a cell need to be in different order
   # here than in slice3 to get the correct results. Hence _node_edges3.
   mask = find_mask (below, _node_edges3 [itype])

   return construct3 (mask, itype)

# ------------------------------------------------------------------------

_poly_permutations4 = _construct3 (0)
_poly_permutations5 = _construct3 (1)
_poly_permutations6 = _construct3 (2)
_poly_permutations8 = _construct3 (3)

_poly_permutations = [_poly_permutations4, _poly_permutations5, 
                     _poly_permutations6, _poly_permutations8]

_ContourError = "ContourError"

# ------------------------------------------------------------------------

def plzcont (nverts, xyzverts, contours = 8, scale = "lin", clear = 1,
   edges = 0, color = None, cmin = None, cmax = None, 
   zaxis_min = None, zaxis_max = None, split = 0) :

   """
   plzcont (nverts, xyzverts, contours = 8, scale = "lin", clear = 1,
   edges = 0, color = None, cmin = None, cmax = None, split = 0
   zaxis_min = None, zaxis_max = None, )

     Plot filled z contours on the specified surface. NVERTS and
     XYZVERTS arrays specify the polygons for the surface being
     drawn. CONTOURS can be one of the following:
        N, an integer: Plot N contours (therefore, N+1 colored
        components of the surface)
        CVALS, a vector of floats: draw the contours at the
        specified levels.
     SCALE can be "lin", "log", or "normal" specifying the
     contour scale. (Only applicable if contours = N, of course).
     If CLEAR = 1, clear the display list first.
     If EDGES = 1, plot the edges.
     The algorithm is to apply slice2x repeatedly to the surface.
     If color == None, then bytscl the palette into N + 1 colors
     and send each of the slices to pl3tree with the appropriate color.
     If color == "bg", will plot only the edges.
     If CMIN is given, use it instead of the minimum z actually
     being plotted in the computation of contour levels. If CMAX is given,
     use it instead of the maximum z actually being plotted in the
     computation of contour levels. This is done so that a component
     of a larger graph will have the same colors at the same levels
     as every other component, rather than its levels being based
     on its own max and min, which may lie inside those of the
     rest of the graph.
     ZAXIS_MIN and ZAXIS_MAX represent axis limits on z as expressed
     by the user. If present, ZAXIS_MIN will inhibit plotting of all
     lesser z values, and ZAXIS_MAX will inhibit the plotting of all
     greater z values.
   """

   # 1. Get contour colors
   if type (contours) == IntType :
      n = contours
      if cmin != None :
         vcmin = cmin
         minz = min (xyzverts [:, 2])
      else :
         vcmin = min (xyzverts [:, 2])
         minz = vcmin
      if cmax != None :
         vcmax = cmax
         maxz = max (xyzverts [:, 2])
      else :
         vcmax = max (xyzverts [:, 2])
         maxz = vcmax
      if scale == "lin" :
          vc = vcmin + arange (1, n + 1, typecode = Float) * \
             (vcmax - vcmin) / (n + 1)
      elif scale == "log" :
          vc = vcmin + exp (arange (1, n + 1, typecode = Float) * \
             log (vcmax - vcmin) / (n + 1))
      elif scale == "normal" :
          zlin = xyzverts [:, 2]
          lzlin = len (zlin)
          zbar = add.reduce (zlin) / lzlin
          zs = sqrt ( (add.reduce (zlin ** 2) - lzlin * zbar ** 2) /
              (lzlin - 1))
          z1 = zbar - 2. * zs
          z2 = zbar + 2. * zs
          diff = (z2 - z1) / (n - 1)
          vc = z1 + arange (n) * diff
      else :
          raise _ContourError, "Incomprehensible scale parameter."
   elif type (contours) == ArrayType and contours.typecode () == Float :
      n = len (contours)
      vc = sort (contours)
   else :
      raise _ContourError, "Incorrect contour specification."
   if split == 0 :
      colors = (arange (n + 1, typecode = Float) * (199. / n)).astype ('b')
   else :
      colors = (arange (n + 1, typecode = Float) * (99. / n)).astype ('b')
   # 2. Loop through slice2x calls
   nv = array (nverts, copy = 1)
   xyzv = array (xyzverts, copy = 1)
   if clear == 1 :
      clear3 ( ) # Clear any previous plot or we're in trouble
   # find imin--contours below this number need not be computed,
   # and imax--contours at this level and above need not be computed.
   imin = imax = 0
   for i in range (n) :
      if vc [i] <= minz :
         imin = i + 1
      if vc [i] >= maxz :
         imax = i
         break
      if i == n - 1 :
         imax = n
   # now make sure that the minimum and maximum contour levels computed
   # are not outside the axis limits.
   if zaxis_min != None and zaxis_min > vc [imin] :
      for i in range (imin, imax) :
         if i + 1 < imax and zaxis_min > vc [i + 1] :
            imin = i + 1
         else :
            break
      vc [imin] = zaxis_min
   if zaxis_max != None and zaxis_max < vc [imax - 1] :
      for i in range (imax - imin) :
         if imax - 2 >= imin and zaxis_max < vc [imax - 2] :
            imax = imax - 1
         else :
            break
      vc [imax - 1] = zaxis_max
   for i in range (imin, imax) :
      [nv, xyzv, d1, nvb, xyzvb, d2] = \
         slice2x (array ( [0., 0., 1., vc [i]], Float) , nv, xyzv, None)
      if i == imin and zaxis_min != None and zaxis_min == vc [i]:
         # Don't send the "back" surface if it's below zaxis_min.
         continue
      else:
         if color == None :
            pl3tree (nvb, xyzvb, (ones (len (nvb)) * colors [i]).astype ('b'),
               split = 0, edges = edges)
         else :
            # N. B. Force edges to be on, otherwise the graph is empty.
            pl3tree (nvb, xyzvb, "bg", split = 0, edges = 1)
   if zaxis_max == None or vc [imax - 1] < zaxis_max:
      # send "front" surface if it's not beyond zaxis_max
      if color == None :
         pl3tree (nv, xyzv, (ones (len (nv)) * colors [i]).astype ('b'),
            split = 0, edges = edges)
      else :
         pl3tree (nv, xyzv, "bg", split = 0, edges = 1)

def pl4cont (nverts, xyzverts, values, contours = 8, scale = "lin", clear = 1,
   edges = 0, color = None, cmin = None, cmax = None,
   caxis_min = None, caxis_max = None, split = 0) :

   """
   pl4cont (nverts, xyzverts, values, contours = 8, scale = "lin", clear = 1,
   edges = 0, color = None, cmin = None, cmax = None,
   caxis_min = None, caxis_max = None, split = 0)

     Plot filled z contours on the specified surface. VALUES is
     a node-centered array the same length as SUM (NVERTS) whose
     contours will be drawn. NVERTS and
     XYZVERTS arrays specify the polygons for the surface being
     drawn. CONTOURS can be one of the following:
        N, an integer: Plot N contours (therefore, N+1 colored
        components of the surface)
        CVALS, a vector of floats: draw the contours at the
        specified levels.
     SCALE can be "lin", "log", or "normal" specifying the
     contour scale. (Only applicable if contours = N, of course).
     If CLEAR == 1, clear the display list first.
     If EDGES == 1, plot the edges.
     The algorithm is to apply slice2x repeatedly to the surface.
     If color == None, then bytscl the palette into N + 1 colors
     and send each of the slices to pl3tree with the appropriate color.
     If color == "bg", will plot only the edges.
     If CMIN is given, use it instead of the minimum c actually
     being plotted in the computation of contour levels. If CMAX is given,
     use it instead of the maximum c actually being plotted in the
     computation of contour levels. This is done so that a component
     of a larger graph will have the same colors at the same levels
     as every other component, rather than its levels being based
     on its own max and min, which may lie inside those of the
     rest of the graph.
     CAXIS_MIN and CAXIS_MAX represent axis limits on c as expressed
     by the user. If present, CAXIS_MIN will inhibit plotting of all
     lesser c values, and CAXIS_MAX will inhibit the plotting of all
     greater c values.
   """

   # 1. Get contour colors
   if type (contours) == IntType :
      n = contours
      if cmin != None :
          vcmin = cmin
          minz = min (values)
      else :
          vcmin = min (values)
          minz = vcmin
      if cmax != None :
          vcmax = cmax
          maxz = max (values)
      else :
          vcmax = max (values)
          maxz = vcmax
      if scale == "lin" :
          vc = vcmin + arange (1, n + 1, \
             typecode = Float) * \
             (vcmax - vcmin) / (n + 1)
      elif scale == "log" :
          vc = vcmin + exp (arange (1, n + 1, \
             typecode = Float) * \
             log (vcmax - vcmin) / (n + 1))
      elif scale == "normal" :
          zbar = add.reduce (values) / lzlin
          zs = sqrt ( (add.reduce (values ** 2) - lzlin * zbar ** 2) /
              (lzlin - 1))
          z1 = zbar - 2. * zs
          z2 = zbar + 2. * zs
          diff = (z2 - z1) / (n - 1)
          vc = z1 + arange (n) * diff
      else :
          raise _ContourError, "Incomprehensible scale parameter."
   elif type (contours) == ArrayType and contours.typecode () == Float :
      n = len (contours)
      vc = sort (contours)
   else :
      raise _ContourError, "Incorrect contour specification."
   if split == 0 :
      colors = (arange (n + 1, typecode = Float) * (199. / n)).astype ('b')
   else :
      colors = (arange (n + 1, typecode = Float) * (99. / n)).astype ('b')
   # 2. Loop through slice2x calls
   nv = array (nverts, copy = 1)
   xyzv = array (xyzverts, copy = 1)
   vals = array (values, copy = 1)
   if clear == 1 :
      clear3 ( ) # Clear any previous plot or we're in trouble
   # find imin--contours below this number need not be computed,
   # and imax--contours at this level and above need not be computed.
   imin = imax = 0
   for i in range (n) :
      if vc [i] <= minz :
         imin = i + 1
      if vc [i] >= maxz :
         imax = i
         break
      if i == n - 1 :
         imax = n
   # now make sure that the minimum and maximum contour levels computed
   # are not outside the axis limits.
   if caxis_min != None and caxis_min > vc [imin] :
      for i in range (imin, imax) :
         if i + 1 < imax and caxis_min > vc [i + 1] :
            imin = i + 1
         else :
            break
      vc [imin] = caxis_min
   if caxis_max != None and caxis_max < vc [imax - 1] :
      for i in range (imax - imin) :
         if imax - 2 >= imin and caxis_max < vc [imax - 2] :
            imax = imax - 1
         else :
            break
      vc [imax - 1] = caxis_max
   for i in range (n) :
      if vc [i] <= minz :
         continue
      if vc [i] >= maxz :
         break
      [nv, xyzv, vals, nvb, xyzvb, d2] = \
         slice2x (vc [i], nv, xyzv, vals)
      if i == imin and caxis_min != None and caxis_min == vc [i]:
         # Don't send the "back" surface if it's below caxis_min.
         continue
      else:
         if color == None :
            pl3tree (nvb, xyzvb, (ones (len (nvb)) * colors [i]).astype ('b'),
               split = 0, edges = edges)
         else :
            # N. B. Force edges to be on, otherwise the graph is empty.
            pl3tree (nvb, xyzvb, "bg", split = 0, edges = 1)
   if caxis_max == None or vc [imax - 1] < caxis_max:
      # send "front" surface if it's not beyond caxis_max
      if color == None :
         pl3tree (nv, xyzv, (ones (len (nv)) * colors [i]).astype ('b'),
            split = 0, edges = edges)
      else :
         pl3tree (nv, xyzv, "bg", split = 0, edges = 1)

def slice2x (plane, nverts, xyzverts, values = None) :

   """
   slice2x (plane, nverts, xyzverts, values)

     Slice a polygon list, retaining only those polygons or
     parts of polygons on the positive side of PLANE, that is,
     the side where xyz(+)*PLANE(+:1:3)-PLANE(4) > 0.0.
     The NVERTS, VALUES, and XYZVERTS arrays have the meanings of
     the return values from the slice3 function.
     Python returns a sextuple
     [nverts, xyzverts, values, nvertb, xyzvertb, valueb]
     with None in the place of missing or None input arguments.

   slice2_precision= precision
     Controls how slice2 (or slice2x) handles points very close to
     the slicing plane.  PRECISION should be a positive number or zero.
     Zero PRECISION means to clip exactly to the plane, with points
     exactly on the plane acting as if they were slightly on the side
     the normal points toward.  Positive PRECISION means that edges
     are clipped to parallel planes a distance PRECISION on either
     side of the given plane.  (Polygons lying entirely between these
     planes are completely discarded.)

     Default value is 0.0.

   """

#    Note (ZCM 2/24/97) Reomved _slice2x as a global and added
#    it as a final argument to slice2.

   retval = slice2 (plane, nverts, xyzverts, values, 1)
   retval = retval + [None] * (6 - len (retval))
   return retval


_Pl3surfError = "Pl3surfError"

def pl3surf(nverts, xyzverts = None, values = None, cmin = None, cmax = None,
            lim = None, edges = 0) :
   """
   pl3surf (nverts, xyzverts)
         or pl3surf (nverts, xyzverts, values)

     Perform simple 3D rendering of an object created by slice3
     (possibly followed by slice2).  NVERTS and XYZVERTS are polygon
     lists as returned by slice3, so XYZVERTS is sum(NVERTS)-by-3,
     where NVERTS is a list of the number of vertices in each polygon.
     If present, the VALUES should have the same length as NVERTS;
     they are used to color the polygon.  If VALUES is not specified,
     the 3D lighting calculation set up using the light3 function
     will be carried out.  Keywords cmin= and cmax= as for plf, pli,
     or plfp are also accepted.  (If you do not supply VALUES, you
     probably want to use the ambient= keyword to light3 instead of
     cmin= here, but cmax= may still be useful.)
   """

   _draw3 = get_draw3_ ( )
   if type (nverts) == ListType :
      list = nverts
      nverts = list [0]
      xyzverts = array (list [1], copy = 1)
      values = list [2]
      cmin = list [3]
      cmax = list [4]
      edges = list [6]
      ## Scale xyzverts to avoid loss of accuracy
      minx = min (xyzverts [:, 0])
      maxx = max (xyzverts [:, 0])
      miny = min (xyzverts [:, 1])
      maxy = max (xyzverts [:, 1])
      minz = min (xyzverts [:, 2])
      maxz = max (xyzverts [:, 2])
      xyzverts [:, 0] = (xyzverts [:, 0] - minx) / (maxx - minx)
      xyzverts [:, 1] = (xyzverts [:, 1] - miny) / (maxy - miny)
      xyzverts [:, 2] = (xyzverts [:, 2] - minz) / (maxz - minz)
      xyztmp = get3_xy (xyzverts, 1)
      x = xyztmp [:, 0]
      y = xyztmp [:, 1]
      z = xyztmp [:, 2]
      if values == None :
#        xyzverts [:, 0] = x
#        xyzverts [:, 1] = y
#        xyzverts [:, 2] = z
         values = get3_light (xyztmp, nverts)
      [list, vlist] = sort3d (z, nverts)
      nverts = take (nverts, list)
      values = take (values, list)
      x = take (x, vlist)
      y = take (y, vlist)
      _square = get_square_ ( )
      [_xfactor, _yfactor] = get_factors_ ()
      xmax = max (x)
      xmin = min (x)
      ymax = max (y)
      ymin = min (y)
      xdif = xmax - xmin
      ydif = ymax - ymin
      if _xfactor != 1. :
         xmax = xmax + (_xfactor - 1.) * xdif /2.
         xmin = xmin - (_xfactor - 1.) * xdif /2.
      if _yfactor != 1. :
         ymax = ymax + (_yfactor - 1.) * ydif /2.
         ymin = ymin - (_yfactor - 1.) * ydif /2.
      if _square :
         xdif = xmax - xmin
         ydif = ymax - ymin
         if xdif > ydif :
            dif = (xdif - ydif) / 2.
            ymin = ymin - dif
            ymax = ymax + dif
         elif ydif > xdif :
            dif = (ydif - xdif) / 2.
            xmin = xmin - dif
            xmax = xmax + dif

      plfp (values, y, x, nverts, cmin = cmin, cmax = cmax, legend = "",
         edges = edges)
      return [xmin, xmax, ymin, ymax]

   nverts = array (nverts, Int)
   xyzverts = array (xyzverts, Float )

   if shape (xyzverts) [0] != sum (nverts) or sum (less (nverts, 3)) or \
      nverts.typecode () != Int :
      raise _Pl3surfError, "illegal or inconsistent polygon list"
   if values != None and len (values) != len (nverts) :
      raise _Pl3surfError, "illegal or inconsistent polygon color values"

   if values != None :
      values = array (values, Float )

   clear3 ( )
   set3_object ( pl3surf, [nverts, xyzverts, values, cmin, cmax, lim, edges])
   if (_draw3) :
      # Plot the current list if _draw3 has been set.
      call_idler ( )
   if lim :
      tmp = get3_xy (xyzverts, 1)
      return max ( max (abs (tmp [:,0:2])))
   else :
      return None


# ------------------------------------------------------------------------

_Pl3treeError = "Pl3treeError"

def pl3tree (nverts, xyzverts = None, values = None, plane = None,
             cmin = None, cmax = None, split = 1, edges = 0) :

   """
   pl3tree (nverts, xyzverts = None, values = None, plane = None,
      cmin = None, cmax = None)

     Add the polygon list specified by NVERTS (number of vertices in
     each polygon) and XYZVERTS (3-by-sum(NVERTS) vertex coordinates)
     to the currently displayed b-tree.  If VALUES is specified, it
     must have the same dimension as NVERTS, and represents the color
     of each polygon.  If VALUES is not specified, the polygons
     are assumed to form an isosurface which will be shaded by the
     current 3D lighting model; the isosurfaces are at the leaves of
     the b-tree, sliced by all of the planes.  If PLANE is specified,
     the XYZVERTS must all lie in that plane, and that plane becomes
     a new slicing plane in the b-tree.

     Each leaf of the b-tree consists of a set of sliced isosurfaces.
     A node of the b-tree consists of some polygons in one of the
     planes, a b-tree or leaf entirely on one side of that plane, and
     a b-tree or leaf on the other side.  The first plane you add
     becomes the root node, slicing any existing leaf in half.  When
     you add an isosurface, it propagates down the tree, getting
     sliced at each node, until its pieces reach the existing leaves,
     to which they are added.  When you add a plane, it also propagates
     down the tree, getting sliced at each node, until its pieces
     reach the leaves, which it slices, becoming the nodes closest to
     the leaves.

     tree is a 4-element list like this:
      [plane, back_tree, inplane_leaf, front_tree]
      plane= tree [0]  is None if this is just a leaf
                       in which case, only inplane_leaf is not None
      back_tree= tree [1]    is the part behind plane
      inplane_leaf= tree [2] is the part in the plane itself
      front_tree= tree [3]   is the part in front of plane

     This structure is relatively easy to plot, since from any
     viewpoint, a node can always be plotted in the order from one
     side, then the plane, then the other side.

     This routine assumes a "split palette"; the colors for the
     VALUES will be scaled to fit from color 0 to color 99, while
     the colors from the shading calculation will be scaled to fit
     from color 100 to color 199.  (If VALUES is specified as a char
     array, however, it will be used without scaling.)
     You may specifiy a cmin= or cmax= keyword to affect the
     scaling; cmin is ignored if VALUES is not specified (use the
     ambient= keyword from light3 for that case).
   """

#    (ZCM 4/23/97) Add the split keyword. This will determine
#    whether or not to split the palette (half to the isosurfaces
#    for shading and the other half to plane sections for contouring).

#    (ZCM 7/17/97) Add a calculation of the maximum and minimum
#    of everything that is put into the tree. This cures distortion
#    caused by loss of accuracy in orientation calculations.
#    What is now put on the display list is pl3tree and [tree, minmax];
#    both components are passed to _pl3tree to normalize results.

   # avoid overhead of local variables for _pl3tree and _pl3leaf
   # -- I don't know if this is such a big deal
   _draw3 = get_draw3_ ()
   if type (nverts) == ListType :
      _nverts = []
      for i in range (len (nverts)) :
         _nverts.append (nverts [i])
      return _pl3tree (_nverts [0], nverts [1])

   # We need copies of everything, or else arrays get clobbered.
   nverts = array (nverts, Int)
   xyzverts = array (xyzverts, Float )
   if values == "background" :
      values = "bg"
   elif values != None and values != "bg" :
      values = array (values, values.typecode ())
   if plane != None :
      plane = plane.astype (Float)

   if shape (xyzverts) [0] != sum (nverts) or sum (less (nverts, 3)) > 0 or \
      type (nverts [0]) != IntType :
      print "Dim1 of xyzverts ", shape (xyzverts) [0], " sum (nverts) ",\
         sum (nverts), " sum (less (nverts, 3)) ", sum (less (nverts, 3)), \
         " type (nverts [0]) ", `type (nverts [0])`
      raise _Pl3treeError, "illegal or inconsistent polygon list."
   if type (values) == ArrayType and len (values) != len (nverts) and \
      len (values) != sum (nverts) :
      raise _Pl3treeError, "illegal or inconsistent polygon color values"
   if type (values) == ArrayType and len (values) == sum (nverts) :
      # We have vertex-centered values, which for Gist must be
      # averaged over each cell
      list = zeros (sum (nverts), Int)
      array_set (list, cumsum (nverts) [0:-1], ones (len (nverts), Int))
      tpc = values.typecode ()
      values = (histogram (cumsum (list), values) / nverts).astype (tpc)
   if plane != None :
      if (len (shape (plane)) != 1 or shape (plane) [0] != 4) :
         raise _Pl3treeError, "illegal plane format, try plane3 function"

   # Note: a leaf is going to be a list of lists.
   leaf = [ [nverts, xyzverts, values, cmin, cmax, split, edges]]

   ## max and min of current leaf
   minmax = array ( [min (xyzverts [:, 0]), max (xyzverts [:, 0]),
                     min (xyzverts [:, 1]), max (xyzverts [:, 1]),
                     min (xyzverts [:, 2]), max (xyzverts [:, 2])])

   # retrieve current b-tree (if any) from 3D display list
   _draw3_list = get_draw3_list_ ()
   _draw3_n = get_draw3_n_ ()
   try :
      tree = _draw3_list [_draw3_n:]
   except :
      tree = []
   if tree == [] or tree [0] != pl3tree :
      tree = [plane, [], leaf, []]
   else :
      oldminmax = tree [1] [1]
      tree = tree [1] [0]
      ## Find new minmax for whole tree
      minmax = array ( [min (minmax [0], oldminmax [0]),
                        max (minmax [1], oldminmax [1]),
                        min (minmax [2], oldminmax [2]),
                        max (minmax [3], oldminmax [3]),
                        min (minmax [4], oldminmax [4]),
                        max (minmax [5], oldminmax [5])])
      _pl3tree_add (leaf, plane, tree)
      set_multiple_components (1)

   tmp = has_multiple_components ()
   clear3 ()
   set_multiple_components (tmp)
#  plist (tree)
   set3_object (pl3tree, [tree, minmax])
   if (_draw3) :
      ## Plot the current list
      call_idler ( )

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

def split_palette ( * name) :

   """
   split_palette
         or split_palette ("palette_name.gp")
     split the current palette or the specified palette into two
     parts; colors 0 to 99 will be a compressed version of the
     original, while colors 100 to 199 will be a gray scale.
   """

   if len (name) > 0 :
      dum = palette (name [0])
      del dum
   r = zeros (240, 'b')
   g = zeros (240, 'b')
   b = zeros (240, 'b')
   dum = palette (r, g, b, query = 1)
   del dum
   try : # r may be all zeros, in which case the following will fail:
      n = max (nonzero (r)) + 1 # (Figure out how many entries there are)
   except :
      n = 0
   if n < 100 :
      dum = palette ("earth.gp")
      dum = palette (r, g, b, query = 1)
      del dum
      n = max (max (nonzero (r)), max (nonzero (g)),
               max (nonzero (b))) + 1
   newr = zeros (200, 'b')
   newg = zeros (200, 'b')
   newb = zeros (200, 'b')
   newr [0:100] = interp (r [0:n].astype (Float), arange (n, typecode = Float ),
      arange (100, typecode = Float ) * n / 100).astype ('b')
   newg [0:100] = interp (g [0:n].astype (Float), arange (n, typecode = Float ),
      arange (100, typecode = Float ) * n / 100).astype ('b')
   newb [0:100] = interp (b [0:n].astype (Float), arange (n, typecode = Float ),
      arange (100, typecode = Float ) * n / 100).astype ('b')
   newr [100:200] = (arange (100, typecode = Int) * 255 / 99).astype ('b')
   newg [100:200] = (arange (100, typecode = Int) * 255 / 99).astype ('b')
   newb [100:200] = (arange (100, typecode = Int) * 255 / 99).astype ('b')
   palette (newr, newg, newb)

def split_bytscl (x, upper, cmin = None, cmax = None) :

   """
   split_bytscl (x, upper, cmin = None, cmax = None)
     as bytscl function, but scale to the lower half of a split
     palette (0-99, normally the color scale) if the second parameter
     is zero or nil, or the upper half (100-199, normally the gray
     scale) if the second parameter is non-zero.
   """

   x = bytscl (x, cmin = cmin, cmax = cmax, top = 99).astype('b')

   if upper :
      x = x + 100
   return x

def _pl3tree (tree, minmax) :
   #  tree is a 4-element list like this:
   #  [plane, back_tree, inplane_leaf, front_tree]
   #   plane= tree [0]  is None if this is just a leaf
   #                    in which case, only inplane_leaf is not None
   #   back_tree= tree [1]    is the part behind plane
   #   inplane_leaf= tree [2] is the part in the plane itself
   #   front_tree= tree [3]   is the part in front of plane
   if tree == None or tree == [] :
      return None
   if tree [0] == None or tree [0] == [] :
      # only the leaf is non-nil (but not a plane)
      return _pl3leaf ( tree [2], 1, minmax)

   # apply the 3D coordinate transform to two points along the
   # normal of the splitting plane to judge which is in front
   xyz = get3_xy (array ( [ [0., 0., 0.],
                  [tree [0] [0], tree [0] [1], tree [0] [2]]], Float), 1)
   [x, y, z] = [xyz [:, 0], xyz [:, 1], xyz [:, 2]]

   # plot the parts in order toward the current viewpoint
   if z [1] >= z [0] :
      q1 = _pl3tree (tree [1], minmax)
      q2 = _pl3leaf (tree [2], 0, minmax)
      q3 = _pl3tree (tree [3], minmax)
   else :
      q1 = _pl3tree (tree [3], minmax)
      q2 = _pl3leaf (tree [2], 0, minmax)
      q3 = _pl3tree (tree [1], minmax)
   if q1 != None :
      if q2 != None and q3 == None :
         return [min (q2 [0], q1 [0]),
                 max (q2 [1], q1 [1]),
                 min (q2 [2], q1 [2]),
                 max (q2 [3], q1 [3])]
      elif q2 == None and q3 != None :
         return [min (q3 [0], q1 [0]),
                 max (q3 [1], q1 [1]),
                 min (q3 [2], q1 [2]),
                 max (q3 [3], q1 [3])]
      elif q2 != None and q3 != None :
         return [min (q3 [0], q2 [0], q1 [0]),
                 max (q3 [1], q2 [1], q1 [1]),
                 min (q3 [2], q2 [2], q1 [2]),
                 max (q3 [3], q2 [3], q1 [3])]
      else :
         return q1
   elif q2 != None :
      if q3 == None :
         return q2
      else :
         return [min (q2 [0], q3 [0]),
                 max (q2 [1], q3 [1]),
                 min (q2 [2], q3 [2]),
                 max (q2 [3], q3 [3])]
   elif q3 != None :
      return q3
   else :
      return None

## from lp import *

def _pl3leaf (leaf, not_plane, minmax) :
   
   # count number of polys, number of vertices
   _nverts = _xyzverts = 0
   if type (leaf) == ListType and type (leaf [0]) == ListType :
       for i in range (len (leaf)) :
         [_nverts, _xyzverts] = _pl3tree_count ( leaf [i], _nverts, _xyzverts )
   else :
      [_nverts, _xyzverts] = _pl3tree_count ( leaf , _nverts, _xyzverts)

   # accumulate polys and vertices into a single polygon list
   # The type of array required for palettes is "Py_GpColor",
   # which translates to "PyArray_UBYTE", which is selected
   # with a second argument of 'b' to the zeros() function.
## _values = zeros (_nverts, 'b') # See below
   old_nverts = _nverts
   _nverts = zeros (_nverts, Int)
   _x = zeros (_xyzverts, Float )
   _y = zeros (_xyzverts, Float )
   if (not_plane) :
      # Not just straight assignment; make _z a separate copy
      _z = zeros (_xyzverts, Float )
   else :
      _z = None
   _list = 1
   _vlist = 1
   if type (leaf) == ListType and type (leaf [0]) == ListType :
       if leaf [0] [2] != "bg" :
          _values = zeros (old_nverts, 'b')
       else :
          _values = "bg"
       for i in range (len (leaf)) :
         [_list, _vlist, _edges] = _pl3tree_accum ( leaf [i] , not_plane,
            _x, _y, _z, _list, _vlist, _values, _nverts, minmax)
   else :
      if leaf [2] != "bg" :
         _values = zeros (old_nverts, 'b')
      else :
         _values = "bg"
      [_list, _vlist, _edges] = _pl3tree_accum ( leaf , not_plane,
         _x, _y, _z, _list, _vlist, _values, _nverts, minmax)

   # sort the single polygon list
   if not_plane :
      [_list, _vlist] = sort3d (_z, _nverts)
      _nverts = take (_nverts, _list)
      if _values != "bg" :
         _values = take (_values, _list)
      _x = take (_x, _vlist)
      _y = take (_y, _vlist)

   _square = get_square_ ( )
   [_xfactor, _yfactor] = get_factors_ ()
   xmax = max (_x)
   xmin = min (_x)
   ymax = max (_y)
   ymin = min (_y)
   xdif = xmax - xmin
   ydif = ymax - ymin
   if _xfactor != 1. :
      xmax = xmax + (_xfactor - 1.) * xdif /2.
      xmin = xmin - (_xfactor - 1.) * xdif /2.
   if _yfactor != 1. :
      ymax = ymax + (_yfactor - 1.) * ydif /2.
      ymin = ymin - (_yfactor - 1.) * ydif /2.
   if _square :
      xdif = xmax - xmin
      ydif = ymax - ymin
      if xdif > ydif :
         dif = (xdif - ydif) / 2.
         ymin = ymin - dif
         ymax = ymax + dif
      elif ydif > xdif :
         dif = (ydif - xdif) / 2.
         xmin = xmin - dif
         xmax = xmax + dif

   if _values == "bg" :
      _values = None
   plfp (_values, _y, _x, _nverts, legend = "", edges = _edges)
   return [xmin, xmax, ymin, ymax]

def _pl3tree_count (item, _nverts, _xyzverts) :
   return [_nverts + len (item [0]), _xyzverts  + len (ravel (item [1])) / 3]

def _pl3tree_accum (item, not_plane, _x, _y, _z, _list, _vlist, _values,
   _nverts, minmax) :
   # (ZCM 2/24/97) Add _x, _y, _z , _list, _vlist, _nverts as parameters to
   # avoid use of globals. return the new values of _list, _vlist.
   # (ZCM 7/16/97) Return item [6] (whether to show edges)
   # (ZCM 7/17/97) Add parameter minmax to normalize values

   # N. B.:
   # item [0] is nverts
   # item [1] is xyzverts
   # item [2] is values   (if present)
   # item [3] is cmin   (if present)
   # item [4] is cmax   (if present)
   # item [5] is split (1 = split the palette, 0 = do not split)
   # item [6] is edges (1 = show edges, 0 = do not show edges)
   # I have cleaned up what I think is extremely obscure Yorick code
   # apparently designed to avoid some overhead.
   # N. B. avoid splitting the palette if split is 0. (ZCM 4/23/97)

   _xyzverts = array (item [1], copy = 1) # protect copy in tree
   # Normalize copy  (I'm only going to do this if it's an
   # isosurface or is not a plane or has multiple components. 
   # There is a real problem here.
   # You get bad distortion without doing this if one coordinate
   # is many orders of magnitude larger than the others but the
   # others have significant figues. You also get bad distortion
   # by doing this in the case of a single plane section
   # when one coordinate is insignificant with
   # respect to the others and doesn't have significant digits.
   # It is awfully hard to come up with a numerical criterion for this.)
   if item [2] == None or not_plane or has_multiple_components ():
      minx = minmax [0]
      maxx = minmax [1]
      miny = minmax [2]
      maxy = minmax [3]
      minz = minmax [4]
      maxz = minmax [5]
      _xyzverts [:, 0] = (_xyzverts [:, 0] - minx) / (maxx - minx)
      _xyzverts [:, 1] = (_xyzverts [:, 1] - miny) / (maxy - miny)
      _xyzverts [:, 2] = (_xyzverts [:, 2] - minz) / (maxz - minz)
   if  item [2] == None :
      # this is an isosurface to be shaded (no values specified)
      _xyzverts = get3_xy (_xyzverts, 1)
      # accumulate nverts and values
      incr = len (item [0])
      _nverts [ _list - 1: _list - 1 + incr] = item [0]
      if item [5] != 0 :
         _values [ _list - 1: _list - 1 + incr] = split_bytscl (
            get3_light (_xyzverts, item [0]), 1, cmin = 0.0,
            cmax = item [4]).astype ('b')
      else : # no split
         _values [ _list - 1: _list - 1 + incr] = bytscl (
            get3_light (_xyzverts, item [0]), cmin = 0.0,
            cmax = item [4]).astype ('b')
      _list = _list + incr
      # accumulate x, y, and z
      incr = shape (_xyzverts) [0]
      _x [_vlist - 1:_vlist - 1 + incr] = _xyzverts [:, 0]
      _y [_vlist - 1:_vlist - 1 + incr] = _xyzverts [:, 1]
      if not_plane :
         _z [_vlist - 1:_vlist - 1 + incr] = _xyzverts [:, 2]
      _vlist = _vlist + incr
   else :
      # this is to be pseudo-colored since values are given
      if (not_plane) :
         __xyz = get3_xy (_xyzverts, 1)
      else :
         __xyz = get3_xy (_xyzverts, 0)
      # accumulate nverts and values
      incr = len (item [0])
      _nverts [ _list - 1: _list - 1 + incr] = item [0]
      if item [2] != "bg" :
         if (item [2]).typecode () != 'b' :
            if item [5] != 0 :
               _values [ _list - 1: _list - 1 + incr] = split_bytscl (
                  item [2], 0, cmin = item [3], cmax = item [4]).astype ('b')
            else :
               _values [ _list - 1: _list - 1 + incr] = bytscl (
                  item [2], cmin = item [3], cmax = item [4]).astype ('b')
         else :
            _values [ _list - 1: _list - 1 + incr] = item [2]
      _list = _list + incr
      # accumulate x, y, and z
      incr = shape (__xyz) [0]
      _x [_vlist - 1:_vlist - 1 + incr] = __xyz [:, 0]
      _y [_vlist - 1:_vlist - 1 + incr] = __xyz [:, 1]
      if not_plane :
         _z [_vlist - 1:_vlist - 1 + incr] = __xyz [:, 2]
      _vlist = _vlist + incr

   return [_list, _vlist, item [6]]

def _pl3tree_add (leaf, plane, tree) :
   if tree != None and tree != [] and \
      not is_scalar (tree) and tree [0] != None :
      # tree has slicing plane, slice new leaf or plane and descend
      [back, leaf1] = _pl3tree_slice (tree [0], leaf)
      if back :
         if len (tree) >= 2 and tree [1] != None and tree [1] != [] :
            _pl3tree_add (back, plane, tree [1])
         else :
            tree [1] = [None, [], back, []]
      if (leaf1) :
         if len (tree) >= 4 and tree [3] != None and tree [3] != [] :
            _pl3tree_add (leaf1, plane, tree [3])
         else :
            tree [3] = [None, [], leaf1, []]

   elif plane != None :
      # tree is just a leaf, but this leaf has slicing plane
      tree [0] = plane
      tmp = tree [2]
      tree [2] = leaf
      leaf = tmp   # swap new leaf with original leaf
      [back, leaf1] = _pl3tree_slice (plane, leaf)
      if (back) :
         tree [1] = [None, [], back, []]
      if (leaf1) :
         tree [3] = [None, [], leaf1, []]
   else :
      # tree is just a leaf and this leaf has no slicing plane
      tree [2] = leaf + tree [2]
   return

def _pl3tree_slice (plane, leaf) :
   back = frnt = None
   for ll in leaf :
      # each item in the leaf list is itself a list
      nvf = ll [0]
      if nvf != None :
         nvb = array (nvf, copy = 1)
      else :
         nvb = None
      xyzf = ll [1]
      if xyzf != None :
         xyzb = array (xyzf, copy = 1)
      else :
         xyzb = None
      valf = ll [2]
      if valf != None :
         tpc = valf.typecode()
         valb = array (valf, copy = 1)
      else :
         valb = None
      if len (ll) > 4 :
         ll4 = ll [4]
      else :
         ll4 = None
      if len (ll) > 5 :
         ll5 = ll [5]
      else :
         ll5 = 1
      if len (ll) > 6 :
         ll6 = ll [6]
      else :
         ll6 = 0
      [nvf, xyzf, valf, nvb, xyzb, valb] = \
         slice2x (plane, nvf, xyzf, valf)
      if valf != None:
         valf = valf.astype (tpc)
      if valb != None:
         valb = valb.astype (tpc)
      if nvf != None :
         if frnt != None :
            frnt = [ [nvf, xyzf, valf, ll [3], ll4, ll5, ll6]] + frnt
         else :
            frnt = [ [nvf, xyzf, valf, ll [3], ll4, ll5, ll6]]
      if nvb != None :
         if back != None :
            back = [ [nvb, xyzb, valb, ll [3], ll4, ll5, ll6]] + back
         else :
            back = [ [nvb, xyzb, valb, ll [3], ll4, ll5, ll6]]
   return [back, frnt]

_Pl3tree_prtError = "Pl3tree_prtError"

def pl3tree_prt () :
   _draw3_list = get_draw3_list_ ()
   _draw3_n = get_draw3_n_ ()
   if len (_draw3_list) >= _draw3_n :
      tree = _draw3_list [_draw3_n:]
      if tree == None or tree == [] or tree [0] != pl3tree :
         print "<current 3D display not a pl3tree>"
#        raise _Pl3tree_prtError, "<current 3D display not a pl3tree>"
      else :
         tree = tree [1] [0]
         _pl3tree_prt (tree, 0)

def pl3_other_prt(tree = None):
   if tree == None:
      pl3tree_prt ()
   else :
      if tree == None or tree == []:
         print "<current 3D display not a pl3tree>"
      else :
         _pl3tree_prt (tree, 0)

def _pl3tree_prt (tree, depth) :
   if tree == None or tree == [] :
      return
   indent = (" " * (1 + 2 * depth)) [0:-1]
   print indent + "+DEPTH= " + `depth`
   if len (tree) != 4 :
      print indent + "***error - not a tree"
   print indent + "plane= " + `tree [0]`
   back = tree [1]
   list = tree [2]
   frnt = tree [3]
   if back == None or back == [] :
      print indent + "back = []"
   else :
      _pl3tree_prt (back, depth + 1)

   for leaf in list :
      print indent + "leaf length= " + `len (leaf)`
      print indent + "npolys= " + `len (leaf [0])` + \
         ", nverts= " + `sum (leaf [0])` + ", max= " + `max (leaf [0])`
      print indent + "nverts= " + `shape (leaf [1]) [0]` + \
         ", nvals= " + `len (leaf [2])`

   if frnt == None or frnt == [] :
      print  indent + "frnt = []"
   else :
         _pl3tree_prt (frnt, depth + 1)
   print indent + "+DEPTH= " + `depth`

# ------------------------------------------------------------------------

def _isosurface_slicer (m3, chunk, iso_index, _value) :
#  Have to remember here that getv3 can either return an array of 
#  values, or a 2-list consisting of values and the corresponding cell
#  indices, the latter in the case of an irregular grid.
# Note: (ZCM 2/24/97) I have fixed slicers to return the vertex
# information and what used to be the global _xyz3, or None. Hence
# returning the tuple [tmp, None].

   tmp = getv3 (iso_index, m3, chunk)
   if type (tmp) == ListType :
      tmp[0] = tmp [0] - _value
   else :
      tmp = tmp - _value
   return [tmp, None]

def _plane_slicer (m3, chunk, normal, projection) :
   # (ZCM 2/24/97) In all cases, return x as the last element of
   # the tuple. This eliminates the global _xyz3.

   x = xyz3(m3,chunk)
   irregular = type (chunk) == ListType and len (chunk) == 2 \
      or type (chunk) == ArrayType and len (shape (chunk)) == 1 \
      and type (chunk [0]) == IntType
   if irregular :
      # Need to return a list, the first component being the x's,
      # the second being the relative cell list, and the third an offset
      verts = m3 [1] [0]
      cell_offset = 0
      if type (verts) == ListType :
         totals = m3 [1] [3]
         if type (chunk) == ListType :
            fin = chunk [0] [1]
         else :
            fin = chunk [-1]
         for j in range (len (verts)) :
            if fin <= totals [j] :
               break
         if j > 0 :
            cell_offset = totals [j - 1]
      if type (chunk) == ListType :
         clist = arange (0, chunk [0] [1] - chunk [0] [0], typecode = Int)
      else :
         clist = chunk - cell_offset
      # In the irregular case we know x is ncells by 3 by something
      return [ [x [:,0] * normal [0] + x [:,1] * normal [1] + \
         x [:,2] * normal [2] - projection, clist, cell_offset], x]
   elif len (shape (x)) == 5 : # It's ncells X 3 X 2 X 2 X 2
      return [x [:,0] * normal [0] + x [:,1] * normal [1] + \
         x [:,2] * normal [2] - projection, x]
   else :                    # It's 3 X ni X nj X nk
      return [x [0] * normal [0] + x [1] * normal [1] + x [2] * normal [2] -\
         projection, x]

def xyz3 (m3, chunk) :

   """
   xyz3 (m3, chunk)

     return vertex coordinates for CHUNK of 3D mesh M3.  The CHUNK
     may be a list of cell indices, in which case xyz3 returns a
     (dimsof(CHUNK))x3x2x2x2 list of vertex coordinates.  CHUNK may
     also be a mesh-specific data structure used in the slice3
     routine, in which case xyz3 may return a 3x(ni)x(nj)x(nk)
     array of vertex coordinates.  For meshes which are logically
     rectangular or consist of several rectangular patches, this
     is up to 8 times less data, with a concomitant performance
     advantage.  Use xyz3 when writing slicing functions or coloring
     functions for slice3.
   """

   xyz = m3 [0] [0] (m3, chunk)
   return xyz

def xyz3_rect (m3, chunk) :
   m3 = m3 [1]
   if len (shape (chunk)) != 1 :
      c = chunk
      # The difference here is that our arrays are 0-based, while
      # yorick's are 1-based; and the last element in a range is not
      # included in the result array.
      return m3 [1] [:,c [0, 0] - 1:1 + c [1, 0], c [0, 1] - 1:1 + c [1, 1] ,
                     c [0, 2] - 1:1 + c [1, 2]]
   else :
      # Need to create an array of m3 [1] values the same size and shape
      # as what to_corners3 returns.
      # To avoid exceedingly arcane calculations attempting to
      # go backwards to a cell list, this branch returns the list
      # [<function values>, chunk]
      # ???????????? ^^^^^^^^^^^^
      # Then it is trivial for slice3 to find a list of cell
      # numbers in which fi has a negative value.
      dims = m3 [0]
      indices = to_corners3 (chunk, dims [1] + 1, dims [2] + 1)
      no_cells = shape (indices) [0]
      indices = ravel (indices)
      retval = zeros ( (no_cells, 3, 2, 2, 2), Float )
      m30 = ravel (m3 [1] [0, ...])
      retval [:, 0, ...] = reshape (take (m30, indices), (no_cells, 2, 2, 2))
      m31 = ravel (m3 [1] [1, ...])
      retval [:, 1, ...] = reshape (take (m31, indices), (no_cells, 2, 2, 2))
      m32 = ravel (m3 [1] [2, ...])
      retval [:, 2, ...] = reshape (take (m32, indices), (no_cells, 2, 2, 2))
      return retval

_xyz3Error = "xyz3Error"

def xyz3_irreg (m3, chunk) :
   xyz = m3 [1] [1]
   if type (chunk) == ListType and len (chunk) == 2 :
      no_cells = chunk [0] [1] - chunk [0] [0]
      if type (m3 [1] [0]) == ListType :
         totals = m3 [1] [3]
         start = chunk [0] [0]
         fin = chunk [0] [1]
         for i in range (len (totals)) :
            if fin <= totals [i] :
               break
         verts = m3 [1] [0] [i]
         if i > 0 :
            start = start - totals [i - 1]
            fin = fin - totals [i - 1]
         ns = verts [start:fin]
         shp = shape (verts)
      else :
         ns = m3 [1] [0] [chunk [0] [0]:chunk [0] [1]]   # This is a kXnv array
         shp = shape (m3 [1] [0])
   elif type (chunk) == ArrayType and len (shape (chunk)) == 1 and \
      type (chunk [0]) == IntType :
      # chunk is an absolute cell list
      no_cells = len (chunk)
      if type (m3 [1] [0]) == ListType :
         totals = m3 [1] [3]
         fin = max (chunk)
         for i in range (len (totals)) :
            if fin <= totals [i] :
               break
         verts = m3 [1] [0] [i]
         if i > 0 :
            start = totals [i - 1]
         else :
            start = 0
         verts = m3 [1] [0] [i]
         ns = take (verts, chunk - start)
         shp = shape (verts)
      else :
         ns = take (m3 [1] [0], chunk)
         shp = shape (m3 [1] [0])
   else :
      raise _xyz3Error, "chunk parameter has the wrong type."
   if shp [1] == 8 : # hex
      retval = zeros ( (no_cells, 3, 2, 2, 2), Float)
      retval [:, 0] = \
         reshape (take (xyz [0], ravel (ns)), (no_cells, 2, 2, 2))
      retval [:, 1] = \
         reshape (take (xyz [1], ravel (ns)), (no_cells, 2, 2, 2))
      retval [:, 2] = \
         reshape (take (xyz [2], ravel (ns)), (no_cells, 2, 2, 2))
   elif shp [1] == 6 : # prism
      retval = zeros ( (no_cells, 3, 3, 2), Float)
      retval [:, 0] = \
         reshape (take (xyz [0], ravel (ns)), (no_cells, 3, 2))
      retval [:, 1] = \
         reshape (take (xyz [1], ravel (ns)), (no_cells, 3, 2))
      retval [:, 2] = \
         reshape (take (xyz [2], ravel (ns)), (no_cells, 3, 2))
   elif shp [1] == 5 : # pyramid
      retval = zeros ( (no_cells, 3, 5), Float)
      retval [:, 0] = \
         reshape (take (xyz [0], ravel (ns)), (no_cells, 5))
      retval [:, 1] = \
         reshape (take (xyz [1], ravel (ns)), (no_cells, 5))
      retval [:, 2] = \
         reshape (take (xyz [2], ravel (ns)), (no_cells, 5))
   elif shp [1] == 4 : # tet
      retval = zeros ( (no_cells, 3, 4), Float)
      retval [:, 0] = \
         reshape (take (xyz [0], ravel (ns)), (no_cells, 4))
      retval [:, 1] = \
         reshape (take (xyz [1], ravel (ns)), (no_cells, 4))
      retval [:, 2] = \
         reshape (take (xyz [2], ravel (ns)), (no_cells, 4))
   else :
      raise _xyz3Error, "Funny number of cell faces: " + `shp [1]`
   return retval

def xyz3_unif (m3, chunk) :
   m3 = m3 [1]
   n = m3 [1]
   if len (chunk.shape) != 1 :
      c = chunk
      i = c [0] - 1
      dn = c [1] + 1 - i
      xyz = zeros ( (3, dn [0], dn [1], dn [2]), Float)
   else :
      dims = m3 [0]
      nj = dims [1]
      nk = dims [2]
      njnk = nj * nk
      zchunk = chunk % nk
      ychunk = chunk / nk % nj
      xchunk = chunk / njnk
      xyz = zeros ( (len (chunk), 3, 2, 2, 2), Float )
      ijk0 = array ( [zeros ( (2, 2), Int ), ones ( (2, 2), Int )])
      ijk1 = array ( [ [0, 0], [1, 1]], Int )
      ijk1 = array ( [ijk1, ijk1] , Int )
      ijk2 = array ( [ [0, 1], [0, 1]], Int )
      ijk2 = array ( [ijk2, ijk2] , Int )
   if len (n) == 2: # we have dxdydz and x0y0z0
      dxdydz = n [0]
      x0y0z0 = n [1]
      if len (shape (chunk)) != 1:
         # Convert the increment and size into array coordinates
         # -- consecutive values
         xx = arange (dn [0], typecode = Float ) * dxdydz [0] / (dn [0] - 1)
         yy = arange (dn [1], typecode = Float ) * dxdydz [1] / (dn [1] - 1)
         zz = arange (dn [2], typecode = Float ) * dxdydz [2] / (dn [2] - 1)
         xyz [0] = x0y0z0 [0] + i [0] * dxdydz [0] + multiply.outer (
            multiply.outer ( xx, ones (dn [1], Float )),
            ones (dn [2], Float ))
         xyz [1] = x0y0z0 [1] + i [1] * dxdydz [0] + \
            multiply.outer ( ones (dn [0], Float ), \
            multiply.outer ( yy, ones (dn [2], Float )))
         xyz [2] = x0y0z0 [2] + i [2] * dxdydz [0] + \
            multiply.outer ( ones (dn [0], Float ), \
            multiply.outer ( ones (dn [1], Float ), zz))
      else :
         # -- nonconsecutive values
         xyz [:, 0] = add.outer ( xchunk, ijk0) * dxdydz [0] + x0y0z0 [0]
         xyz [:, 1] = add.outer ( ychunk, ijk1) * dxdydz [1] + x0y0z0 [1]
         xyz [:, 2] = add.outer ( zchunk, ijk2) * dxdydz [2] + x0y0z0 [2]
   elif type (n) == ListType and len (n) == 3: # We have three one-dimensional arrays.
      xx = n [0]
      yy = n [1]
      zz = n [2]
      n0 = len (xx)
      n1 = len (yy)
      n2 = len (zz)
      if len (shape (chunk)) != 1:
         # Convert the increment and size into array coordinates
         # -- consecutive values
         xyz [0] = multiply.outer (
            multiply.outer ( xx [i [0]:i [0] + n0],  ones (n1, Float )), \
            ones (n2, Float ))
         xyz [1] =  multiply.outer ( ones (n0, Float ), \
            multiply.outer ( yy [i [1]: i[1] + n1], ones (n2, Float )))
         xyz [2] = multiply.outer ( ones (n0, Float ), \
            multiply.outer ( ones (n1, Float ), zz [i [2]: i[2] + n2]))
      else :
         # -- nonconsecutive values
         xyz [:, 0] = reshape (take (xx, ravel (add.outer (xchunk, ijk0))),
            (len (chunk),  2, 2, 2))
         xyz [:, 1] = reshape (take (yy, ravel (add.outer (ychunk, ijk1))),
            (len (chunk),  2, 2, 2))
         xyz [:, 2] = reshape (take (zz, ravel (add.outer (zchunk, ijk2))),
            (len (chunk),  2, 2, 2))
   return xyz

def to_corners3 (list, nj, nk) :

   """
   to_corners3(list, nj, nk)
     convert an array of cell indices in an (ni-1)-by-(NJ-1)-by-(NK-1)
     logically rectangular grid of cells into the list of
     len(LIST)-by-2-by-2-by-2 cell corner indices in the
     corresponding ni-by-NJ-by-NK array of vertices.
     Note that this computation in Yorick gives an absolute offset
     for each cell quantity in the grid. In Yorick it is legal to
     index a multidimensional array with an absolute offset. In
     Python it is not. However, an array can be flattened if
     necessary.
     Other changes from Yorick were necessitated by row-major
     order and 0-origin indices, and of course the lack of
     Yorick array facilities.
   """

   njnk = nj * nk
   kk = list / (nk - 1)
   list = list + kk + nk * (kk / (nj - 1))
   adder = array ( [ [ [0, 1], [nk, nk + 1]],
                     [ [njnk, njnk + 1], [nk + njnk, nk + njnk + 1]]])
   res = zeros ( (len (list), 2, 2, 2), Int)
   for i in range (len(list)):
      res [i] = adder + list [i]
   return res
