## Automatically adapted for scipy Oct 31, 2005 by

# $Id: slice3.py 2183 2006-08-29 10:30:44Z oliphant $
# 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 2183 2006-08-29 10:30:44Z oliphant $
#

# 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 scipy import *
from shapetest import *
from types import *
from pl3d import *
from numpy import *
from gistC import *
from gistfuncs import *
from numpy.oldnumeric import nonzero

 #
 # 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,axis=0))
    newnorm [3] = sum (multiply (normal, point),axis=0)
    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.typed  == y.typed == \
       z.typed == 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,axis=0) 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,axis=0) != 0) :
                    clist = take (fs [1], nonzero (critical_cells),axis=0)
                    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,axis=0) != 0) :
                    clist = take (fs [1], nonzero (critical_cells),axis=0)
                    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,axis=0) != 0) :
                    clist = take (fs [1], nonzero (critical_cells),axis=0)
                    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,axis=0) != 0) :
                    clist = take (fs [1], nonzero (critical_cells),axis=0)
                    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,axis=0) > 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,axis=0)
                if got_xyz :
                    _xyz3 = take (_xyz3, clist,axis=0)
                if col :
                    col = take (col, clist,axis=0)
            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,axis=0),\
                   (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,axis=0), (no_cells, 2, 2, 2))
                    new_xyz3 [:, 1, ...] = reshape (take (ravel (_xyz3 [1, ...]),\
                       indices,axis=0), (no_cells, 2, 2, 2))
                    new_xyz3 [:, 2, ...] = reshape (take (ravel (_xyz3 [2, ...]),\
                       indices,axis=0), (no_cells, 2, 2, 2))
                    _xyz3 = new_xyz3
                    del new_xyz3
                if col != None :
                    col = reshape (take (ravel (col), indices,axis=0), (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,axis=0) + _no_verts [i] * cells
        upper = take (_upper_vert [i], edges,axis=0) + _no_verts [i] * cells
        fsl = take (ravel (fs), lower,axis=0)
        fsu = take (ravel (fs), upper,axis=0)
        # 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,axis=0) * fsu - \
           take (ravel (xyz [:, 0]), upper,axis=0) * fsl), (len (lower),))
        new_xyz [:, 1] = reshape ( (take (ravel (xyz [:, 1]), lower,axis=0) * fsu - \
           take (ravel (xyz [:, 1]), upper,axis=0) * fsl), (len (lower),))
        new_xyz [:, 2] = reshape ( (take (ravel (xyz [:, 2]), lower,axis=0) * fsu - \
           take (ravel (xyz [:, 2]), upper,axis=0) * 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,axis=0) * fsu - \
               take (ravel (col), upper,axis=0) * 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)),axis=0))

        # broadcast the cell's pattern onto each of its sliced edges
        pattern = take (pattern, list / _no_edges [i],axis=0)
        # 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,axis=0) + 4 * _no_edges [i] * cells
        order = argsort (pattern)
        xyz1 = zeros ( (len (order), 3), Float )
        xyz1 [:,0] = take (ravel (xyz [:,0]), order,axis=0)
        xyz1 [:,1] = take (ravel (xyz [:,1]), order,axis=0)
        xyz1 [:,2] = take (ravel (xyz [:,2]), order,axis=0)
        xyz = xyz1
        if col != None :
            col = take (col, order,axis=0)
        edges = take (edges, order,axis=0)
        pattern = take (pattern, order,axis=0)
        # 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),axis=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, axis=0),axis=0)
        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,axis=0) 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, dtype = Float )
            y = arange (ncy, dtype = Float )
    elif len (args) == 3 :
        # must be the (nxny, dxdy, x0y0, z...) form
        ncx = xyz [0] + 1
        ncy = xyz [1] + 1
        x = arange (ncx, dtype = Float ) * args [0] [0] + args [1] [0]
        y = arange (ncy, dtype = 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, dtype = Int) * (ncy)
    ncyy = arange (ncy - 1, dtype = 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)),axis=0) +
               take (col, ravel (add.outer ( ncxx, ncyy + 1)),axis=0) +
               take (col, ravel (add.outer ( ncxx + ncy, ncyy + 1)),axis=0) +
               take (col, ravel (add.outer ( ncxx + ncy, ncyy)),axis=0))
        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)),axis=0),
               take (col, ravel (add.outer ( ncxx, ncyy + 1)),axis=0),
               take (col, ravel (add.outer ( ncxx + ncy, ncyy + 1)),axis=0),
               take (col, ravel (add.outer ( ncxx + ncy, ncyy)),axis=0)])))
        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)),axis=0),
           take (newx, ravel (add.outer ( ncxx, ncyy + 1)),axis=0),
           take (newx, ravel (add.outer ( ncxx + ncy, ncyy + 1)),axis=0),
           take (newx, ravel (add.outer ( ncxx + ncy, ncyy)),axis=0)])))
        newy = ravel (y)
        xyzverts [:, 1] = ravel (transpose (array ( [
           take (newy, ravel (add.outer ( ncxx, ncyy)),axis=0),
           take (newy, ravel (add.outer ( ncxx, ncyy + 1)),axis=0),
           take (newy, ravel (add.outer ( ncxx + ncy, ncyy + 1)),axis=0),
           take (newy, ravel (add.outer ( ncxx + ncy, ncyy)),axis=0)])))
    newz = ravel (z)
    xyzverts [:, 2] = ravel (transpose (array ( [
       take (newz, ravel (add.outer ( ncxx, ncyy)),axis=0),
       take (newz, ravel (add.outer ( ncxx, ncyy + 1)),axis=0),
       take (newz, ravel (add.outer ( ncxx + ncy, ncyy + 1)),axis=0),
       take (newz, ravel (add.outer ( ncxx + ncy, ncyy)),axis=0)])))

    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),axis=0) + 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),
                     dtype = 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),
                         dtype = 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),
                     dtype = 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),
                            dtype = 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,axis=0), (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,axis=0), (no_cells, 2, 2, 2)),
                arange (0, no_cells, dtype = Int), oldstart]
    elif tc == 6 : # pyramids
        return [ reshape (take (fi [i], indices,axis=0), (no_cells, 3, 2)),
                arange (0, no_cells, dtype = Int), oldstart]
    else : # tetrahedron or pyramid
        return [ reshape (take (fi [i], indices,axis=0), (no_cells, tc)),
                arange (0, no_cells, dtype = 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,axis=0),
               (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,axis=0)
            if l == None :
                return 0.125 * sum (transpose (reshape (corners, (no_cells, 8))),axis=0)
            else :
                # interpolate corner values to get edge values
                corners = (take (corners, l,axis=0) * fsu -
                   take (corners, u,axis=0) * 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,axis=0)
        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,axis=0)
        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,axis=0)
        else :
            raise _Getc3Error, "chunk argument is incomprehensible."

    shp = shape (indices)
    no_cells = shp [0]
    indices = ravel (indices)
    corners = take (fi [i - 1], indices,axis=0)
    if l == None :
        return (1. / shp [1]) * transpose ((sum (transpose (reshape (corners,
           (no_cells, shp [1]))) [0:shp [1]],axis=0)))
    else :
        # interpolate corner values to get edge values
        corners = (take (corners, l,axis=0) * fsu -
           take (corners, u,axis=0) * 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, dtype = 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, dtype = Float) * \
               (vcmax - vcmin) / (n + 1)
        elif scale == "log" :
            vc = vcmin + exp (arange (1, n + 1, dtype = 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.dtype == Float :
        n = len (contours)
        vc = sort (contours)
    else :
        raise _ContourError, "Incorrect contour specification."
    if split == 0 :
        colors = (arange (n + 1, dtype = Float) * (199. / n)).astype ('B')
    else :
        colors = (arange (n + 1, dtype = 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, \
               dtype = Float) * \
               (vcmax - vcmin) / (n + 1)
        elif scale == "log" :
            vc = vcmin + exp (arange (1, n + 1, \
               dtype = 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.dtype == Float :
        n = len (contours)
        vc = sort (contours)
    else :
        raise _ContourError, "Incorrect contour specification."
    if split == 0 :
        colors = (arange (n + 1, dtype = Float) * (199. / n)).astype ('B')
    else :
        colors = (arange (n + 1, dtype = 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,axis=0)-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,axis=0)
        values = take (values, list,axis=0)
        x = take (x, vlist,axis=0)
        y = take (y, vlist,axis=0)
        _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,axis=0) or sum (less (nverts, 3),axis=0) or \
       nverts.dtype != 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,axis=0) 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.dtype)
    if plane != None :
        plane = plane.astype (Float)

    if shape (xyzverts) [0] != sum (nverts,axis=0) or sum (less (nverts, 3),axis=0) > 0 or \
       type (nverts [0]) != IntType :
        print "Dim1 of xyzverts ", shape (xyzverts) [0], " sum (nverts,axis=0) ",\
           sum (nverts,axis=0), " sum (less (nverts, 3),axis=0) ", sum (less (nverts, 3),axis=0), \
           " 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,axis=0) :
        raise _Pl3treeError, "illegal or inconsistent polygon color values"
    if type (values) == ArrayType and len (values) == sum (nverts,axis=0) :
        # We have vertex-centered values, which for Gist must be
        # averaged over each cell
        list = zeros (sum (nverts,axis=0), Int)
        array_set (list, cumsum (nverts,axis=0) [0:-1], ones (len (nverts), Int))
        tpc = values.dtype
        values = (histogram (cumsum (list,axis=0), 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, dtype = Float ),
       arange (100, dtype = Float ) * n / 100).astype ('B')
    newg [0:100] = interp (g [0:n].astype (Float), arange (n, dtype = Float ),
       arange (100, dtype = Float ) * n / 100).astype ('B')
    newb [0:100] = interp (b [0:n].astype (Float), arange (n, dtype = Float ),
       arange (100, dtype = Float ) * n / 100).astype ('B')
    newr [100:200] = (arange (100, dtype = Int) * 255 / 99).astype ('B')
    newg [100:200] = (arange (100, dtype = Int) * 255 / 99).astype ('B')
    newb [100:200] = (arange (100, dtype = 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,axis=0)
        if _values != "bg" :
            _values = take (_values, _list,axis=0)
        _x = take (_x, _vlist,axis=0)
        _y = take (_y, _vlist,axis=0)

    _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]).dtype.char != '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.dtype.char
            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],axis=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], dtype = 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,axis=0), (no_cells, 2, 2, 2))
        m31 = ravel (m3 [1] [1, ...])
        retval [:, 1, ...] = reshape (take (m31, indices,axis=0), (no_cells, 2, 2, 2))
        m32 = ravel (m3 [1] [2, ...])
        retval [:, 2, ...] = reshape (take (m32, indices,axis=0), (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,axis=0)
            shp = shape (verts)
        else :
            ns = take (m3 [1] [0], chunk,axis=0)
            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),axis=0), (no_cells, 2, 2, 2))
        retval [:, 1] = \
           reshape (take (xyz [1], ravel (ns),axis=0), (no_cells, 2, 2, 2))
        retval [:, 2] = \
           reshape (take (xyz [2], ravel (ns),axis=0), (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),axis=0), (no_cells, 3, 2))
        retval [:, 1] = \
           reshape (take (xyz [1], ravel (ns),axis=0), (no_cells, 3, 2))
        retval [:, 2] = \
           reshape (take (xyz [2], ravel (ns),axis=0), (no_cells, 3, 2))
    elif shp [1] == 5 : # pyramid
        retval = zeros ( (no_cells, 3, 5), Float)
        retval [:, 0] = \
           reshape (take (xyz [0], ravel (ns),axis=0), (no_cells, 5))
        retval [:, 1] = \
           reshape (take (xyz [1], ravel (ns),axis=0), (no_cells, 5))
        retval [:, 2] = \
           reshape (take (xyz [2], ravel (ns),axis=0), (no_cells, 5))
    elif shp [1] == 4 : # tet
        retval = zeros ( (no_cells, 3, 4), Float)
        retval [:, 0] = \
           reshape (take (xyz [0], ravel (ns),axis=0), (no_cells, 4))
        retval [:, 1] = \
           reshape (take (xyz [1], ravel (ns),axis=0), (no_cells, 4))
        retval [:, 2] = \
           reshape (take (xyz [2], ravel (ns),axis=0), (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], dtype = Float ) * dxdydz [0] / (dn [0] - 1)
            yy = arange (dn [1], dtype = Float ) * dxdydz [1] / (dn [1] - 1)
            zz = arange (dn [2], dtype = 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)),axis=0),
               (len (chunk),  2, 2, 2))
            xyz [:, 1] = reshape (take (yy, ravel (add.outer (ychunk, ijk1)),axis=0),
               (len (chunk),  2, 2, 2))
            xyz [:, 2] = reshape (take (zz, ravel (add.outer (zchunk, ijk2)),axis=0),
               (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