Ticket #10132: trac_10132_final.patch

File trac_10132_final.patch, 25.5 KB (added by mikarm, 10 years ago)
  • sage/riemann/parametrized_surface3d.py

    exporting patch:
    # HG changeset patch
    # User Mikhail Malakhaltsev
    # Date 1299032807 18000
    # Node ID d1259922e233eae5a77e57e4b808734c3ce9d30a
    # Parent  8d868706a5852b3b10c536f21d9e0ea577bcbd65
    Trac 10132: final version
    
    diff -r 8d868706a585 -r d1259922e233 sage/riemann/parametrized_surface3d.py
    a b  
    11"""
    22Differential Geometry of Parametrized Surfaces
    33
    4 jvkersch: we need to have some more examples of how to use this class here.
    54
    65AUTHORS:
    76        - Mikhail Malakhaltsev (2010-09-25): initial version
    8 
     7        - Joris Vankerschaver  (2010-10-25): a lot of principal improvements and changes in the code, and in documentation as well
    98"""   
    109#*****************************************************************************
    1110#       Copyright (C) 2010  Mikhail Malakhaltsev <mikarm@gmail.com>
     11#       Copyright (C) 2010  Joris Vankerschaver <joris.vankerschaver@gmail.com>
    1212#
    1313#  Distributed under the terms of the GNU General Public License (GPL)
    1414#                  http://www.gnu.org/licenses/
     
    1818# The following is maybe too specific ...
    1919from vector_functions import *
    2020
    21 # mikarm: This is for the moment, at the final version will be replaced by import
    22 #attach('/home/prince/MY_SAGE_MODULES/vector_functions.py')
    23 
    2421from sage.structure.sage_object import SageObject
    2522from sage.calculus.functional import diff
    2623from sage.functions.other import sqrt
     
    4542
    4643class ParametrizedSurface3D(GenericManifold):
    4744    r"""
    48     This is a class for finding basic invariants of surfaces.
     45    This class includes methods for calculation the main geometrical objects related
     46    to a parametrized surface such as the first and the second fundamental form, the total (gaussian)
     47    and the mean curvature, the geodesic curves, the parallel transport, etc.
    4948
    5049    The orientation of the surface is determined by the parametrization, that is, the
    5150    natural frame with positive orientation is given by $\partial_1 \vec r$, $\partial_2 \vec r$.         
    5251
    53     jvkersch: This description, as well as the doctests, should be expanded.
    5452
    5553    INPUT:
    5654
     
    6058
    6159     - ``name_of_surface`` - string with the name of the surface (optional).
    6260       
    63     EXAMPLES::
     61    EXAMPLES:
    6462
    65       sage: u, v = var('u,v')
    66       sage: paraboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'paraboloid'); paraboloid
    67       Parametrized surface ('paraboloid') with equation [u, v, u^2 + v^2]
     63      Let us examplify the class by solving several standard tasks of differential geometry of surfaces.
     64      First, look how to set a surface.
     65     
     66      The elliptic paraboloid::
     67           
     68            sage: u, v = var('u,v')
     69            sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
     70            sage: eparaboloid
     71            Parametrized surface ('elliptic paraboloid') with equation [u, v, u^2 + v^2]
     72     
     73      An ellipsoid with axes $a$, $b$, $c$::
    6874
    69     """
     75            sage: a, b,c  = var('a,b,c');
     76            sage: u1, u2 = var ('u1,u2');
     77            sage: u = [u1,u2]
     78            sage: ellipsoid_parametric_equation = vector([a*cos(u1)*cos(u2),b*sin(u1)*cos(u2),c*sin(u2)])
     79            sage: ellipsoid = ParametrizedSurface3D(ellipsoid_parametric_equation,u,'ellipsoid')
     80            sage: print ellipsoid
     81            Parametrized surface ('ellipsoid') with equation (a*cos(u1)*cos(u2), b*sin(u1)*cos(u2), c*sin(u2))
     82     
     83      One may obtain the latex output of the equation of the surface::
    7084
     85            sage: u, v = var('u, v')
     86            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     87            sage: print latex(sphere)
     88            \left[\cos\left(u\right) \cos\left(v\right), \sin\left(u\right) \cos\left(v\right), \sin\left(v\right)\right]
     89            sage: print sphere._latex_()
     90            \left[\cos\left(u\right) \cos\left(v\right), \sin\left(u\right) \cos\left(v\right), \sin\left(v\right)\right]
     91            sage: print sphere
     92            Parametrized surface ('sphere') with equation [cos(u)*cos(v), sin(u)*cos(v), sin(v)]
     93     
     94      We can plot the surface using the inner method plot::
     95           
     96            sage: ellipsoid_parametric_equation_abc = ellipsoid_parametric_equation.substitute(a=2,b=1.5,c=1)
     97            sage: ellipsoid_abc = ParametrizedSurface3D(ellipsoid_parametric_equation_abc,[u1,u2],'ellipsoid_abc')
     98            sage: ellipsoid_abc_plot = ellipsoid_abc.plot((u1,0,2*pi),(u2,-pi/2,pi/2))
     99            sage: ellipsoid_abc_plot.show(aspect_ratio=(1,1,1))
     100
     101      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::
     102
     103            sage: ellipsoid.natural_frame()
     104            {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))}     
     105
     106      Now find the normal vector field.   
     107      The normal vector field which is the vector product of the vectors of the natural frame is::
     108
     109            sage: ellipsoid.normal_vector()
     110            (b*c*cos(u1)*cos(u2)^2, a*c*sin(u1)*cos(u2)^2, a*b*sin(u2)*cos(u2))
     111
     112      And the unit normal vector field of the elliptic paraboloid is::
     113     
     114            sage: u, v = var('u,v')
     115            sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
     116            sage: eparaboloid.normal_vector(1)
     117            (-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)) 
     118     
     119      Now let us find the coefficients of the first  fundamental form of the torus::
     120
     121            u, v = var('u,v')
     122            a, b = var('a,b')
     123            sage: torus = ParametrizedSurface3D(((a + b*cos(u))*cos(v),(a + b*cos(u))*sin(v), b*sin(u)),[u,v],'torus')
     124            sage: torus.first_fundamental_form_coefficients()
     125            {(1, 2): 0, (1, 1): b^2, (2, 1): 0, (2, 2): b^2*cos(u)^2 + 2*a*b*cos(u) + a^2}
     126     
     127      We can find the length of a curve on the surface.
     128      For example, let us find the length of the curve $u^1 = t$, $u^2 =t$, $t \in [0,2\pi]$,
     129      on the ellipsoid with axes $a=1$, $b=1.5$ and $c=1$. So we take the curve::
     130 
     131            sage: t = var('t')
     132            sage: uu1 = t
     133            sage: uu2 = t
     134
     135      Then find the tangent vector::
     136
     137            sage: duu1 = diff(uu1,t)
     138            sage: duu2 = diff(uu2,t)
     139            sage: duu = vector([duu1,duu2])
     140
     141      Then calculate the symbolic expression for the length::
     142 
     143            sage: L=sqrt(ellipsoid.first_fundamental_form(duu,duu).substitute(u1=uu1,u2=uu2))
     144            sage: integrate(L,t,0,2*pi)
     145            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)
     146
     147      And finally find the numerical value::
     148
     149            sage: print numerical_integral(L.substitute(a=2,b=1.5,c=1),0,1)[0]
     150            2.00127905972
     151
     152      Find the area of the sphere of radius $R$::
     153
     154            sage: R = var('R');
     155            sage: u, v = var('u,v');
     156            sage: assume(R>0)
     157            sage: assume(cos(v)>0)
     158            sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     159            sage: integral(integral(sphere.area_form(),u,0,2*pi),v,-pi/2,pi/2)
     160            4*pi*R^2
     161
     162      We can find an orthonormal frame field $\{e_1, e_2\}$ of a surface and calculate it structure functions $c^k_{ij}$,
     163      where $[e_i,e_j] = \sum_k c^k_{ij} e_k$.
     164      Let us first find the orthonormal frame field on the elliptic paraboloid::
     165
     166            sage: u, v = var('u,v')
     167            sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
     168            sage: eparaboloid.orthonormal_frame()
     169            {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)))}   
     170     
     171
     172      Now find the inner (local) coordinates of the orthonormal frame field on the elliptic paraboloid::
     173           
     174            sage: EE = eparaboloid.orthonormal_frame(coordinates='int')
     175            sage: E1 = EE[1]; E2 = EE[2]
     176            sage: CC = eparaboloid.frame_structure_functions(E1,E2)
     177            sage: CC[1,2,1].simplify_full()
     178            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))
     179
     180
     181      We find the Gaussian and Mean Curvatures of the sphere::
     182
     183            sage: u, v = var('u,v')
     184            sage: sphere = ParametrizedSurface3D((u,v,sqrt(1-u^2-v^2)),[u,v],'sphere')
     185            sage: sphere.gauss_curvature()
     186            1
     187            sage: sphere.mean_curvature()
     188            1
     189
     190      Another example is to plot gaussian curvature of ellipsoid by color::
     191
     192            sage: u1, u2 = var('u1,u2');
     193            sage: u = [u1,u2]
     194            sage: ellipsoid_equation(u1,u2) = [2*cos(u1)*cos(u2),1.5*cos(u1)*sin(u2),sin(u1)]
     195            sage: ellipsoid = ParametrizedSurface3D(ellipsoid_equation(u1,u2),u,'ellipsoid')
     196            sage: # set intervals for variables and the number of division points
     197            sage: interval_1 = (-1.5,1.5)
     198            sage: interval_2 = (0,6.28)
     199            sage: plot_points = [10,20]
     200            sage: # make the arguments array
     201            sage: u1_array = [interval_1[0] + cnt*(interval_1[1]-interval_1[0])/plot_points[0] for cnt in range(0,plot_points[0])]
     202            sage: u2_array = [interval_2[0] + cnt*(interval_2[1]-interval_2[0])/plot_points[1] for cnt in range(0,plot_points[1])]
     203            sage: u_array = [ (uu1,uu2) for uu1 in u1_array for uu2 in u2_array]
     204            sage: # Find the gaussian curvature
     205            sage: K(u1,u2) = ellipsoid.gauss_curvature()
     206            sage: # Make array of K values
     207            sage: K_array = [K(uu[0],uu[1]) for uu in u_array]
     208            sage: # Find minimum and max of the gauss curvature
     209            sage: K_max = max(K_array)
     210            sage: K_min = min(K_array)
     211            sage: # Make the array of color coefficients
     212            sage: cc_array = [ (ccc - K_min)/(K_max - K_min) for ccc in K_array ]
     213            sage: points_array = [ellipsoid_equation(u_array[counter][0],u_array[counter][1]) for counter in range(0,len(u_array)) ]
     214            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)) )
     215            sage: curvature_ellipsoid_plot.show(aspect_ratio=1)
     216
     217      We can find the principal curvatures and principal directions of the elliptic paraboloid::
     218 
     219            sage: u, v = var('u,v')
     220            sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
     221            sage: pd = eparaboloid.principal_directions()
     222
     223      We extract the principal curvatures::
     224
     225            sage: k1 = pd[0][0].simplify_full()
     226            sage: k1
     227            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)
     228            sage: k2 = pd[1][0].simplify_full()
     229            sage: k2
     230            2/sqrt(4*u^2 + 4*v^2 + 1)
     231   
     232      and let us check it::
     233
     234            sage: K = eparaboloid.gauss_curvature().simplify_full()
     235            sage: K
     236            4/(16*u^4 + 8*u^2 + 16*v^4 + 8*(4*u^2 + 1)*v^2 + 1) 
     237            sage: H = eparaboloid.mean_curvature().simplify_full()
     238            sage: H
     239            2*(2*u^2 + 2*v^2 + 1)/(4*u^2 + 4*v^2 + 1)^(3/2)
     240            sage: (K - k1*k2).simplify_full()
     241            0
     242            sage: (2*H - k1 - k2).simplify_full()
     243            0
     244
     245      And we can find the intrinsic (local coordinates) of the principal directions::
     246         
     247           sage: pd[0][1]
     248           [(1, v/u)]
     249           sage: pd[1][1]
     250           [(1, -u/v)]
     251           
     252      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).
     253          """
     254     
    71255    def __init__(self,equation,variables,*name):
    72256        """
    73257        See ``ParametrizedSurface3D`` for full documentation.
     
    136320               {'designation': name, 'eq': str(self.equation)}
    137321
    138322
    139     def plot3d(self, urange, vrange):
     323    def plot(self, urange, vrange):
    140324        """
    141325        Enable easy plotting directly from the surface class.
    142326
     
    150334            sage: u, v = var('u, v')
    151335            sage: eq = [3*u + 3*u*v^2 - u^3, 3*v + 3*u^2*v - v^3, 3*(u^2-v^2)]
    152336            sage: enneper = ParametrizedSurface3D(eq,[u,v],'enneper_surface')
    153             sage: enneper.plot3d((u, -5, 5), (v, -5, 5))
     337            sage: enneper.plot((u, -5, 5), (v, -5, 5))
    154338        """
    155339        from sage.plot.plot3d.parametric_plot3d import parametric_plot3d
    156340        P = parametric_plot3d(self.equation, urange, vrange)
     
    238422        Compute a single component $g_{ij}$ of the first fundamental form.  If the
    239423        parametric representation of the surface is given by the vector function
    240424        $\vec r(u^i)$, where $u^i$, $i = 1, 2$ are curvilinear coordinates, then
    241 
    242          .. math::
    243 
    244            g_{ij} = \frac{\partial \vec r}{\partial u^i} \cdot
    245              \frac{\partial \vec r}{\partial u^j}.
     425        $g_{ij} = \frac{\partial \vec r}{\partial u^i} \cdot \frac{\partial \vec r}{\partial u^j}$.
    246426
    247427        INPUT:
    248428
     
    343523        Returns the square of the coefficient of the area form on the surface.
    344524        In terms of the coefficients $g_{ij}$ (where $i, j = 1, 2$) of the first fundamental
    345525        form, this invariant is given by
    346 
    347         .. math::
    348 
    349            A^2 = g_{11}g_{22} - g_{12}^2.
     526        $A^2 = g_{11}g_{22} - g_{12}^2$.
    350527
    351528        See also ``area_form``.
    352529       
     
    372549        Returns the coefficient of the area form on the surface.  In terms of the
    373550        coefficients $g_{ij}$ (where $i, j = 1, 2$) of the first fundamental form, the
    374551        coefficient of the area form is given by
    375 
    376         .. math::
    377 
    378            A = \sqrt{g_{11}g_{22} - g_{12}^2}.
     552        $A = \sqrt{g_{11}g_{22} - g_{12}^2}$.
    379553
    380554        See also ``area_form_squared``.
    381555       
     
    399573        r"""
    400574        Returns the coefficients $g^{ij}$ of the inverse of the fundamental form,
    401575        as a dictionary.  The inverse coefficients are defined by
    402 
    403          .. math::
    404 
    405                \sum_j g^{ij} g_{jk} = \delta^i_k
    406 
     576        $\sum_j g^{ij} g_{jk} = \delta^i_k$
    407577        with $\delta^i_k$ the Kronecker delta.
    408578
    409579        OUTPUT:
     
    650820        Returns the structure functions $c^k_{ij}$ for a frame field $e_1, e_2$, i.e.
    651821        a pair of vector fields on the surface which are linearly independent at
    652822        each point.  The structure functions are defined using the Lie bracket by
    653 
    654          .. math::
    655 
    656               [e_i,e_j] = c^k_{ij}e_k.
     823        $[e_i,e_j] = c^k_{ij}e_k$.
    657824         
    658825        INPUT:
    659826
     
    801968        r"""
    802969        Returns the coefficient $h_{ij}$ of the second fundamental form corresponding to the
    803970        index $(i, j)$.  If the equation of the surface is $\vec{r}(u^1, u^2)$, then
    804 
    805          .. math::
    806 
    807             h_{ij} = \vec{r}_{u^i u^j} \cdot \vec{n},
    808 
     971        $h_{ij} = \vec{r}_{u^i u^j} \cdot \vec{n}$,
    809972        where $\vec{n}$ is the unit normal.
    810973
    811974        INPUT:
     
    9541117        return (enum/denom).simplify_full()
    9551118
    9561119
    957     @cached_method
    958     def principal_curvatures(self):
    959         """
    960         Finds the principal curvatures of the surface
    961 
    962         OUTPUT:
    963         The two principal curvatures, given as a dictionary.
    964 
    965         EXAMPLES::
    966                      
    967            sage: R = var('R')
    968            sage: assume(R>0)
    969            sage: u, v = var('u,v')
    970            sage: assume(cos(v)>0)
    971            sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
    972            sage: sphere.principal_curvatures()
    973            {1: -1/R, 2: -1/R}
    974 
    975            sage: u, v = var('u,v')
    976            sage: R, r = var('R, r')
    977            sage: assume(R>r,r>0)
    978            sage: torus = ParametrizedSurface3D([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
    979            sage: torus.principal_curvatures()
    980            {1: -cos(v)/(r*cos(v) + R), 2: -1/r}
    981 
    982         """
    983 
    984         from sage.symbolic.assumptions import assume
    985         from sage.symbolic.relation import solve
    986         from sage.calculus.var import var
    987 
    988         KK = self.gauss_curvature()
    989         HH = self.mean_curvature()
    990 
    991         # jvkersch: when this assumption is uncommented, Sage raises an error stating that the assumption
    992         # is redundant... Can we safely omit this, based on some geometric reasoning?
    993        
    994         # assume(HH**2-KK>=0)
    995         # mikarm: This is a problem, I had a lot of trobles here. Of course, this inequality always hold true.
    996         # I insert this assumption because sage sometimes, in simplification, uses the complex numbers, though the roots
    997         # are, for sure, real. 
    998         # This, in turn, causes problems when we substitute coordinates into the expression of principal curvatures.
    999         # Unfortunately, I did not manage to tell Sage that they are real (to declare the variables as real).
    1000         # Moreover, in a neighborhood of an umbilic point we even cannot "smoothly" order the set of principal curvatures.
    1001         # So, in general, at present this method is far from the final form.
    1002          
    1003 
    1004  
    1005         x = var('x')
    1006         sol = solve(x**2 -2*HH*x + KK == 0,x)
    1007        
    1008         #k1=var('k1')
    1009         #k2=var('k2')
    1010 
    1011         # 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.
    1012 
    1013         k1 = (x.substitute(sol[0])).simplify_full()
    1014         if len(sol) == 1:
    1015             k2 = k1
    1016         else:
    1017             k2 = (x.substitute(sol[1])).simplify_full()
    1018        
    1019         return {1:k1,2: k2}
    1020 
    10211120
    10221121    @cached_method
    10231122    def shape_operator_coefficients(self):
     
    10531152
    10541153
    10551154
    1056     # jvkersch: I've changed the definition of the shape operator to return the matrix
    1057     # of the shape operator, rather than a matrix times a vector.  This will make it easier
    1058     # to compute eigenvalues and eigenvectors of the shape operator.
    10591155   
    10601156    def shape_operator(self):
    10611157        r"""
     
    10871183            [       -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)]
    10881184            sage: S.eigenvalues()
    10891185            [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)]
    1090             sage: paraboloid.principal_curvatures()
    1091             {1: 2/(4*u^2 + 4*v^2 + 1)^(3/2), 2: 2/sqrt(4*u^2 + 4*v^2 + 1)}
    10921186
    10931187        """
    10941188       
    1095         #vv = vector([xx for xx in v])
    10961189        shop = self.shape_operator_coefficients()
    10971190        shop_matrix=matrix([[shop[(1,1)],shop[(1,2)]],[shop[(2,1)],shop[(2,2)]]])
    10981191        return simplify_matrix(shop_matrix)
    10991192
    11001193
    1101     @cached_method
     1194
     1195
     1196
    11021197    def principal_directions(self):
    1103         """
    1104         Finds the principal curvatures and principal directions of the surface
    1105 
    1106         OUTPUT:
    1107 
    1108          - The dictionary of lists [principal curvature, corresponding principal direction]         
    1109 
    1110         If principal curvatures coincide, gives the warning that the surface is a sphere.
    1111 
    1112         EXAMPLES::
    1113 
    1114            sage: u, v = var('u, v')
    1115            sage: R, r = var('R,r')
    1116            sage: assume(R>r,r>0)
    1117            sage: torus = ParametrizedSurface3D([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
    1118            sage: pdd = torus.principal_directions()
    1119            sage: pdd[1]
    1120            [-cos(v)/(r*cos(v) + R), (1, 0)]
    1121            sage: pdd[2]
    1122            [-1/r, (0, -(R*r*cos(v) + R^2)/r)]
    1123 
    1124            sage: R = var('R')
    1125            sage: assume(R>0)
    1126            sage: assume(cos(v)>0)
    1127            sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
    1128            sage: sphere.principal_directions()
    1129            'This is a sphere, so any direction is principal'
    1130 
    1131            sage: catenoid = ParametrizedSurface3D([R*cosh(v)*cos(u), R*cosh(v)*sin(u),v],[u,v],'catenoid')
    1132            sage: pd = catenoid.principal_directions()
    1133            sage: pd[1][1]
    1134            (1, 0)
    1135            sage: pd[2][1]
    1136            (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))
    1137            sage: pd[1][1]*pd[2][1]
    1138            0
    1139         """
    1140         gg = self.first_fundamental_form_coefficients()
    1141         hh = self.second_fundamental_form_coefficients()
    1142         kk = self.principal_curvatures()
    1143         if kk[1]==kk[2]:
    1144             return "This is a sphere, so any direction is principal"
    1145         pd1 = simplify_vector([hh[(1,2)]-kk[1]*gg[(1,2)],-hh[(1,1)]+kk[1]*gg[(1,1)]])
    1146         if pd1==vector([0,0]):
    1147            pd1 = vector([1,0])
    1148         pd2 = simplify_vector([hh[(1,2)]-kk[2]*gg[(1,2)],-hh[(1,1)]+kk[2]*gg[(1,1)]])
    1149         if pd2==vector([0,0]):
    1150            pd2 = vector([1,0])
    1151         return {1:[kk[1],pd1],2:[kk[2],pd2]}
    1152 
    1153 
    1154     ### Jvkersch: doctests checked up to this point.
    1155 
    1156     # jvkersch: alternative definition of principal_directions.  The behavior of this
    1157     # function is consistent with the eigenvector  function in Sage (since it is just
    1158     # looking for eigenvectors of the shape operator).  In particular, for a spherical
    1159     # surface it just returns an arbitrary pair of vectors.
    1160     def principal_directions_new(self):
    11611198        r"""
    11621199        Finds the principal curvatures and principal directions of the surface
    11631200
     
    11731210           sage: R, r = var('R,r')
    11741211           sage: assume(R>r,r>0)
    11751212           sage: torus = ParametrizedSurface3D([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
    1176            sage: torus.principal_directions_new()
     1213           sage: torus.principal_directions()
    11771214           [(-cos(v)/(r*cos(v) + R), [(1, 0)], 1), (-1/r, [(0, 1)], 1)]
    11781215
    11791216        """
     
    11851222        r"""
    11861223        Computes the connection coefficients or Christoffel symbols $\Gamma^k_{ij}$ of the surface.
    11871224        If the coefficients of the first fundamental form are given by $g_{ij}$
    1188         (where $i, j = 1, 2$), then
    1189        
    1190          .. math::
    1191 
    1192            \Gamma^k_{ij} = \frac{1}{2} g^{kl} \left( \frac{\partial g_{li}}{\partial x^j}
    1193              - \frac{\partial g_{ij}}{\partial x^l}
    1194              + \frac{\partial g_{lj}}{\partial x^i} \right).
    1195 
     1225        (where $i, j = 1, 2$), then
     1226        $\Gamma^k_{ij} = \frac{1}{2} g^{kl} \left( \frac{\partial g_{li}}{\partial x^j}
     1227        - \frac{\partial g_{ij}}{\partial x^l}
     1228        + \frac{\partial g_{lj}}{\partial x^i} \right)$.
    11961229        Here, $(g^{kl})$ is the inverse of the matrix $(g_{ij})$, with $i, j = 1, 2$.
    11971230
    11981231        OUTPUT:
     
    12631296    def geodesics_numerical(self, p0, v0, tinterval):
    12641297        r"""
    12651298        Numerical integration of the geodesic equations.  Explicitly, the
    1266         geodesic equations are given by
    1267 
    1268          .. math::
    1269 
    1270            \frac{d^2 u^i}{dt^2} + \Gamma^i_{jk} \frac{d u^j}{dt} \frac{d u^k}{dt} = 0.
     1299        geodesic equations are given by
     1300        $\frac{d^2 u^i}{dt^2} + \Gamma^i_{jk} \frac{d u^j}{dt} \frac{d u^k}{dt} = 0$.
    12711301
    12721302        Solving these equations gives the coordinates $(u^1, u^2)$ on the surface
    12731303        of the geodesic.  The coordinates in space can then be found by substituting
     
    13061336           sage: assume(cos(q)>0)
    13071337           sage: sphere = ParametrizedSurface3D([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
    13081338           sage: gg_array = sphere.geodesics_numerical([0.0,0.0],[1.0,1.0],[0,2*pi,5])
    1309            sage: gg_array[0]
    1310            [0, [0.000000000000000, 0.000000000000000], [1.00000000000000, 1.00000000000000], [1, 0, 0]]
    1311            sage: gg_array[1]
    1312            [1.2566370614359172, [0.76440092815484362, 1.8586224292405702], [-0.28386842687533731, 1.9194187166389509], [-0.204895409519981, 0.692104714174527, 0.692104714458033]]
     1339           sage: [N(xx[0],digits=4) for xx in gg_array]
     1340           [0.0000, 1.257, 2.513, 3.770, 5.027, 6.283] 
     1341           sage: [ [ N(x,digits=4) for x in xx[1] ] for xx in gg_array ]
     1342           [[0.0000, 0.0000], [0.7644, 1.859], [-0.2876, 3.442], [-0.6137, 5.502], [0.5464, 6.937], [0.3714, 9.025]]   
     1343           sage: [ [ N(x,digits=4) for x in xx[3] ] for xx in gg_array ]
     1344           [[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]] 
    13131345        """
    13141346
    13151347        u1 = self.variables[1]
     
    13781410        r"""
    13791411        Numerically solves the equations for parallel translation of a vector along a curve
    13801412        on the surface.  Explicitly, the equations for parallel translation are given by
    1381 
    1382          .. math::
    1383 
    1384            \frac{d u^i}{dt} + u^j \frac{d c^k}{dt} \Gamma^i_{jk} = 0,
    1385 
     1413        $\frac{d u^i}{dt} + u^j \frac{d c^k}{dt} \Gamma^i_{jk} = 0$,
    13861414        where $\Gamma^i_{jk}$ are the connection coefficients of the surface, the vector to
    13871415        be transported has components $u^j$ and the curve along which to transport has
    13881416        components $c^k$.