PSphDistance - Distance on projective sphere.
   
   USAGE:

   dist = PSphDistance(pt, ptlist)

   INPUT:

   pt     is 3 x 1, 
          a point on the unit sphere (S^2)
   ptlist is 3 x n, 
          a list of points on the unit sphere

   OUTPUT:

   dist is 1 x n, 
        the distance from `pt' to each point in the list

   NOTES:

   *  The distance between two points on the sphere is the angle
      in radians between the two vectors.  On the projective
      sphere, antipodal points are considered equal, so the 
      distance is the minimum of the distances obtained by
      including the negative of pt as well.

 PSphGaussian - Gaussian distribution for smoothing on projective sphere.
   
   USAGE:

   fsm = PSphGaussian(center, pts, stdev)

   INPUT:

   center is 3 x 1, 
          the center of the distribution
   pts    is 3 x n, 
          a list of points on the sphere; antipodal
          points are considered equal
   stdev  is 1 x 1, 
          the (1D) standard deviation

   OUTPUT:

   fsm is 1 x n, 
       the list of values at each point of pts

   Notes:  

   *  The result is not normalized, so this may have to be
      done after the fact.
   *  The distribution is a 1-D normal distribution applied
      to the distance function on the projective sphere.
   *  The actual scaling factor to give unit integral over 
      the projective sphere is not computed; the result
      is not normalized.

 SphBaseMesh - Generate base mesh for spheres.

   USAGE:

   mesh = SphBaseMesh(dim)
   mesh = SphBaseMesh(dim, 'param', 'value')

   INPUT:

   dim is a positive integer,
       the dimension of the sphere (2 for the usual sphere S^2)

   These arguments can be followed by a list of
   parameter/value pairs which specify keyword options.
   Available options include:

   'Hemisphere'  'on'|'off'
                 to mesh only the upper hemisphere
   

   OUTPUT:

   mesh is a MeshStructure,
        on the sphere of the specified dimension

   NOTES:

 SphCrdMesh - Generate a hemisphere mesh based on spherical coordinates.
   
   USAGE:

   smesh = SphCrdMesh(ntheta, nphi)

   INPUT:

   ntheta is a positive integer,
          the number of subdivisions in theta
   nphi   is a positive integer,
          the number of subdivisions in phi

   OUTPUT:

   smesh is a MeshStructure,
         on the hemisphere (H^2)

 SphDifferential - Compute differential of mapping to sphere.
   
   USAGE:

   diff = SphDifferential(mesh, refpts)
   
   INPUT:

   mesh   is a mesh,
          on a sphere of any dimension
   refpts is d x n, 
          a list of points in the reference element,
          usually the quadrature points, given in barycentric
          coordinates

   OUTPUT:

   diff is d x (d-1) x nq, 
        a list of tangent vectors at each reference
        point for each element; nq is the number of global 
        quadrature points, that is n x ne, where ne is the
        number of elements
        

 SphDistance - Find distance on sphere between a fixed point and others.
   
   USAGE:

   dist = SphDistance(pt, ptlist)

   INPUT:

   pt     is 3 x 1, 
          a point on the unit sphere (S^2)
   ptlist is 3 x n, 
          a list of points on the unit sphere (S^2)

   OUTPUT:

   dist is 1 x n, 
        the distance from `pt' to each point on the list

   NOTES:

   *  The distance is just the angle between the two vectors.

 SphDistanceFunc - Return half sum of squared distances on sphere.
   
   USAGE:

   f           = SphDistanceFunc(x, pts, @Sofx)
   [f, gf]     = SphDistanceFunc(x, pts, @Sofx)
   [f, gf, Hf] = SphDistanceFunc(x, pts, @Sofx)

   INPUT:

   x    is d x 1, 
        a point in parameter space 
   pts  is (d+1) x n, 
        a list of n points on the sphere
   Sofx is a function handle, 
        returning parameterization component quantities 
        (function, gradient and Hessian)

   OUTPUT:

   f  is a scalar, 
      the objective function at x
   gf is a vector, 
      the gradient of f at x
   Hf is a matrix, 
      the Hessian of f at x

   NOTES:

   *  See MisorientationStats

 SphGQRule - Global quadrature rule for sphere.
   
   USAGE:

   [gqrule, l2ip] = SphGQRule(mesh, qrule)

   INPUT:

   mesh  is a MeshStructure,
         on a sphere of any dimension
   qrule is a QRuleStructure,
         on the reference element of the mesh

   OUTPUT:

   gqrule is a QRuleStructure,
          for the entire mesh
   l2ip   is n x n, (sparse)
          it gives the l2 inner product in terms of function
          values at the nodal points

 SphGaussian - Gaussian distribution on angular distance.
   
   USAGE:

   fsm = SphGaussian(cen, pts, stdev)

   INPUT:

   cen   is 3 x 1, 
         the center of the distribution
   pts   is 3 x n, 
         a list of points on the sphere at which to
         evaluate the result
   stdev is a scalar, 
         the standard deviation

   OUTPUT:

   fsm is 1 x n, 
       the list of values at each point of `pts'

   NOTES:

   *  This returns the values of a 1D Gaussian applied
      to spherical distance from the center point
   *  The result is not normalized to have unit integral
      over the sphere.

 SphH1SIP -- H^1 semi-inner product on sphere

   USAGE:

   h1sip = SphH1SIP(mesh, qrule)

   INPUT:

   mesh  is a MeshStructure,
         a mesh on a sphere of any dimension
   qrule is a QRuleStructure,
         a quadrature rule for the reference simplex

   OUTPUT:

   h1sip is n x n, (sparse)
         it is the matrix of the H^1 semi-inner product,
         where n is the number of independent degrees of
         freedom associated with the mesh
        

 SphereAverage - Find `average' of list of points on sphere.
   
   USAGE:

   avg = SphereAverage(pts)
   [avg, optdat] = SphereAverage(pts, Pzation, nlopts)

   INPUT:

   pts     is m x n, 
           a list of n points in R^m of unit length
   Pzation is a function handle,
           (see note below)
   x0      is the initial guess,
           in given parameterization
   nlopts  are options to be passed to the nonlinear minimizer.

   OUTPUT:

   avg is m x 1, 
          is a unit vector representing the "average" of `pts'
   optdat is a cell array, 
          with three members,  {fval, exitflag, output}
          (see documentation for `fminunc')

   NOTES:

   *  If only one argument is given, the average returned is 
      the arithmetic average of the points.  If all three arguments
      are given, then the average is computed using unconstrained
      minimization of the sum of squared angles from the data points,
      using the parameterization specified and the options given.
      
   *  See the matlab builtin `fminunc' for details.

   *  This routine needs to be fixed.  Currently it uses the
      parameterization given by `SpherePZ' instead of the 
      function handle `PZation'.

 SpherePZ - Generate point, gradients and Hessians of map to sphere.
   
   USAGE:

   sk           = SpherePZ(x)
   [sk, gk]     = SpherePZ(x)
   [sk, gk, Hk] = SpherePZ(x)

   INPUT:

   x is d x 1, 
     a vector with norm <= 1

   OUTPUT:

   sk is e x 1, 
      a point on the sphere (sqrt(1-x^2), x)
   gk is d x e, 
      the gradients of each component of sk
   Hk is d x d x e, 
      the Hessians of each component of sk

 XYZOfThetaPhi - Map spherical coordinates to sphere.
   
   USAGE:

   xyz = XYZOfThetaPhi(thetaphi)

   INPUT:

   thetaphi is 2 x n, 
            the spherical coordinates for a list of n points; 
            theta is the angle that the projection onto x-y plane 
            makes with the x-axis, and phi is the angle with z-axis

   OUTPUT:

   xyz is 3 x n, 
       the Cartesian coordinates of the points described by (theta, phi)

   NOTES: 

   *  The matlab builtin `sph2cart' could also be used, but 
      the convention for phi is different (90 degrees minus thisone).