Ticket #10132: trac_10132_final.patch

sage/riemann/parametrized_surface3d.py

@@ -1 +1 @@
 """
 Differential Geometry of Parametrized Surfaces
 
-jvkersch: we need to have some more examples of how to use this class here.
 
 AUTHORS: 
         - Mikhail Malakhaltsev (2010-09-25): initial version
-
+        - Joris Vankerschaver  (2010-10-25): a lot of principal improvements and changes in the code, and in documentation as well
 """    
 #*****************************************************************************
 #       Copyright (C) 2010  Mikhail Malakhaltsev <mikarm@gmail.com>
+#       Copyright (C) 2010  Joris Vankerschaver <joris.vankerschaver@gmail.com>
 #
 #  Distributed under the terms of the GNU General Public License (GPL)
 #                  http://www.gnu.org/licenses/
@@ -18 +18 @@
 # The following is maybe too specific ...
 from vector_functions import *
 
-# mikarm: This is for the moment, at the final version will be replaced by import 
-#attach('/home/prince/MY_SAGE_MODULES/vector_functions.py') 
-
 from sage.structure.sage_object import SageObject
 from sage.calculus.functional import diff
 from sage.functions.other import sqrt
@@ -45 +42 @@
 
 class ParametrizedSurface3D(GenericManifold):
     r"""
-    This is a class for finding basic invariants of surfaces.
+    This class includes methods for calculation the main geometrical objects related 
+    to a parametrized surface such as the first and the second fundamental form, the total (gaussian)
+    and the mean curvature, the geodesic curves, the parallel transport, etc.
 
     The orientation of the surface is determined by the parametrization, that is, the
     natural frame with positive orientation is given by $\partial_1 \vec r$, $\partial_2 \vec r$.         
 
-    jvkersch: This description, as well as the doctests, should be expanded.
 
     INPUT: 
 
@@ -60 +58 @@
 
      - ``name_of_surface`` - string with the name of the surface (optional).
         
-    EXAMPLES::
+    EXAMPLES:
 
-      sage: u, v = var('u,v')
-      sage: paraboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'paraboloid'); paraboloid
-      Parametrized surface ('paraboloid') with equation [u, v, u^2 + v^2]
+      Let us examplify the class by solving several standard tasks of differential geometry of surfaces.
+      First, look how to set a surface. 
+      
+      The elliptic paraboloid:: 
+            
+            sage: u, v = var('u,v')
+            sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
+            sage: eparaboloid
+            Parametrized surface ('elliptic paraboloid') with equation [u, v, u^2 + v^2]
+      
+      An ellipsoid with axes $a$, $b$, $c$:: 
 
-    """
+            sage: a, b,c  = var('a,b,c');
+            sage: u1, u2 = var ('u1,u2');
+            sage: u = [u1,u2]
+            sage: ellipsoid_parametric_equation = vector([a*cos(u1)*cos(u2),b*sin(u1)*cos(u2),c*sin(u2)])
+            sage: ellipsoid = ParametrizedSurface3D(ellipsoid_parametric_equation,u,'ellipsoid')
+            sage: print ellipsoid
+            Parametrized surface ('ellipsoid') with equation (a*cos(u1)*cos(u2), b*sin(u1)*cos(u2), c*sin(u2))
+      
+      One may obtain the latex output of the equation of the surface:: 
 
+            sage: u, v = var('u, v')
+            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
+            sage: print latex(sphere)
+            \left[\cos\left(u\right) \cos\left(v\right), \sin\left(u\right) \cos\left(v\right), \sin\left(v\right)\right]
+            sage: print sphere._latex_()
+            \left[\cos\left(u\right) \cos\left(v\right), \sin\left(u\right) \cos\left(v\right), \sin\left(v\right)\right]
+            sage: print sphere
+            Parametrized surface ('sphere') with equation [cos(u)*cos(v), sin(u)*cos(v), sin(v)]
+      
+      We can plot the surface using the inner method plot:: 
+            
+            sage: ellipsoid_parametric_equation_abc = ellipsoid_parametric_equation.substitute(a=2,b=1.5,c=1) 
+            sage: ellipsoid_abc = ParametrizedSurface3D(ellipsoid_parametric_equation_abc,[u1,u2],'ellipsoid_abc')
+            sage: ellipsoid_abc_plot = ellipsoid_abc.plot((u1,0,2*pi),(u2,-pi/2,pi/2))
+            sage: ellipsoid_abc_plot.show(aspect_ratio=(1,1,1))
+
+      Let us find the natural frame of the parametrization of the ellipsoid. We get the dictionary of vector fields of the natural frame of the ellipsoid:: 
+
+            sage: ellipsoid.natural_frame()
+            {1: (-a*sin(u1)*cos(u2), b*cos(u1)*cos(u2), 0), 2: (-a*sin(u2)*cos(u1), -b*sin(u1)*sin(u2), c*cos(u2))}     
+
+      Now find the normal vector field.   
+      The normal vector field which is the vector product of the vectors of the natural frame is::
+
+            sage: ellipsoid.normal_vector()
+            (b*c*cos(u1)*cos(u2)^2, a*c*sin(u1)*cos(u2)^2, a*b*sin(u2)*cos(u2)) 
+
+      And the unit normal vector field of the elliptic paraboloid is::
+      
+            sage: u, v = var('u,v')
+            sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
+            sage: eparaboloid.normal_vector(1)
+            (-2*u/sqrt(4*u^2 + 4*v^2 + 1), -2*v/sqrt(4*u^2 + 4*v^2 + 1), 1/sqrt(4*u^2 + 4*v^2 + 1))  
+      
+      Now let us find the coefficients of the first  fundamental form of the torus:: 
+
+            u, v = var('u,v') 
+            a, b = var('a,b')
+            sage: torus = ParametrizedSurface3D(((a + b*cos(u))*cos(v),(a + b*cos(u))*sin(v), b*sin(u)),[u,v],'torus')
+            sage: torus.first_fundamental_form_coefficients()
+            {(1, 2): 0, (1, 1): b^2, (2, 1): 0, (2, 2): b^2*cos(u)^2 + 2*a*b*cos(u) + a^2} 
+      
+      We can find the length of a curve on the surface. 
+      For example, let us find the length of the curve $u^1 = t$, $u^2 =t$, $t \in [0,2\pi]$, 
+      on the ellipsoid with axes $a=1$, $b=1.5$ and $c=1$. So we take the curve::
+ 
+            sage: t = var('t')
+            sage: uu1 = t
+            sage: uu2 = t
+
+      Then find the tangent vector::
+
+            sage: duu1 = diff(uu1,t)
+            sage: duu2 = diff(uu2,t)
+            sage: duu = vector([duu1,duu2])
+
+      Then calculate the symbolic expression for the length::
+ 
+            sage: L=sqrt(ellipsoid.first_fundamental_form(duu,duu).substitute(u1=uu1,u2=uu2))
+            sage: integrate(L,t,0,2*pi)
+            integrate(sqrt(2*(a^2 - b^2)*sin(t)^2*cos(t)^2 + c^2*cos(t)^2 + (a^2*cos(t)^2 + b^2*sin(t)^2)*sin(t)^2 + (a^2*sin(t)^2 + b^2*cos(t)^2)*cos(t)^2), t, 0, 2*pi)
+
+      And finally find the numerical value::
+
+            sage: print numerical_integral(L.substitute(a=2,b=1.5,c=1),0,1)[0]
+            2.00127905972
+
+      Find the area of the sphere of radius $R$:: 
+
+            sage: R = var('R');
+            sage: u, v = var('u,v');
+            sage: assume(R>0)
+            sage: assume(cos(v)>0)
+            sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
+            sage: integral(integral(sphere.area_form(),u,0,2*pi),v,-pi/2,pi/2)
+            4*pi*R^2
+
+      We can find an orthonormal frame field $\{e_1, e_2\}$ of a surface and calculate it structure functions $c^k_{ij}$, 
+      where $[e_i,e_j] = \sum_k c^k_{ij} e_k$.
+      Let us first find the orthonormal frame field on the elliptic paraboloid::
+
+            sage: u, v = var('u,v')
+            sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
+            sage: eparaboloid.orthonormal_frame()
+            {1: (1/sqrt(4*u^2 + 1), 0, 2*u/sqrt(4*u^2 + 1)), 2: (-4*u*v/(sqrt(4*u^2 + 1)*sqrt(4*u^2 + 4*v^2 + 1)), sqrt(4*u^2 + 1)/sqrt(4*u^2 + 4*v^2 + 1), 2*v/(sqrt(4*u^2 + 1)*sqrt(4*u^2 + 4*v^2 + 1)))}   
+      
+
+      Now find the inner (local) coordinates of the orthonormal frame field on the elliptic paraboloid::
+            
+            sage: EE = eparaboloid.orthonormal_frame(coordinates='int')
+            sage: E1 = EE[1]; E2 = EE[2]
+            sage: CC = eparaboloid.frame_structure_functions(E1,E2)
+            sage: CC[1,2,1].simplify_full()
+            4*sqrt(4*u^2 + 4*v^2 + 1)*v/(sqrt(4*u^2 + 1)*(16*u^4 + 4*(4*u^2 + 1)*v^2 + 8*u^2 + 1))
+
+
+      We find the Gaussian and Mean Curvatures of the sphere:: 
+
+            sage: u, v = var('u,v')
+            sage: sphere = ParametrizedSurface3D((u,v,sqrt(1-u^2-v^2)),[u,v],'sphere')
+            sage: sphere.gauss_curvature()
+            1
+            sage: sphere.mean_curvature()
+            1
+
+      Another example is to plot gaussian curvature of ellipsoid by color::
+
+            sage: u1, u2 = var('u1,u2');
+            sage: u = [u1,u2]
+            sage: ellipsoid_equation(u1,u2) = [2*cos(u1)*cos(u2),1.5*cos(u1)*sin(u2),sin(u1)]
+            sage: ellipsoid = ParametrizedSurface3D(ellipsoid_equation(u1,u2),u,'ellipsoid')
+            sage: # set intervals for variables and the number of division points
+            sage: interval_1 = (-1.5,1.5)
+            sage: interval_2 = (0,6.28)
+            sage: plot_points = [10,20]
+            sage: # make the arguments array
+            sage: u1_array = [interval_1[0] + cnt*(interval_1[1]-interval_1[0])/plot_points[0] for cnt in range(0,plot_points[0])]
+            sage: u2_array = [interval_2[0] + cnt*(interval_2[1]-interval_2[0])/plot_points[1] for cnt in range(0,plot_points[1])]
+            sage: u_array = [ (uu1,uu2) for uu1 in u1_array for uu2 in u2_array]
+            sage: # Find the gaussian curvature 
+            sage: K(u1,u2) = ellipsoid.gauss_curvature()
+            sage: # Make array of K values
+            sage: K_array = [K(uu[0],uu[1]) for uu in u_array]
+            sage: # Find minimum and max of the gauss curvature
+            sage: K_max = max(K_array)
+            sage: K_min = min(K_array)
+            sage: # Make the array of color coefficients
+            sage: cc_array = [ (ccc - K_min)/(K_max - K_min) for ccc in K_array ]
+            sage: points_array = [ellipsoid_equation(u_array[counter][0],u_array[counter][1]) for counter in range(0,len(u_array)) ]
+            sage: curvature_ellipsoid_plot = sum( point([xx for xx in points_array[counter]],color=hue(cc_array[counter]/2))  for counter in range(0,len(u_array)) )
+            sage: curvature_ellipsoid_plot.show(aspect_ratio=1)
+
+      We can find the principal curvatures and principal directions of the elliptic paraboloid::
+ 
+            sage: u, v = var('u,v')
+            sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
+            sage: pd = eparaboloid.principal_directions()
+
+      We extract the principal curvatures::
+
+            sage: k1 = pd[0][0].simplify_full()
+            sage: k1
+            2*sqrt(4*u^2 + 4*v^2 + 1)/(16*u^4 + 8*u^2 + 16*v^4 + 8*(4*u^2 + 1)*v^2 + 1) 
+            sage: k2 = pd[1][0].simplify_full()
+            sage: k2
+            2/sqrt(4*u^2 + 4*v^2 + 1)
+    
+      and let us check it::
+
+            sage: K = eparaboloid.gauss_curvature().simplify_full()
+            sage: K
+            4/(16*u^4 + 8*u^2 + 16*v^4 + 8*(4*u^2 + 1)*v^2 + 1)  
+            sage: H = eparaboloid.mean_curvature().simplify_full()
+            sage: H
+            2*(2*u^2 + 2*v^2 + 1)/(4*u^2 + 4*v^2 + 1)^(3/2)
+            sage: (K - k1*k2).simplify_full()
+            0
+            sage: (2*H - k1 - k2).simplify_full()
+            0
+
+      And we can find the intrinsic (local coordinates) of the principal directions::
+          
+           sage: pd[0][1] 
+           [(1, v/u)]
+           sage: pd[1][1] 
+           [(1, -u/v)]
+            
+      Also with the ParametrizedSurface3D class on can find the coefficients of the second fundamental form, the shape operator, the rotation on the surface at a given angle, the connection coefficients, also one can calculate numerically the geodesics and the parallel translation along a curve (see the documentation below). 
+          """
+      
     def __init__(self,equation,variables,*name):
         """
         See ``ParametrizedSurface3D`` for full documentation.
@@ -136 +320 @@
                {'designation': name, 'eq': str(self.equation)}
 
 
-    def plot3d(self, urange, vrange):
+    def plot(self, urange, vrange):
         """
         Enable easy plotting directly from the surface class.
 
@@ -150 +334 @@
             sage: u, v = var('u, v')
             sage: eq = [3*u + 3*u*v^2 - u^3, 3*v + 3*u^2*v - v^3, 3*(u^2-v^2)]
             sage: enneper = ParametrizedSurface3D(eq,[u,v],'enneper_surface')
-            sage: enneper.plot3d((u, -5, 5), (v, -5, 5))
+            sage: enneper.plot((u, -5, 5), (v, -5, 5))
         """
         from sage.plot.plot3d.parametric_plot3d import parametric_plot3d 
         P = parametric_plot3d(self.equation, urange, vrange)
@@ -238 +422 @@
         Compute a single component $g_{ij}$ of the first fundamental form.  If the
         parametric representation of the surface is given by the vector function
         $\vec r(u^i)$, where $u^i$, $i = 1, 2$ are curvilinear coordinates, then
-
-         .. math::
-
-           g_{ij} = \frac{\partial \vec r}{\partial u^i} \cdot
-             \frac{\partial \vec r}{\partial u^j}.
+        $g_{ij} = \frac{\partial \vec r}{\partial u^i} \cdot \frac{\partial \vec r}{\partial u^j}$.
 
         INPUT:
 
@@ -343 +523 @@
         Returns the square of the coefficient of the area form on the surface.
         In terms of the coefficients $g_{ij}$ (where $i, j = 1, 2$) of the first fundamental
         form, this invariant is given by 
-
-        .. math::
-
-           A^2 = g_{11}g_{22} - g_{12}^2.
+        $A^2 = g_{11}g_{22} - g_{12}^2$.
 
         See also ``area_form``.
         
@@ -372 +549 @@
         Returns the coefficient of the area form on the surface.  In terms of the
         coefficients $g_{ij}$ (where $i, j = 1, 2$) of the first fundamental form, the
         coefficient of the area form is given by 
-
-        .. math::
-
-           A = \sqrt{g_{11}g_{22} - g_{12}^2}.
+        $A = \sqrt{g_{11}g_{22} - g_{12}^2}$.
 
         See also ``area_form_squared``.
         
@@ -399 +573 @@
         r"""
         Returns the coefficients $g^{ij}$ of the inverse of the fundamental form, 
         as a dictionary.  The inverse coefficients are defined by
-
-         .. math::
-
-               \sum_j g^{ij} g_{jk} = \delta^i_k 
-
+        $\sum_j g^{ij} g_{jk} = \delta^i_k$ 
         with $\delta^i_k$ the Kronecker delta.
 
         OUTPUT:
@@ -650 +820 @@
         Returns the structure functions $c^k_{ij}$ for a frame field $e_1, e_2$, i.e.
         a pair of vector fields on the surface which are linearly independent at
         each point.  The structure functions are defined using the Lie bracket by 
-
-         .. math::
-
-              [e_i,e_j] = c^k_{ij}e_k.
+        $[e_i,e_j] = c^k_{ij}e_k$.
          
         INPUT:
 
@@ -801 +968 @@
         r"""
         Returns the coefficient $h_{ij}$ of the second fundamental form corresponding to the
         index $(i, j)$.  If the equation of the surface is $\vec{r}(u^1, u^2)$, then
-
-         .. math::
-
-            h_{ij} = \vec{r}_{u^i u^j} \cdot \vec{n},
-
+        $h_{ij} = \vec{r}_{u^i u^j} \cdot \vec{n}$,
         where $\vec{n}$ is the unit normal.
 
         INPUT:
@@ -954 +1117 @@
         return (enum/denom).simplify_full()
 
 
-    @cached_method
-    def principal_curvatures(self):
-        """
-        Finds the principal curvatures of the surface 
-
-        OUTPUT:
-        The two principal curvatures, given as a dictionary.
-
-        EXAMPLES::
-                      
-           sage: R = var('R')
-           sage: assume(R>0)
-           sage: u, v = var('u,v')
-           sage: assume(cos(v)>0)
-           sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
-           sage: sphere.principal_curvatures()
-           {1: -1/R, 2: -1/R}
-
-           sage: u, v = var('u,v')
-           sage: R, r = var('R, r')
-           sage: assume(R>r,r>0)
-           sage: torus = ParametrizedSurface3D([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
-           sage: torus.principal_curvatures()
-           {1: -cos(v)/(r*cos(v) + R), 2: -1/r}
-
-        """
-
-        from sage.symbolic.assumptions import assume
-        from sage.symbolic.relation import solve
-        from sage.calculus.var import var
-
-        KK = self.gauss_curvature()
-        HH = self.mean_curvature()
-
-        # jvkersch: when this assumption is uncommented, Sage raises an error stating that the assumption
-        # is redundant... Can we safely omit this, based on some geometric reasoning?
-        
-        # assume(HH**2-KK>=0)
-        # mikarm: This is a problem, I had a lot of trobles here. Of course, this inequality always hold true.
-        # I insert this assumption because sage sometimes, in simplification, uses the complex numbers, though the roots
-        # are, for sure, real.  
-        # This, in turn, causes problems when we substitute coordinates into the expression of principal curvatures.
-        # Unfortunately, I did not manage to tell Sage that they are real (to declare the variables as real). 
-        # Moreover, in a neighborhood of an umbilic point we even cannot "smoothly" order the set of principal curvatures.
-        # So, in general, at present this method is far from the final form.
-          
-
- 
-        x = var('x')
-        sol = solve(x**2 -2*HH*x + KK == 0,x)
-        
-        #k1=var('k1')
-        #k2=var('k2')
-
-        # jvkersch: when I tried to run the previous version of the code, I ran into the problem that if the equation for the principal curvatures had a double root (as in the case of the sphere example in the worksheet), solve returned only one root.  Maybe this is a difference due to having different versions of sage.
-
-        k1 = (x.substitute(sol[0])).simplify_full()
-        if len(sol) == 1:
-            k2 = k1
-        else:
-            k2 = (x.substitute(sol[1])).simplify_full()
-        
-        return {1:k1,2: k2}
-
 
     @cached_method
     def shape_operator_coefficients(self):
@@ -1053 +1152 @@
 
 
 
-    # jvkersch: I've changed the definition of the shape operator to return the matrix
-    # of the shape operator, rather than a matrix times a vector.  This will make it easier
-    # to compute eigenvalues and eigenvectors of the shape operator.
     
     def shape_operator(self):
         r"""
@@ -1087 +1183 @@
             [       -8*u*v/(4*u^2 + 4*v^2 + 1)^(3/2) 2*(4*u^2 + 1)/(4*u^2 + 4*v^2 + 1)^(3/2)]
             sage: S.eigenvalues()
             [2*sqrt(4*u^2 + 4*v^2 + 1)/(16*u^4 + 8*u^2 + 16*v^4 + 8*(4*u^2 + 1)*v^2 + 1), 2/sqrt(4*u^2 + 4*v^2 + 1)]
-            sage: paraboloid.principal_curvatures()
-            {1: 2/(4*u^2 + 4*v^2 + 1)^(3/2), 2: 2/sqrt(4*u^2 + 4*v^2 + 1)}
 
         """
         
-        #vv = vector([xx for xx in v])
         shop = self.shape_operator_coefficients()
         shop_matrix=matrix([[shop[(1,1)],shop[(1,2)]],[shop[(2,1)],shop[(2,2)]]])
         return simplify_matrix(shop_matrix)
 
 
-    @cached_method
+
+
+
     def principal_directions(self):
-        """
-        Finds the principal curvatures and principal directions of the surface 
-
-        OUTPUT:
-
-         - The dictionary of lists [principal curvature, corresponding principal direction]         
-
-        If principal curvatures coincide, gives the warning that the surface is a sphere.
-
-        EXAMPLES::
-
-           sage: u, v = var('u, v')
-           sage: R, r = var('R,r')
-           sage: assume(R>r,r>0)
-           sage: torus = ParametrizedSurface3D([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
-           sage: pdd = torus.principal_directions()
-           sage: pdd[1]
-           [-cos(v)/(r*cos(v) + R), (1, 0)]
-           sage: pdd[2]
-           [-1/r, (0, -(R*r*cos(v) + R^2)/r)]
-
-           sage: R = var('R')
-           sage: assume(R>0)
-           sage: assume(cos(v)>0)
-           sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
-           sage: sphere.principal_directions()
-           'This is a sphere, so any direction is principal'
-
-           sage: catenoid = ParametrizedSurface3D([R*cosh(v)*cos(u), R*cosh(v)*sin(u),v],[u,v],'catenoid')
-           sage: pd = catenoid.principal_directions()
-           sage: pd[1][1]
-           (1, 0)
-           sage: pd[2][1]
-           (0, 1/2*(sqrt(R^2*sinh(v)^2 + 1)*sqrt((4*sinh(v)^2*cosh(v)^2 + 1)*R^4 + 2*(2*cosh(v)^2 - 1)*R^2 + 1)*R*cosh(v) + sqrt(R^2*sinh(v)^2 + 1)*((2*sinh(v)^2*cosh(v) + cosh(v))*R^3 + R*cosh(v)))/(R^4*sinh(v)^4 + 2*R^2*sinh(v)^2 + 1))
-           sage: pd[1][1]*pd[2][1]
-           0
-        """
-        gg = self.first_fundamental_form_coefficients()
-        hh = self.second_fundamental_form_coefficients()
-        kk = self.principal_curvatures()
-        if kk[1]==kk[2]:
-            return "This is a sphere, so any direction is principal"
-        pd1 = simplify_vector([hh[(1,2)]-kk[1]*gg[(1,2)],-hh[(1,1)]+kk[1]*gg[(1,1)]])
-        if pd1==vector([0,0]): 
-           pd1 = vector([1,0]) 
-        pd2 = simplify_vector([hh[(1,2)]-kk[2]*gg[(1,2)],-hh[(1,1)]+kk[2]*gg[(1,1)]])
-        if pd2==vector([0,0]):
-           pd2 = vector([1,0]) 
-        return {1:[kk[1],pd1],2:[kk[2],pd2]}
-
-
-    ### Jvkersch: doctests checked up to this point.
-
-    # jvkersch: alternative definition of principal_directions.  The behavior of this
-    # function is consistent with the eigenvector  function in Sage (since it is just
-    # looking for eigenvectors of the shape operator).  In particular, for a spherical
-    # surface it just returns an arbitrary pair of vectors.
-    def principal_directions_new(self):
         r"""
         Finds the principal curvatures and principal directions of the surface 
 
@@ -1173 +1210 @@
            sage: R, r = var('R,r')
            sage: assume(R>r,r>0)
            sage: torus = ParametrizedSurface3D([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
-           sage: torus.principal_directions_new()
+           sage: torus.principal_directions()
            [(-cos(v)/(r*cos(v) + R), [(1, 0)], 1), (-1/r, [(0, 1)], 1)]
 
         """
@@ -1185 +1222 @@
         r"""
         Computes the connection coefficients or Christoffel symbols $\Gamma^k_{ij}$ of the surface.
         If the coefficients of the first fundamental form are given by $g_{ij}$
-        (where $i, j = 1, 2$), then
-        
-         .. math::
-
-           \Gamma^k_{ij} = \frac{1}{2} g^{kl} \left( \frac{\partial g_{li}}{\partial x^j}
-             - \frac{\partial g_{ij}}{\partial x^l}
-             + \frac{\partial g_{lj}}{\partial x^i} \right).
-
+        (where $i, j = 1, 2$), then 
+        $\Gamma^k_{ij} = \frac{1}{2} g^{kl} \left( \frac{\partial g_{li}}{\partial x^j}
+        - \frac{\partial g_{ij}}{\partial x^l}
+        + \frac{\partial g_{lj}}{\partial x^i} \right)$.
         Here, $(g^{kl})$ is the inverse of the matrix $(g_{ij})$, with $i, j = 1, 2$.
 
         OUTPUT:
@@ -1263 +1296 @@
     def geodesics_numerical(self, p0, v0, tinterval):
         r""" 
         Numerical integration of the geodesic equations.  Explicitly, the
-        geodesic equations are given by
-
-         .. math::
-
-           \frac{d^2 u^i}{dt^2} + \Gamma^i_{jk} \frac{d u^j}{dt} \frac{d u^k}{dt} = 0.
+        geodesic equations are given by 
+        $\frac{d^2 u^i}{dt^2} + \Gamma^i_{jk} \frac{d u^j}{dt} \frac{d u^k}{dt} = 0$.
 
         Solving these equations gives the coordinates $(u^1, u^2)$ on the surface
         of the geodesic.  The coordinates in space can then be found by substituting
@@ -1306 +1336 @@
            sage: assume(cos(q)>0)
            sage: sphere = ParametrizedSurface3D([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
            sage: gg_array = sphere.geodesics_numerical([0.0,0.0],[1.0,1.0],[0,2*pi,5])
-           sage: gg_array[0]
-           [0, [0.000000000000000, 0.000000000000000], [1.00000000000000, 1.00000000000000], [1, 0, 0]]
-           sage: gg_array[1]
-           [1.2566370614359172, [0.76440092815484362, 1.8586224292405702], [-0.28386842687533731, 1.9194187166389509], [-0.204895409519981, 0.692104714174527, 0.692104714458033]]
+           sage: [N(xx[0],digits=4) for xx in gg_array]
+           [0.0000, 1.257, 2.513, 3.770, 5.027, 6.283]  
+           sage: [ [ N(x,digits=4) for x in xx[1] ] for xx in gg_array ]
+           [[0.0000, 0.0000], [0.7644, 1.859], [-0.2876, 3.442], [-0.6137, 5.502], [0.5464, 6.937], [0.3714, 9.025]]   
+           sage: [ [ N(x,digits=4) for x in xx[3] ] for xx in gg_array ]
+           [[1.000, 0.0000, 0.0000], [-0.2049, 0.6921, 0.6921], [-0.9160, -0.2836, -0.2836], [0.5803, -0.5759, -0.5759], [0.6782, 0.5196, 0.5196], [-0.8582, 0.3629, 0.3629]]  
         """ 
 
         u1 = self.variables[1]
@@ -1378 +1410 @@
         r"""
         Numerically solves the equations for parallel translation of a vector along a curve
         on the surface.  Explicitly, the equations for parallel translation are given by
-
-         .. math::
-
-           \frac{d u^i}{dt} + u^j \frac{d c^k}{dt} \Gamma^i_{jk} = 0,
-
+        $\frac{d u^i}{dt} + u^j \frac{d c^k}{dt} \Gamma^i_{jk} = 0$,
         where $\Gamma^i_{jk}$ are the connection coefficients of the surface, the vector to
         be transported has components $u^j$ and the curve along which to transport has
         components $c^k$.