Ticket #10132: trac_10132_differential_geometry_sage_v3.patch

sage/all.py

@@ -139 +139 @@
 from sage.ext.fast_eval      import fast_float
 
 from sage.tensor.all     import *
+from sage.riemann.all    import *
 
 from copy import copy, deepcopy
 

new file sage/riemann/__init__.py

@@ +1 @@
+# This comment is here so the file is non-empty (so Mercurial will check it in).

new file sage/riemann/all.py

@@ +1 @@
+from parametrized_surface3d import ParametrizedSurface3D

new file sage/riemann/parametrized_surface3d.py

@@ +1 @@
+"""
+Basic invariants of parametrized surfaces.
+
+
+AUTHORS: 
+        - Mikhail Malakhaltsev (2010-09-25): initial version
+
+"""    
+#*****************************************************************************
+#       Copyright (C) 2010  Mikhail Malakhaltsev <mikarm@gmail.com>
+#
+#  Distributed under the terms of the GNU General Public License (GPL)
+#                  http://www.gnu.org/licenses/
+#*****************************************************************************
+
+
+# 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
+
+from sage.misc.cachefunc import cached_method
+
+
+
+# TODO (distant future): update class...
+class GenericManifold(SageObject):
+    pass
+
+
+class ParametrizedSurface3D(GenericManifold):
+    """
+    This is a class for finding basic invariants of surfaces.
+
+    USAGE:
+    parametrized_surface3d(surface_equation,variables,name_of_surface)
+    where
+      surface_equation is a vector or list of three components
+      variables is a list of two variables
+      name_of_surface is a string (optional)
+
+    Warning: The orientation on the surface is given by the parametrization, that is the positive frame is 
+    $\partial_1 \\vec r$, $\partial_2 \\vec r$.
+        
+    This class includes the following methods:
+      natural_frame 
+      normal_vector
+      first_fundamental_form_coefficients
+      first_fundamental_form
+      first_fundamental_form_inverse_coefficients
+      first_fundamental_form_inverse_coefficients
+      area_form_squared
+      area_form
+      rotation
+      orthonormal_frame
+      lie_bracket
+      frame_structure_functions
+      second_order_natural_frame
+      second_fundamental_form_coefficients
+      second_fundamental_form
+      gauss_curvature
+      mean_curvature
+      principal curvatures
+      principal directions
+      connection_coefficients
+      geodesics_numerical
+      parallel_translation_numerical
+      
+
+    EXAMPLES::
+
+      sage: var('u,v')
+      sage: paraboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'paraboloid')
+      sage: paraboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v])
+           
+    """
+
+    def __init__(self,equation,variables,*name):
+
+        self.equation = equation 
+        self.variables_list = variables
+        self.variables = {1:self.variables_list[0],2:self.variables_list[1]}
+        self.name = name
+
+        # Callable version of the underlying equation
+        def eq_callable(u, v):
+            u1, u2 = self.variables_list
+            return [fun.subs({u1: u, u2: v}) for fun in self.equation]
+            
+        self.eq_callable = eq_callable
+
+        # Various index tuples        
+        self.index_list = [(i,j) for i in (1,2) for j in (1,2)]
+        self.index_list_3=[(i,j,k) for i in (1,2) for j in (1,2) for k in (1,2)]
+
+
+    @cached_method            
+    def natural_frame(self,vector_number=0):
+        """ 
+        This functions find the natural frame of the parametrized surface
+        INPUT: 
+        the argument can be empty, or equal to 1, or to 2 
+
+        OUTPUT: 
+        if argument is equal to 1 or 2, the output is the corresponding vector of the natural frame 
+        if argument is empty, the output is the dictionary of the vector fields of the natural frame
+
+        EXAMPLES::  
+
+            sage:  eparaboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
+            sage:  eparaboloid.natural_frame()
+            {1: (1, 0, 2*u), 2: (0, 1, 2*v)}
+            sage:  eparaboloid.natural_frame(1)
+            (1, 0, 2*u)
+            sage:  eparaboloid.natural_frame(2)
+            (0, 1, 2*v)
+        """ 
+
+        dr1 = vector([diff(f,self.variables[1]).simplify_full() for f in self.equation]) 
+        dr2 = vector([diff(f,self.variables[2]).simplify_full() for f in self.equation])
+
+        return {1:dr1, 2:dr2}
+
+
+    @cached_method
+    def normal_vector(self, normalized=False):
+        """
+        This functions find the normal vector of the parametrized surface
+        
+        INPUT: 
+        empty argument, of a nonzero real number l
+        
+        OUTPUT: 
+        without arguments gives the normal vector which is the vector product of the 
+        natural frame vectors;
+
+        with the argument l gives the normal vector of length l (the orientation determined 
+        by the sign of l) 
+        
+        EXAMPLES::
+
+          sage: eparaboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
+          sage: eparaboloid.normal_vector()
+          (-2*u,-2*v,1)
+          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))
+
+        """
+
+        dr = self.natural_frame()
+        normal = dr[1].cross_product(dr[2])
+
+        if normalized:
+            normal /= normal.norm()
+        return simplify_vector(normal)
+        
+
+
+    # The following can become part of the docstring at some point...
+
+    # I've defined a cached (private) method which does not check its arguments and
+    # computes one component at a time.  Note that caching is just a matter of
+    # prepending "@cached_method" to the method signature.  For this to work, the
+    # arguments of the method must be hashable (e.g tuple, string, ...) and the method
+    # must not change the internal state of the object.  For all of the methods in this
+    # class, this is the case.
+    #
+    # I then added a method that does simple argument checking, converting the arguments to
+    # tuples, and invokes the private method to do the work.  Note that this method should
+    # not be cached, as the user might specify the index of the component as a list, which is
+    # not hashable.
+    #
+    # A second method returns a dictionary of all the components, again invoking the hidden
+    # helper method.
+ 
+    @cached_method
+    def _compute_first_fundamental_form_coefficient(self, index):
+        dr = self.natural_frame()
+        return (dr[index[0]]*dr[index[1]]).simplify_full()        
+    
+    def first_fundamental_form_coefficient(self, index):
+        index = tuple(sorted(index))
+        if index not in self.index_list:
+            raise ValueError, "Index %s out of bounds." % str(index)
+        return self._compute_first_fundamental_form_coefficient(index)
+
+    @cached_method
+    def first_fundamental_form_coefficients(self):
+
+        """
+        This function gives the coefficients of the first fundamental form
+        
+        INPUT:
+        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
+
+        OUTPUT:
+        with empty argument the output is the dictionary of coefficients of the first fundamental form
+
+        with given indices the output is the corresponding coefficient of the first fundamental form, 
+        for example index1 = 1, index2 = 2 gives the coefficient $g_{12}$ 
+
+        EXAMPLES::
+
+            sage: var('u,v')
+            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
+            sage: sphere.first_fundamental_form_coefficients() 
+            {(1, 2): 0, (1, 1): cos(v)^2, (2, 1): 0, (2, 2): 1}
+
+            sage: sphere.first_fundamental_form_coefficients(1,1)
+            cos(v)^2
+
+            sage: sphere.first_fundamental_form_coefficients(2,1)
+            0
+
+            sage: sphere.first_fundamental_form_coefficients(3,2)
+            'The argument is not appropriate. Read doc' 
+        """        
+        coefficients = {}
+        for index in self.index_list:
+            sorted_index = list(sorted(index))
+            coefficients[index] = \
+                self._compute_first_fundamental_form_coefficient(index)
+        return coefficients
+
+
+
+    def first_fundamental_form(self,vector1,vector2):
+        """
+        Finds the value of first fundamental form on two vectors 
+
+        INPUT: 
+        Two vectors $v=(v^1,v^2)$ and $w=(w^1,w^2)$
+
+        OUTPUT:
+        $g_{11} v^1 w^1 + g_{12}(v^1 w^2 + v^2 w^1) + g_{22} v^2 w^2$ 
+
+        #EXAMPLES::
+         
+            sage: var('u,v')
+            sage: var ('v1,v2,w1,w2')
+            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
+
+            sage: sphere.first_fundamental_form(vector([v1,v2]),vector([w1,w2])) 
+            v1*w1*cos(v)^2 + v2*w2
+
+            sage: vv = vector([1,2])
+            sage: sphere.first_fundamental_form(vv,vv)
+            cos(v)^2 + 4
+
+            sage: sphere.first_fundamental_form([1,1],[2,1])
+            2*cos(v)^2 + 1
+        """
+        gg = self.first_fundamental_form_coefficients()
+        return sum(gg[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list)
+
+
+    @cached_method
+    def area_form_squared(self):
+        """
+        Gives $g_{11}g_{22} - g_{12}^2$, which is the square of the coefficient of the area form
+        
+        INPUT: 
+        No arguments
+
+        OUTPUT:  
+        $g_{11}g_{22} - g_{12}^2$
+
+        EXAMPLES::
+
+            sage: var('u,v')
+            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
+            sage: sphere.area_form_squared()
+            cos(v)^2
+
+        """
+        gg = self.first_fundamental_form_coefficients()
+        return  (gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2).simplify_full()
+
+
+    @cached_method
+    def area_form(self):
+        """
+        Gives $\sqrt{g_{11}g_{22} - g_{12}^2}$, which is the coefficient of the area form
+        
+        INPUT: 
+        No arguments
+
+        OUTPUT:  
+        $\sqrt{g_{11}g_{22} - g_{12}^2}$
+
+        EXAMPLES::
+
+            sage: var('u,v')
+            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
+            sage: sphere.area_form_squared()
+            abs(cos(v))
+
+        """
+        return sqrt(self.area_form_squared()).simplify_full()
+
+
+    @cached_method
+    def first_fundamental_form_inverse_coefficients(self):
+        """
+        Gives $g^{ij}$
+        
+        INPUT: 
+        No arguments
+
+        OUTPUT:  
+        dictionary {(1,1):$g^{11}$, (1,2):$g^{12}$, (2,1):$g^{21}$, (2,2):$g^{22}$
+
+        EXAMPLES::
+
+            sage: var('u,v')
+            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
+            sage: sphere.first_fundamental_form_inverse_coefficients()
+            {(1, 2): 0, (1, 1): cos(v)^(-2), (2, 1): 0, (2, 2): 1}
+        """
+
+        gg = self.first_fundamental_form_coefficients()
+        DD = gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2
+        
+        gi11 = (gg[(2,2)]/DD).simplify_full()
+        gi12 = (-gg[(1,2)]/DD).simplify_full()
+        gi21 = gi12
+        gi22 = (gg[(1,1)]/DD).simplify_full()
+
+        return {(1,1):gi11,(1,2):gi12,(2,1):gi21,(2,2):gi22}
+
+    def first_fundamental_form_inverse_coefficient(self, index):
+        index = tuple(sorted(index))
+        if index not in self.index_list:
+            raise ValueError, "Index %s out of bounds." % str(index)
+        return self._compute_first_fundamental_form_inverse_coefficients()[index]
+
+
+    @cached_method
+    def rotation(self,theta):
+        """
+        Gives the matrix of the operator of rotation on the given angle $\\theta$ with respect to the natural frame 
+        
+        INPUT: 
+        given angle $\\theta$
+
+        OUTPUT:  
+        matrix of the operator of rotation on $\\theta$ with respect to the natural frame 
+
+        ALGORITHM:
+        The operator of rotation on $\pi/2$ is $J^i_j = g^{ik}\omega_{jk}$, where $\omega$ is the area form
+        The operator of rotation on angle $\\theta$ is $\cos(\\theta) I + sin(\\theta) J$
+
+        EXAMPLES::
+
+            sage: var('u,v')
+            sage: assume(cos(v)>0)
+            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
+            sage: rotation = sphere.rotation(pi/3)
+            sage: rotation^3
+            [-1  0]
+            [ 0 -1]
+            # it is true because the rotation at $\pi/3$ applied three times gives rotation a $\pi$
+        """
+
+        from sage.functions.trig import sin, cos
+
+        gi = self.first_fundamental_form_inverse_coefficients()
+        w12 = self.area_form() 
+        RR11 = (cos(theta) + sin(theta)*gi[1,2]*w12).simplify_full()
+        RR12 = (- sin(theta)*gi[1,1]*w12).simplify_full()
+        RR21 = (sin(theta)*gi[2,2]*w12).simplify_full()
+        RR22 = (cos(theta) - sin(theta)*gi[2,1]*w12).simplify_full()        
+        rotation = matrix([[RR11,RR12],[RR21,RR22]])
+        return rotation 
+
+
+    # XXXXXX
+
+
+
+
+
+    @cached_method
+    def orthonormal_frame_vector(self, vector_number=0, coordinates='ext'):
+        """
+        Gives the orthonormal frame field of the surface  
+        
+        INPUT: 
+        (vector_number,coordinates), where 
+        vector number is 0 (default), 1 or 2
+        coordinates is 'ext'(default) or 'int' 
+               
+
+        OUTPUT: 
+        If vector_number is equal to 1 or 2, the output is the corresponding vector of an orthonormal frame. 
+        If argument vector_number is empty, the output is the dictionary of the vector fields of an orthonormal frame.
+        
+        If coordinates='ext', output is coordinates in $\mathbb{R}^3$
+        If coordinates='int', output is coordinates with respect to the natural frame
+
+        ALGORITHM:
+        We normalize the first vector $\\vec e_1$ of the natural frame and then get the second frame vector 
+        $\\vec e_2 = [\\vec n, \\vec e_1]$. 
+
+        EXAMPLES::
+
+            sage: var('u,v')
+            sage: assume(cos(v)>0)
+            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
+            sage: EE = sphere.orthonormal_frame()
+            sage: (EE[1]*EE[1]).simplify_full()
+            1
+            sage: (EE[1]*EE[2]).simplify_full()
+            0
+            sage: simplify_vector(vector_product(EE[1],EE[2])-sphere.normal_vector(1))
+            [0,0,0]
+
+            sage: EE_int = sphere.orthonormal_frame(coordinates='int')
+            sage: EE_int
+            {1: (1/cos(v), 0), 2: (0, 1)}
+            sage: sphere.first_fundamental_form(EE_int[1],EE_int[1])
+            1
+            sage: sphere.first_fundamental_form(EE_int[1],EE_int[2])
+            0 
+            sage: sphere.first_fundamental_form(EE_int[2],EE_int[2])
+            1
+            
+        """
+
+        if vector_number not in [1, 2]:
+            raise ValueError, "Basis vector number out of range."
+        return self.orthonormal_frame(coordinates)[vector_number]
+
+
+    @cached_method
+    def orthonormal_frame(self, coordinates='ext'):
+
+        from sage.symbolic.constants import pi
+
+        if coordinates not in ['ext', 'int']:
+            raise ValueError, \
+                  r"Coordinate system must be exterior ('ext') or interior ('int')."
+
+        c  = self.first_fundamental_form_coefficient([1,1])
+        if coordinates == 'ext':
+            f1 = self.natural_frame()[1]
+
+            E1 = simplify_vector(f1/sqrt(c))
+            E2 = simplify_vector(self.normal_vector(normalized=True).cross_product(E1))
+        else:
+            E1 =  vector([(1/sqrt(c)).simplify_full(), 0])
+            E2 = simplify_vector(self.rotation(pi/2)*E1) 
+        return  {1:E1, 2:E2}
+        
+        
+    def lie_bracket(self,v,w):
+        """
+        Gives the Lie bracket of two vector fields which are given by coordinates with respect to the natural frame
+         
+        
+        INPUT:
+        v and w are vectors 
+               
+
+        OUTPUT: 
+        vector [v,w]
+
+
+        EXAMPLES::
+
+            sage: var('u,v')
+            sage: assume(cos(v)>0)
+            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
+            sage: sphere.lie_bracket([u,v],[-v,u])
+            (0, 0)
+
+            sage: EE_int = sphere.orthonormal_frame(coordinates='int')
+            sage: sphere.lie_bracket(EE_int[1],EE_int[2])
+            (sin(v)/cos(v)^2, 0)
+        """
+        vv = vector([xx for xx in v])
+        ww = vector([xx for xx in w])
+        uuu = [self.variables[1],self.variables[2]]
+        Dvv = matrix([ [ diff(xxx,uu).simplify_full() for uu in uuu ] for xxx in vv])   
+        Dww = matrix([ [ diff(xxx,uu).simplify_full() for uu in uuu ] for xxx in ww])   
+        return simplify_vector(vector(Dvv*ww - Dww*vv)) 
+
+
+
+    def frame_structure_functions(self,e1,e2):
+        """
+        Gives the structure functions $c^k_{ij}$ for a frame field $e_1$, $e_2$, where
+        $[e_i,e_j] = c^k_{ij}e_k$ 
+         
+        
+        INPUT:
+        Frame vectors $e_1$, $e_2$ 
+               
+
+        OUTPUT: 
+        The dictionary of $c^k_{ij}$. 
+        Warning: the tuple (i,j,k) corresponds to $c^k_{ij}$ 
+
+
+        EXAMPLES::
+
+            sage: var('u,v')
+            sage: assume(cos(v)>0)
+            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
+            sage: sphere.frame_structure_functions([u,v],[-v,u])
+            {(1, 2, 1): 0, (2, 1, 2): 0, (2, 2, 2): 0, (1, 2, 2): 0, (1, 1, 1): 0, (2, 1, 1): 0, (2, 2, 1): 0, (1, 1, 2): 0}
+            sage: structfun = sphere.frame_structure_functions(EE_int[1],EE_int[2])
+            sage: structfun
+            {(1, 2, 1): sin(v)/cos(v), (2, 1, 2): 0, (2, 2, 2): 0, (1, 2, 2): 0, (1, 1, 1): 0, (2, 1, 1): -sin(v)/cos(v), (2, 2, 1): 0, (1, 1, 2): 0}
+            sage: sphere.lie_bracket(EE_int[1],EE_int[2]) - CC[(1,2,1)]*EE_int[1] - CC[(1,2,2)]*EE_int[2] 
+            (0, 0)
+            """
+        ee1 = vector([xx for xx in e1])
+        ee2 = vector([xx for xx in e2])
+        
+        lb = simplify_vector(self.lie_bracket(ee1,ee2))
+        trans_matrix = matrix([[ee1[0],ee2[0]],[ee1[1],ee2[1]]])
+        zz = simplify_vector(trans_matrix.inverse()*lb)
+        return {(1,1,1):0, (1,1,2):0, (1,2,1):zz[0], (1,2,2):zz[1], (2,1,1):-zz[0], (2,1,2):-zz[1],(2,2,1):0, (2,2,2):0}
+
+
+    @cached_method
+    def _compute_second_order_frame_element(self, index):
+        variables = [self.variables[i] for i in index]
+        ddr_element = vector([diff(f, variables).simplify_full() for f in self.equation])
+        
+        return ddr_element
+
+    @cached_method
+    def second_order_natural_frame(self):
+        
+        """
+        Gives the second derivatives of the equation $\\vec r = \\vec r(u^1,u^2)$ of parametrized surface 
+         
+        
+        INPUT:
+        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
+               
+
+        OUTPUT: 
+        With empty argument the output is the dictionary of $\partial_{ij}\\vec r(u^1,u^2)$.
+
+        With given indices the output is the corresponding second derivative of the surface parametric equation, 
+        For example, index1 = 1, index2 = 2 gives $\partial_{12}\\vec r(u^1,u^2)$.
+
+
+        EXAMPLES::
+
+            sage: var('u,v')
+            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
+            sage: sphere.second_order_natural_frame()
+            {(1, 2): (sin(u)*sin(v), -sin(v)*cos(u), 0), (1, 1): (-cos(u)*cos(v),
+            -sin(u)*cos(v), 0), (2, 1): (sin(u)*sin(v), -sin(v)*cos(u), 0), (2, 2):
+            (-cos(u)*cos(v), -sin(u)*cos(v), -sin(v))}
+            
+            sage: sphere.second_order_natural_frame(1,1)
+            (-cos(u)*cos(v), -sin(u)*cos(v), 0)
+            sage: sphere.second_order_natural_frame(1,2)
+            (sin(u)*sin(v), -sin(v)*cos(u), 0)
+            sage: sphere.second_order_natural_frame(2,2)
+            (-cos(u)*cos(v), -sin(u)*cos(v), -sin(v)
+            """
+
+        vectors = {}
+        for index in self.index_list:
+            sorted_index = tuple(sorted(index))
+            vectors[index] = \
+                self._compute_second_order_frame_element(sorted_index)
+        return vectors
+
+    def second_order_natural_frame_element(self, index):
+        index = tuple(sorted(index))
+        if index not in self.index_list:
+            raise ValueError, "Index %s out of bounds." % str(index)
+        return self._compute_second_order_frame_element(index)
+
+
+    @cached_method
+    def _compute_second_fundamental_form_coefficient(self, index):
+        NN = self.normal_vector(normalized=True)
+        v  = self.second_order_natural_frame_element_new(index) 
+        return (v*NN).simplify_full()
+
+        
+    def second_fundamental_form_coefficient(self, index):
+        index = tuple(index)
+        if index not in self.index_list:
+            raise ValueError, "Index %s out of bounds." % str(index)
+        return self._compute_second_fundamental_form_coefficient(index)
+
+
+    @cached_method
+    def second_fundamental_form_coefficients(self):
+        """
+        This function gives the coefficients $h_{ij}$ of the second fundamental form.
+        
+        INPUT:
+        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
+
+        OUTPUT:
+        with empty argument the output is the dictionary of coefficients of the second fundamental form
+
+        with given indices the output is the corresponding coefficient of the second fundamental form, 
+        for example index1 = 1, index2 = 2 gives the coefficient $h_{12}$ 
+
+        EXAMPLES::
+
+            sage: var('u,v')
+            sage: assume(cos(v)>0)
+            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
+
+            sage: sphere.second_fundamental_form_coefficients()
+            {(1, 2): 0, (1, 1): -cos(v)^2, (2, 1): 0, (2, 2): -1}
+            sage: sphere.second_fundamental_form_coefficients(1,1)
+            -cos(v)^2
+            sage: sphere.second_fundamental_form_coefficients(2,1)
+            0
+
+            sage: sphere.second_fundamental_form_coefficients(3,2)
+            'The argument is not appropriate. Read doc' 
+        """
+
+        coefficients = {}
+        for index in self.index_list:
+            coefficients[index] = \
+                self._compute_second_fundamental_form_coefficient(index)
+        return coefficients
+
+
+    def second_fundamental_form(self,vector1,vector2):
+        """
+        Finds the value of second fundamental form on two vectors 
+
+        INPUT: 
+        Two vectors $v=(v^1,v^2)$ and $w=(w^1,w^2)$
+
+        OUTPUT:
+        $h_{11} v^1 w^1 + h_{12}(v^1 w^2 + v^2 w^1) + h_{22} v^2 w^2$ 
+
+        EXAMPLES::
+
+           sage: var('u,v')
+           sage: var ('v1,v2,w1,w2')
+           sage: assume(cos(v) > 0)
+           sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
+           sage: sphere.second_fundamental_form(vector([v1,v2]),vector([w1,w2]))
+           -v1*w1*cos(v)^2 - v2*w2
+           sage: vv = vector([1,2])
+           sage: sphere.second_fundamental_form(vv,vv)
+           -cos(v)^2 - 4
+           sage: sphere.second_fundamental_form([1,1],[2,1])
+           -2*cos(v)^2 - 1
+                    
+        """
+        hh = self.second_fundamental_form_coefficients()
+        return sum(hh[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list)
+
+
+    @cached_method
+    def gauss_curvature(self):
+        """
+        Finds the gaussian curvature $K = \\frac{h_{11}h_{22} - h_{12}^2}{g_{11}g_{22} - g_{12}^2}$. 
+
+        INPUT: 
+        No arguments 
+
+        OUTPUT:
+        $K = \\frac{h_{11}h_{22} - h_{12}^2}{g_{11}g_{22} - g_{12}^2}$
+
+        EXAMPLES::
+
+           sage: var('R')
+           sage: assume(R>0)
+           sage: var('u,v')
+           sage: assume(cos(v)>0)
+           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
+           sage: sphere.gauss_curvature()
+           R^(-2)
+                    
+        """
+        hh = self.second_fundamental_form_coefficients()
+        return ((hh[(1,1)]*hh[(2,2)]-hh[(1,2)]**2)/self.area_form_squared()).simplify_full()
+
+
+    @cached_method
+    def mean_curvature(self):
+        """
+        Finds the mean curvature $H = \\frac{1}{2}\\frac{g_{22}h_{11} - 2g_{12}h_{12} + g_{11}h_{22}}{g_{11}g_{22} - g_{12}^2}$. 
+
+        INPUT: 
+        No arguments 
+
+        OUTPUT:
+        $H = \\frac{1}{2}\\frac{g_{22}h_{11} - 2g_{12}h_{12} + g_{11}h_{22}}{g_{11}g_{22} - g_{12}^2}$. 
+
+        EXAMPLES::
+           
+           sage: var('R')
+           sage: assume(R>0)
+           sage: var('u,v')
+           sage: assume(cos(v)>0)
+           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
+           sage: sphere.mean_curvature()
+           -1/R
+
+        """
+        gg = self.first_fundamental_form_coefficients()
+        hh = self.second_fundamental_form_coefficients()
+        denom = 2*self.area_form_squared()
+        enum =(gg[(2,2)]*hh[(1,1)]-2*gg[(1,2)]*hh[(1,2)]+gg[(1,1)]*hh[(2,2)]).simplify_full()
+        return (enum/denom).simplify_full()
+
+
+    @cached_method
+    def principal_curvatures(self):
+        """
+        Finds the principal curvatures of the surface 
+
+        INPUT: 
+        No arguments 
+
+        OUTPUT:
+        The dictionary of principal curvatures         
+
+        EXAMPLES::
+                      
+           sage: var('R')
+           sage: assume(R>0)
+           sage: var('u,v')
+           sage: assume(cos(v)>0)
+           sage: sphere = parametrized_surface3d([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: var('u,v')
+           sage: var('R,r')
+           sage: assume(R>r,r>0)
+           sage: torus = parametrized_surface3d([(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 principal_directions(self):
+        """
+        Finds the principal curvatures and principal directions of the surface 
+
+        INPUT: 
+        No arguments 
+
+        OUTPUT:
+        The dictionary of lists [a principal curvature, the corresponding principal direction]         
+
+        If principal curvatures coincide, gives the warning that the surface is a sphere.
+
+        EXAMPLES::
+
+           sage: var('u,v')
+           sage: var('R,r')
+           sage: assume(R>r,r>0)
+           sage: torus = parametrized_surface3d([(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: var('RR')
+           sage: assume(RR>0)
+           sage: var('u,v')
+           sage: assume(cos(v)>0)
+           sage: sphere = parametrized_surface3d([RR*cos(u)*cos(v),RR*sin(u)*cos(v),RR*sin(v)],[u,v],'sphere')
+           sage: sphere.principal_directions()
+           'This is a sphere, so any direction is principal'
+
+           sage: var('aa')
+           sage: assume(aa>0)
+           sage: catenoid = parametrized_surface3d([aa*cosh(v)*cos(u),aa*cosh(v)*sin(u),v],[u,v],'catenoid')
+           sage: pd = catenoid.principal_directions()
+           sage: pd[1][1]
+           (0, 1/2*(2*aa^3*sinh(v)^2*cosh(v) + sqrt(4*aa^4*sinh(v)^2*cosh(v)^2 + aa^4 + 4*aa^2*cosh(v)^2 - 2*aa^2 + 1)*aa*cosh(v) + (aa^3 + aa)*cosh(v))/(aa^2*sinh(v)^2 + 1)^(3/2))
+           sage: pd[2][1]
+           (1, 0)
+           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]}
+
+
+    @cached_method
+    def connection_coefficients(self):
+        """
+        Finds the connection coefficients of the surface 
+
+        INPUT: 
+        No arguments 
+
+        OUTPUT:
+        The dictionary of connection coefficients.
+
+        Warning: the triple $(i,j,k)$ corresponds to $\Gamma^k_{ij}$.  
+
+        EXAMPLES::
+                                 
+           sage: var('r')
+           sage: assume(r > 0)
+           sage: var('u,v')
+           sage: assume(cos(v)>0)
+           sage: sphere = parametrized_surface3d([r*cos(u)*cos(v),r*sin(u)*cos(v),r*sin(v)],[u,v],'sphere')
+           sage: sphere.connection_coefficients()
+           {(1, 2, 1): -sin(v)/cos(v), (2, 2, 2): 0, (1, 2, 2): 0, (2, 1, 1): -sin(v)/cos(v), (1, 1, 2): sin(v)*cos(v), (2, 2, 1): 0, (2, 1, 2): 0, (1, 1, 1): 0}
+           
+        """
+        x = self.variables
+        gg = self.first_fundamental_form_coefficients()
+        gi = self.first_fundamental_form_inverse_coefficients()
+
+        dg = {}
+        for kkk in self.index_list_3:
+            dg[kkk]=gg[(kkk[1],kkk[2])].differentiate(x[kkk[0]]).simplify_full()
+        structfun={}
+
+        for kkk in self.index_list_3:
+            structfun[kkk]=sum(gi[(kkk[2],s)]*(dg[(kkk[0],kkk[1],s)]+dg[(kkk[1],kkk[0],s)]-dg[(s,kkk[0],kkk[1])])/2 for s in (1,2)).full_simplify()
+        return structfun
+
+
+    # jvkersch: this private method creates an ode_solver object, which can be used
+    # to integrate the geodesic equations numerically.
+    
+    @cached_method
+    def _create_geodesic_ode_system(self):
+        from sage.ext.fast_eval import fast_float
+        from sage.calculus.var import var
+        from sage.gsl.ode import ode_solver
+
+        u1 = self.variables[1]
+        u2 = self.variables[2]
+        v1, v2 = var('v1, v2')
+
+        C = self.connection_coefficients()
+        
+        dv1 = - C[(1,1,1)]*v1**2 - 2*C[(1,2,1)]*v1*v2 - C[(2,2,1)]*v2**2
+        dv2 = - C[(1,1,2)]*v1**2 - 2*C[(1,2,2)]*v1*v2 - C[(2,2,2)]*v2**2
+        fun1 = fast_float(dv1, str(u1), str(u2), str(v1), str(v2))
+        fun2 = fast_float(dv2, str(u1), str(u2), str(v1), str(v2))
+
+        geodesic_ode = ode_solver()
+        geodesic_ode.function = \
+                              lambda t, (u1, u2, v1, v2) : \
+                              [v1, v2, fun1(u1, u2, v1, v2), fun2(u1, u2, v1, v2)]
+        return geodesic_ode
+
+
+    # jvkersch: integrate the geodesic equations numerically
+    # mikarm: Very good. I cheked it restarting each time the worksheet, it works much faster.
+        
+    def geodesics_numerical(self, p0, v0, tinterval):
+        """ 
+        This method gives the numerical solution for the geodesic equations  
+
+        INPUT:
+        p0 is the list of the coordinates of the initial point 
+        
+        v0 is the  list of the coordinates of the initial vector 
+
+        tinterval is the list [a,b,M], where (a,b) is the domain of the geodesic, M is the number of division points
+
+        OUTPUT:
+        The list consisting of lists [t, [u1(t),u2(t)], [v1(t),v2(t)], [x1(t),x2(t),x3(t)]], where
+
+        t is a subdivision point
+
+        [u1(t),u2(t)] is the list of coordinates of the geodesic point 
+         
+        [v1(t),v2(t)] is the list of coordinates of the vector tanget to the geodesic  
+
+        [x1(t),x2(t),x3(t)] is the list of coordinates of the geodesic point in the three-dimensional space
+
+
+        EXAMPLES::
+
+           sage: var('p,q')
+           sage: v = [p,q]
+           sage: assume(cos(q)>0)
+           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
+           sage: gg_array = sphere.geodesics_numerical([0,0],[1,1],[0,2*pi,5])
+           sage: gg_array[0]
+           [0.0, [0.0, 0.0], [1.0, 1.0], (1, 0, 0)]
+           sage: gg_array[1]
+           [1.2566370614359172, [0.76440104189216407, 1.8586223516062499], [-0.2838683571264714, 1.9194187087799912], (-0.204895333443, 0.692104654602, 0.692104796553)]
+                                            
+         
+        """ 
+
+        u1 = self.variables[1]
+        u2 = self.variables[2]
+
+        solver = self._create_geodesic_ode_system()
+            
+        t_interval, n = tinterval[0:2], tinterval[2]
+        solver.y_0 = [p0[0], p0[1], v0[0], v0[1]]
+        solver.ode_solve(t_span=t_interval, num_points=n)
+
+        parsed_solution = \
+          [[vec[0], vec[1][0:2], vec[1][2:], self.eq_callable(vec[1][0], vec[1][1])] \
+           for vec in solver.solution]
+
+        return parsed_solution
+
+
+    @cached_method
+    def _create_pt_ode_system(self, curve, t):
+        from sage.ext.fast_eval import fast_float
+        from sage.calculus.var import var
+        from sage.gsl.ode import ode_solver
+
+        u1 = self.variables[1]
+        u2 = self.variables[2]
+        v1, v2 = var('v1, v2')
+
+        du1 = diff(curve[0], t)
+        du2 = diff(curve[1], t)
+
+        C = self.connection_coefficients()
+        for coef in C:
+            C[coef] = C[coef].subs({u1: curve[0], u2: curve[1]})
+        
+        dv1 = - C[(1,1,1)]*v1*du1 - C[(1,2,1)]*(du1*v2 + du2*v1) - C[(2,2,1)]*du2*v2
+        dv2 = - C[(1,1,2)]*v1*du1 - C[(1,2,2)]*(du1*v2 + du2*v1) - C[(2,2,2)]*du2*v2
+        fun1 = fast_float(dv1, str(t), str(v1), str(v2))
+        fun2 = fast_float(dv2, str(t), str(v1), str(v2))
+
+        pt_ode = ode_solver()
+        pt_ode.function = lambda t, (v1, v2): [fun1(t, v1, v2), fun2(t, v1, v2)]
+        return pt_ode
+
+
+    # mikarm: We should rewrite it like you did for geodesics.
+    # jvkersch: OK, see this and the above
+    def parallel_translation_numerical_new(self,curve,t,v0,tinterval): 
+        """
+        This method gives the numerical solution to the equation of parallel translation of a vector 
+
+        INPUT: 
+        curve equation = list of functions which determine the curve wrt the local coordinate
+
+        t - curve parameter
+
+        v0 - initial vector
+
+        tinterval = [a,b,N], (a,b) is the domain of the curve, N is the number of subdivision points 
+
+        OUTPUT:
+        The list consisting of lists [t, [v1(t),v2(t)]], where
+
+        t is a subdivision point
+
+        [v1(t),v2(t)] is the list of coordinates of the vector translated parallely along the curve 
+
+        EXAMPLES::
+
+           sage: var('p,q')
+           sage: v = [p,q]
+           sage: assume(cos(q)>0)
+           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
+           sage: var('ss')
+           sage: vv_array = sphere.parallel_translation_numerical([ss,ss],ss,[1,1],[0,pi/4,20])
+           sage: vv_array[0]
+           [0.0, [1.0, 1.0]]
+           sage: vv_array[5]
+           [0.19634954084936207, [0.98060182522638917, 1.0389930474425113]]
+
+        """
+
+        u1 = self.variables[1]
+        u2 = self.variables[2]
+
+        solver = self._create_pt_ode_system(tuple(curve), t)
+            
+        t_interval, n = tinterval[0:2], tinterval[2]
+        solver.y_0 = v0
+        solver.ode_solve(t_span=t_interval, num_points=n)
+
+        return solver.solution
+
+
+
+
+
+    def parallel_translation_numerical(self,curve,t,v0,tinterval): 
+        """
+        This method gives the numerical solution to the equation of parallel translation of a vector 
+
+        INPUT: 
+        curve equation = list of functions which determine the curve wrt the local coordinate
+
+        t - curve parameter
+
+        v0 - initial vector
+
+        tinterval = [a,b,N], (a,b) is the domain of the curve, N is the number of subdivision points 
+
+        OUTPUT:
+        The list consisting of lists [t, [v1(t),v2(t)]], where
+
+        t is a subdivision point
+
+        [v1(t),v2(t)] is the list of coordinates of the vector translated parallely along the curve 
+
+        EXAMPLES::
+
+           sage: var('p,q')
+           sage: v = [p,q]
+           sage: assume(cos(q)>0)
+           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
+           sage: var('ss')
+           sage: vv_array = sphere.parallel_translation_numerical([ss,ss],ss,[1,1],[0,pi/4,20])
+           sage: vv_array[0]
+           [0.0, [1.0, 1.0]]
+           sage: vv_array[5]
+           [0.19634954084936207, [0.98060182522638917, 1.0389930474425113]]
+
+        """
+        u1 = self.variables[1]
+        u2 = self.variables[2]
+        #var('uu1,uu2')
+        
+        #f(u1,u2) = self.equation 
+
+        def f(ux, vx):
+            return [self.equation[0].subs({u1: ux, u2: vx}),
+                    self.equation[1].subs({u1: ux, u2: vx}),
+                    self.equation[2].subs({u1: ux, u2: vx})]
+
+        #from sage.calculus.var import var
+        #t = var('t')
+        tt = t
+
+        C111 = self.connection_coefficients()[(1,1,1)] 
+        C121 = self.connection_coefficients()[(1,2,1)] 
+        C221 = self.connection_coefficients()[(2,2,1)] 
+        C112 = self.connection_coefficients()[(1,1,2)] 
+        C122 = self.connection_coefficients()[(1,2,2)] 
+        C222 = self.connection_coefficients()[(2,2,2)] 
+
+        du1=diff(curve[0],tt)
+        du2=diff(curve[1],tt)
+        uu1=curve[0]
+        uu2=curve[1]
+
+
+        def ode_system(unknown_functions,t1):
+            du1p = du1.subs(tt==t1)
+            du2p = du2.subs(tt==t1)
+            uu1p = uu1.subs(tt==t1)
+            uu2p = uu2.subs(tt==t1)
+
+            c111 = C111.subs({u1: uu1p, u2: uu2p})
+            c121 = C121.subs({u1: uu1p, u2: uu2p})
+            c221 = C221.subs({u1: uu1p, u2: uu2p})
+            c112 = C112.subs({u1: uu1p, u2: uu2p})
+            c122 = C122.subs({u1: uu1p, u2: uu2p})
+            c222 = C222.subs({u1: uu1p, u2: uu2p})
+
+            #print c111, c121, c221, c112, c122, c222
+            
+            v1,v2 = unknown_functions
+            f1 = - c111*du1p*v1 - c121*(du1p*v2 + du2p*v1) - c221*du2p*v2 
+            f2 = - c112*du1p*v1 - c122*(du1p*v2 + du2p*v1) - c222*du2p*v2 
+            dv1 = f1.subs(tt==t1)
+            dv2 = f2.subs(tt==t1)
+
+            #print dv1, dv2
+            
+            return float(dv1), float(dv2) 
+
+        step = float((tinterval[1]-tinterval[0])/tinterval[2])
+        ttt =  float(tinterval[0]) 
+        tarray = [ ttt ]
+        for counter in range(1,tinterval[2]+1):
+           ttt = ttt + step 
+           tarray = tarray + [ ttt ] 
+
+
+        # import ode solving routine 
+        import scipy.integrate
+
+        initial_data = (v0[0],v0[1])
+        solution = scipy.integrate.odeint(ode_system,initial_data,tarray)
+        return [[ tarray[counter],[solution[counter,0], solution[counter,1]]]   for counter in range(0,len(tarray))]
+

new file sage/riemann/vector_functions.py

@@ +1 @@
+### TODO: Some of these functions have Sage equivalents -- need to use these instead.
+
+from sage.modules.free_module_element import vector
+from sage.matrix.constructor import matrix
+
+# vector functions 
+#def scalar_product(a,b):
+#    return(a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
+#def norm(a):
+#    return(sqrt(scalar_product(a,a)))
+def vector_product(a,b):
+  return(vector([a[1]*b[2]-a[2]*b[1],a[2]*b[0]-
+  a[0]*b[2],a[0]*b[1]-a[1]*b[0]]))
+def mixed_product(a,b,c):
+    AA=matrix([[a[0],a[1],a[2]],[b[0],b[1],b[2]],[c[0],c[1],c[2]]])
+    return AA.determinant()
+
+# standard vector functions
+def vector_function_e(u):
+    return vector([cos(u),sin(u),0])
+def vector_function_g(u):
+    return vector([-sin(u),cos(u),0])
+vector_k=vector([0,0,1])
+
+    
+# simplify list of functions
+def simplify_list(a):
+    return([i.simplify_full() for i in a])
+    
+def simplify_vector(a):
+    return(vector([i.simplify_full() for i in a]))
+
+def simplify_matrix(A):
+    return matrix([[ xx.simplify_full() for xx in row ] for row in A])
+
+

setup.py

@@ -985 +985 @@
                      'sage.quadratic_forms',
                      'sage.quadratic_forms.genera',
 
+                     'sage.riemann',
+
                      'sage.rings',
                      'sage.rings.finite_rings',
                      'sage.rings.number_field',