BuildOdfPfMatrix - Build ODF/PF matrix in pieces.
   
   USAGE:

   opm = BuildOdfPfMatrix(hkl, mesh, sym, pts, div, invpf, block)

   INPUT:

   hkl,  mesh, sym, pts, div, invpf 
         are the same as in `OdfPfMatrix', except that `invpf' is 
         no longer optional 
   block is a postive integer, 
         the number of pole figure points to be processed at a
         time;  the pole figure points can be processed 
         independently to control memory allocation

   OUTPUT:

   See `OdfPfMatrix'

   NOTES:

   *  See OdfPfMatrix for further documentation.

   *  This is preferrd to OdfPfMatrix for problems with
      many divisions per fiber and many pole figure points.

 FiberCoordinates - Find mesh coordinates for a given fiber.
   
   USAGE:

   [fele, fcrd] = FiberCoordinates(fib, mesh)

   INPUT:

   fib is a 3 x ndiv x npole,
       it is an array with each page giving the fiber over 
       the n'th pole point; usually the result of FiberOfPoint
   mesh is a MeshStructure 
       on the fundamental region

   OUTPUT:

   fele is ndiv x npole, 
        the list of elements containing the fiber points.
   fcrd is 4 x ndiv x npole, 
        the barycentric coordinates of the fiber points

  Notes:

  *  The call to the matlab builtin `tsearchn' can fail, possibly
     due to the fact that the meshes are not necessarily Delaunay
     tesselations.

 FiberOfPoint - Find fiber above point for specified pole figure.
   
   USAGE:

   fib = FiberOfPoint(p, h, ndiv, qsym)
   fib = FiberOfPoint(p, h, ndiv, qsym, invfib)
   
   INPUT:

   p    is 3 x n,  
        an array of points on the sphere 
   h    is 3 x 1,  
        a specified pole direction; if `invfib' is nonzero,
        this is interpreted as a crystal direction; otherwise
        it is interpreted as a sample direction
   ndiv is a positive integer, 
        the number of equally spaced divisions along the fiber
   qsym is 4 x m, 
        the quaternion representation of the symmetry group
   invfib is a scalar, (optional) 
        a flag which, if nonzero, produces the inverse fiber 
   
   OUTPUT:

   fib is 3 x ndiv x n, 
       the column represents the Rodrigues
       vector of a point on the fiber, the row spans the points
       on the fiber, and the page spans the points on the sphere
   
   NOTES:

   * `h' need not be a unit vector; it is normalized here

 FindFiber - Find fiber for given crystal and sample directions.
  
   USAGE:

   fibs = FindFiber(c, s, ndiv)

   INPUT:

   c is 3 x m,
        an array of crystal directions
   s is 3 x n,
        an array of sample directions;
        if n > 1 then m must be 1, and vice versa
   ndiv is a postive integer,
        the number of points to return per fiber

   OUTPUT:

   fibs is 4 x ndiv x nfibs,
        an array of equally spaced points (quaternions) along each
        specified fiber; `nfibs' is either `m' or `n' above

   NOTES:

   * The (c, s) fiber is the collection of rotations R
     such that Rc = s

   * This routine can be called with many crystal directions
     and a single sample (as for inverse pole figures) or
     with many sample directions and a single crystal direction
     (as for pole figures)

 MeshCoordinates - find elements and barycentric coordinates of points
   
   USAGE:

   [els, crds] = MeshCoordinates(mesh, pts)

   INPUT:

   mesh is a MeshStructure,
        it should be Delaunay, but it may be okay anyway if it is
        not too irregular
   pts  is m x n,
        a list of n m-dimensional points

   OUTPUT:

   els is an n-vector, (integers)
       the list of element numbers containing each point
   crds is (m+1) x n,  
       the barycentric coordinates of each point

   NOTES:

   *  The call to the matlab builtin `tsearchn' can fail if 
      the mesh is not Delaunay.

 OdfPfMatrix - Find matrix relating ODF to pole figure.
   
   USAGE:

   opm = OdfPfMatrix(hkl, mesh, sym, pts, div)
   opm = OdfPfMatrix(hkl, mesh, sym, pts, div, invpf)

   INPUT:

   hkl   is a 3-vector, 
         the crystal direction specifying the pole figure
   mesh  is a MeshStructure,
         on orientation space
   sym   is 4 x s, 
         the symmetry group in quaternions
   pts   is 3 x p, 
         a list of p points on the sphere (S^2)
   div   is a positive integer, 
         the number of divisions per fiber to use in
         calculating the fiber integral
   invpf is a scalar, (optional, default = 0) 
         a nonzero value flags causes computation of inverse pole 
         figure matrix in which the `hkl' is interpreted
         as a fixed sample direction and `pts' is an array
         of crystal directions

   OUTPUT:

   opm is p x n, (sparse) 
       the matrix which takes nodal point values on the
       fundamental region to pole figure or inverse
       values at the specified points

   NOTES:

   *  opm acts on the "reduced" set of nodes, the set of nodes
      in which equivalent nodes are combined into a single
      degree of freedom
   *  this routine can be very memory intensive; you may need to
      build the matrix in pieces; if you have a lot of points
      and a large value of `div', then you should break the points
      into smaller groups 
   *  the ODF-PF matrix preserves mean value, so that 
      an ODF in MUD (multiples of uniform) maps to a pole
      figure also in MUD.

 QFiberOfPoint - Find fiber above point in quaternions.
   
   USAGE:

   qfib = QFiberOfPoint(p, h, ndiv, qsym)
   qfib = QFiberOfPoint(p, h, ndiv, qsym, invfib)
   
   INPUT:

   p    is 3 x n,  
        an array of points on the sphere 
   h    is 3 x 1,  
        a specified pole direction; if `invfib' is nonzero,
        this is interpreted as a crystal direction; otherwise
        it is interpreted as a sample direction
   ndiv is a positive integer, 
        the number of equally spaced divisions along the fiber
   qsym is 4 x m, 
        the quaternion representation of the symmetry group
   invfib is a scalar, (optional, default: 0) 
        a flag which, if nonzero, produces the inverse fiber 
   
   OUTPUT:

   qfib is 4 x ndiv x n, 
        each column gives the quaternion representation of a point 
        on the fiber; the row spans the points on the fiber;
        and the page (third index) spans the points on the sphere
   
   NOTES:
      
   * h need not be a unit vector; it is normalized here