Ticket #10132: trac_10132_differential_geometry_sage_v3.patch

File trac_10132_differential_geometry_sage_v3.patch, 43.5 KB (added by jvkersch, 10 years ago)
  • sage/all.py

    # HG changeset patch
    # User Joris Vankerschaver <joris.vankerschaver@gmail.com>
    # Date 1288935163 25200
    # Node ID d7b6542c48b80af1939dd72192a7dd2badd5694f
    # Parent  120c07be6358d93bcff503363d379c26b8342f2b
    Parametrized Surface 3D
    
    diff -r 120c07be6358 -r d7b6542c48b8 sage/all.py
    a b  
    139139from sage.ext.fast_eval      import fast_float
    140140
    141141from sage.tensor.all     import *
     142from sage.riemann.all    import *
    142143
    143144from copy import copy, deepcopy
    144145
  • new file sage/riemann/__init__.py

    diff -r 120c07be6358 -r d7b6542c48b8 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

    diff -r 120c07be6358 -r d7b6542c48b8 sage/riemann/all.py
    - +  
     1from parametrized_surface3d import ParametrizedSurface3D
  • new file sage/riemann/parametrized_surface3d.py

    diff -r 120c07be6358 -r d7b6542c48b8 sage/riemann/parametrized_surface3d.py
    - +  
     1"""
     2Basic invariants of parametrized surfaces.
     3
     4
     5AUTHORS:
     6        - Mikhail Malakhaltsev (2010-09-25): initial version
     7
     8"""   
     9#*****************************************************************************
     10#       Copyright (C) 2010  Mikhail Malakhaltsev <mikarm@gmail.com>
     11#
     12#  Distributed under the terms of the GNU General Public License (GPL)
     13#                  http://www.gnu.org/licenses/
     14#*****************************************************************************
     15
     16
     17# The following is maybe too specific ...
     18from vector_functions import *
     19
     20# mikarm: This is for the moment, at the final version will be replaced by import
     21#attach('/home/prince/MY_SAGE_MODULES/vector_functions.py')
     22
     23from sage.structure.sage_object import SageObject
     24from sage.calculus.functional import diff
     25from sage.functions.other import sqrt
     26
     27from sage.misc.cachefunc import cached_method
     28
     29
     30
     31# TODO (distant future): update class...
     32class GenericManifold(SageObject):
     33    pass
     34
     35
     36class ParametrizedSurface3D(GenericManifold):
     37    """
     38    This is a class for finding basic invariants of surfaces.
     39
     40    USAGE:
     41    parametrized_surface3d(surface_equation,variables,name_of_surface)
     42    where
     43      surface_equation is a vector or list of three components
     44      variables is a list of two variables
     45      name_of_surface is a string (optional)
     46
     47    Warning: The orientation on the surface is given by the parametrization, that is the positive frame is
     48    $\partial_1 \\vec r$, $\partial_2 \\vec r$.
     49       
     50    This class includes the following methods:
     51      natural_frame
     52      normal_vector
     53      first_fundamental_form_coefficients
     54      first_fundamental_form
     55      first_fundamental_form_inverse_coefficients
     56      first_fundamental_form_inverse_coefficients
     57      area_form_squared
     58      area_form
     59      rotation
     60      orthonormal_frame
     61      lie_bracket
     62      frame_structure_functions
     63      second_order_natural_frame
     64      second_fundamental_form_coefficients
     65      second_fundamental_form
     66      gauss_curvature
     67      mean_curvature
     68      principal curvatures
     69      principal directions
     70      connection_coefficients
     71      geodesics_numerical
     72      parallel_translation_numerical
     73     
     74
     75    EXAMPLES::
     76
     77      sage: var('u,v')
     78      sage: paraboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'paraboloid')
     79      sage: paraboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v])
     80           
     81    """
     82
     83    def __init__(self,equation,variables,*name):
     84
     85        self.equation = equation
     86        self.variables_list = variables
     87        self.variables = {1:self.variables_list[0],2:self.variables_list[1]}
     88        self.name = name
     89
     90        # Callable version of the underlying equation
     91        def eq_callable(u, v):
     92            u1, u2 = self.variables_list
     93            return [fun.subs({u1: u, u2: v}) for fun in self.equation]
     94           
     95        self.eq_callable = eq_callable
     96
     97        # Various index tuples       
     98        self.index_list = [(i,j) for i in (1,2) for j in (1,2)]
     99        self.index_list_3=[(i,j,k) for i in (1,2) for j in (1,2) for k in (1,2)]
     100
     101
     102    @cached_method           
     103    def natural_frame(self,vector_number=0):
     104        """
     105        This functions find the natural frame of the parametrized surface
     106        INPUT:
     107        the argument can be empty, or equal to 1, or to 2
     108
     109        OUTPUT:
     110        if argument is equal to 1 or 2, the output is the corresponding vector of the natural frame
     111        if argument is empty, the output is the dictionary of the vector fields of the natural frame
     112
     113        EXAMPLES:: 
     114
     115            sage:  eparaboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
     116            sage:  eparaboloid.natural_frame()
     117            {1: (1, 0, 2*u), 2: (0, 1, 2*v)}
     118            sage:  eparaboloid.natural_frame(1)
     119            (1, 0, 2*u)
     120            sage:  eparaboloid.natural_frame(2)
     121            (0, 1, 2*v)
     122        """
     123
     124        dr1 = vector([diff(f,self.variables[1]).simplify_full() for f in self.equation])
     125        dr2 = vector([diff(f,self.variables[2]).simplify_full() for f in self.equation])
     126
     127        return {1:dr1, 2:dr2}
     128
     129
     130    @cached_method
     131    def normal_vector(self, normalized=False):
     132        """
     133        This functions find the normal vector of the parametrized surface
     134       
     135        INPUT:
     136        empty argument, of a nonzero real number l
     137       
     138        OUTPUT:
     139        without arguments gives the normal vector which is the vector product of the
     140        natural frame vectors;
     141
     142        with the argument l gives the normal vector of length l (the orientation determined
     143        by the sign of l)
     144       
     145        EXAMPLES::
     146
     147          sage: eparaboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
     148          sage: eparaboloid.normal_vector()
     149          (-2*u,-2*v,1)
     150          sage: eparaboloid.normal_vector(1)
     151          (-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))
     152
     153        """
     154
     155        dr = self.natural_frame()
     156        normal = dr[1].cross_product(dr[2])
     157
     158        if normalized:
     159            normal /= normal.norm()
     160        return simplify_vector(normal)
     161       
     162
     163
     164    # The following can become part of the docstring at some point...
     165
     166    # I've defined a cached (private) method which does not check its arguments and
     167    # computes one component at a time.  Note that caching is just a matter of
     168    # prepending "@cached_method" to the method signature.  For this to work, the
     169    # arguments of the method must be hashable (e.g tuple, string, ...) and the method
     170    # must not change the internal state of the object.  For all of the methods in this
     171    # class, this is the case.
     172    #
     173    # I then added a method that does simple argument checking, converting the arguments to
     174    # tuples, and invokes the private method to do the work.  Note that this method should
     175    # not be cached, as the user might specify the index of the component as a list, which is
     176    # not hashable.
     177    #
     178    # A second method returns a dictionary of all the components, again invoking the hidden
     179    # helper method.
     180 
     181    @cached_method
     182    def _compute_first_fundamental_form_coefficient(self, index):
     183        dr = self.natural_frame()
     184        return (dr[index[0]]*dr[index[1]]).simplify_full()       
     185   
     186    def first_fundamental_form_coefficient(self, index):
     187        index = tuple(sorted(index))
     188        if index not in self.index_list:
     189            raise ValueError, "Index %s out of bounds." % str(index)
     190        return self._compute_first_fundamental_form_coefficient(index)
     191
     192    @cached_method
     193    def first_fundamental_form_coefficients(self):
     194
     195        """
     196        This function gives the coefficients of the first fundamental form
     197       
     198        INPUT:
     199        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
     200
     201        OUTPUT:
     202        with empty argument the output is the dictionary of coefficients of the first fundamental form
     203
     204        with given indices the output is the corresponding coefficient of the first fundamental form,
     205        for example index1 = 1, index2 = 2 gives the coefficient $g_{12}$
     206
     207        EXAMPLES::
     208
     209            sage: var('u,v')
     210            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     211            sage: sphere.first_fundamental_form_coefficients()
     212            {(1, 2): 0, (1, 1): cos(v)^2, (2, 1): 0, (2, 2): 1}
     213
     214            sage: sphere.first_fundamental_form_coefficients(1,1)
     215            cos(v)^2
     216
     217            sage: sphere.first_fundamental_form_coefficients(2,1)
     218            0
     219
     220            sage: sphere.first_fundamental_form_coefficients(3,2)
     221            'The argument is not appropriate. Read doc'
     222        """       
     223        coefficients = {}
     224        for index in self.index_list:
     225            sorted_index = list(sorted(index))
     226            coefficients[index] = \
     227                self._compute_first_fundamental_form_coefficient(index)
     228        return coefficients
     229
     230
     231
     232    def first_fundamental_form(self,vector1,vector2):
     233        """
     234        Finds the value of first fundamental form on two vectors
     235
     236        INPUT:
     237        Two vectors $v=(v^1,v^2)$ and $w=(w^1,w^2)$
     238
     239        OUTPUT:
     240        $g_{11} v^1 w^1 + g_{12}(v^1 w^2 + v^2 w^1) + g_{22} v^2 w^2$
     241
     242        #EXAMPLES::
     243         
     244            sage: var('u,v')
     245            sage: var ('v1,v2,w1,w2')
     246            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     247
     248            sage: sphere.first_fundamental_form(vector([v1,v2]),vector([w1,w2]))
     249            v1*w1*cos(v)^2 + v2*w2
     250
     251            sage: vv = vector([1,2])
     252            sage: sphere.first_fundamental_form(vv,vv)
     253            cos(v)^2 + 4
     254
     255            sage: sphere.first_fundamental_form([1,1],[2,1])
     256            2*cos(v)^2 + 1
     257        """
     258        gg = self.first_fundamental_form_coefficients()
     259        return sum(gg[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list)
     260
     261
     262    @cached_method
     263    def area_form_squared(self):
     264        """
     265        Gives $g_{11}g_{22} - g_{12}^2$, which is the square of the coefficient of the area form
     266       
     267        INPUT:
     268        No arguments
     269
     270        OUTPUT: 
     271        $g_{11}g_{22} - g_{12}^2$
     272
     273        EXAMPLES::
     274
     275            sage: var('u,v')
     276            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     277            sage: sphere.area_form_squared()
     278            cos(v)^2
     279
     280        """
     281        gg = self.first_fundamental_form_coefficients()
     282        return  (gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2).simplify_full()
     283
     284
     285    @cached_method
     286    def area_form(self):
     287        """
     288        Gives $\sqrt{g_{11}g_{22} - g_{12}^2}$, which is the coefficient of the area form
     289       
     290        INPUT:
     291        No arguments
     292
     293        OUTPUT: 
     294        $\sqrt{g_{11}g_{22} - g_{12}^2}$
     295
     296        EXAMPLES::
     297
     298            sage: var('u,v')
     299            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     300            sage: sphere.area_form_squared()
     301            abs(cos(v))
     302
     303        """
     304        return sqrt(self.area_form_squared()).simplify_full()
     305
     306
     307    @cached_method
     308    def first_fundamental_form_inverse_coefficients(self):
     309        """
     310        Gives $g^{ij}$
     311       
     312        INPUT:
     313        No arguments
     314
     315        OUTPUT: 
     316        dictionary {(1,1):$g^{11}$, (1,2):$g^{12}$, (2,1):$g^{21}$, (2,2):$g^{22}$
     317
     318        EXAMPLES::
     319
     320            sage: var('u,v')
     321            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     322            sage: sphere.first_fundamental_form_inverse_coefficients()
     323            {(1, 2): 0, (1, 1): cos(v)^(-2), (2, 1): 0, (2, 2): 1}
     324        """
     325
     326        gg = self.first_fundamental_form_coefficients()
     327        DD = gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2
     328       
     329        gi11 = (gg[(2,2)]/DD).simplify_full()
     330        gi12 = (-gg[(1,2)]/DD).simplify_full()
     331        gi21 = gi12
     332        gi22 = (gg[(1,1)]/DD).simplify_full()
     333
     334        return {(1,1):gi11,(1,2):gi12,(2,1):gi21,(2,2):gi22}
     335
     336    def first_fundamental_form_inverse_coefficient(self, index):
     337        index = tuple(sorted(index))
     338        if index not in self.index_list:
     339            raise ValueError, "Index %s out of bounds." % str(index)
     340        return self._compute_first_fundamental_form_inverse_coefficients()[index]
     341
     342
     343    @cached_method
     344    def rotation(self,theta):
     345        """
     346        Gives the matrix of the operator of rotation on the given angle $\\theta$ with respect to the natural frame
     347       
     348        INPUT:
     349        given angle $\\theta$
     350
     351        OUTPUT: 
     352        matrix of the operator of rotation on $\\theta$ with respect to the natural frame
     353
     354        ALGORITHM:
     355        The operator of rotation on $\pi/2$ is $J^i_j = g^{ik}\omega_{jk}$, where $\omega$ is the area form
     356        The operator of rotation on angle $\\theta$ is $\cos(\\theta) I + sin(\\theta) J$
     357
     358        EXAMPLES::
     359
     360            sage: var('u,v')
     361            sage: assume(cos(v)>0)
     362            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     363            sage: rotation = sphere.rotation(pi/3)
     364            sage: rotation^3
     365            [-1  0]
     366            [ 0 -1]
     367            # it is true because the rotation at $\pi/3$ applied three times gives rotation a $\pi$
     368        """
     369
     370        from sage.functions.trig import sin, cos
     371
     372        gi = self.first_fundamental_form_inverse_coefficients()
     373        w12 = self.area_form()
     374        RR11 = (cos(theta) + sin(theta)*gi[1,2]*w12).simplify_full()
     375        RR12 = (- sin(theta)*gi[1,1]*w12).simplify_full()
     376        RR21 = (sin(theta)*gi[2,2]*w12).simplify_full()
     377        RR22 = (cos(theta) - sin(theta)*gi[2,1]*w12).simplify_full()       
     378        rotation = matrix([[RR11,RR12],[RR21,RR22]])
     379        return rotation
     380
     381
     382    # XXXXXX
     383
     384
     385
     386
     387
     388    @cached_method
     389    def orthonormal_frame_vector(self, vector_number=0, coordinates='ext'):
     390        """
     391        Gives the orthonormal frame field of the surface 
     392       
     393        INPUT:
     394        (vector_number,coordinates), where
     395        vector number is 0 (default), 1 or 2
     396        coordinates is 'ext'(default) or 'int'
     397               
     398
     399        OUTPUT:
     400        If vector_number is equal to 1 or 2, the output is the corresponding vector of an orthonormal frame.
     401        If argument vector_number is empty, the output is the dictionary of the vector fields of an orthonormal frame.
     402       
     403        If coordinates='ext', output is coordinates in $\mathbb{R}^3$
     404        If coordinates='int', output is coordinates with respect to the natural frame
     405
     406        ALGORITHM:
     407        We normalize the first vector $\\vec e_1$ of the natural frame and then get the second frame vector
     408        $\\vec e_2 = [\\vec n, \\vec e_1]$.
     409
     410        EXAMPLES::
     411
     412            sage: var('u,v')
     413            sage: assume(cos(v)>0)
     414            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     415            sage: EE = sphere.orthonormal_frame()
     416            sage: (EE[1]*EE[1]).simplify_full()
     417            1
     418            sage: (EE[1]*EE[2]).simplify_full()
     419            0
     420            sage: simplify_vector(vector_product(EE[1],EE[2])-sphere.normal_vector(1))
     421            [0,0,0]
     422
     423            sage: EE_int = sphere.orthonormal_frame(coordinates='int')
     424            sage: EE_int
     425            {1: (1/cos(v), 0), 2: (0, 1)}
     426            sage: sphere.first_fundamental_form(EE_int[1],EE_int[1])
     427            1
     428            sage: sphere.first_fundamental_form(EE_int[1],EE_int[2])
     429            0
     430            sage: sphere.first_fundamental_form(EE_int[2],EE_int[2])
     431            1
     432           
     433        """
     434
     435        if vector_number not in [1, 2]:
     436            raise ValueError, "Basis vector number out of range."
     437        return self.orthonormal_frame(coordinates)[vector_number]
     438
     439
     440    @cached_method
     441    def orthonormal_frame(self, coordinates='ext'):
     442
     443        from sage.symbolic.constants import pi
     444
     445        if coordinates not in ['ext', 'int']:
     446            raise ValueError, \
     447                  r"Coordinate system must be exterior ('ext') or interior ('int')."
     448
     449        c  = self.first_fundamental_form_coefficient([1,1])
     450        if coordinates == 'ext':
     451            f1 = self.natural_frame()[1]
     452
     453            E1 = simplify_vector(f1/sqrt(c))
     454            E2 = simplify_vector(self.normal_vector(normalized=True).cross_product(E1))
     455        else:
     456            E1 =  vector([(1/sqrt(c)).simplify_full(), 0])
     457            E2 = simplify_vector(self.rotation(pi/2)*E1)
     458        return  {1:E1, 2:E2}
     459       
     460       
     461    def lie_bracket(self,v,w):
     462        """
     463        Gives the Lie bracket of two vector fields which are given by coordinates with respect to the natural frame
     464         
     465       
     466        INPUT:
     467        v and w are vectors
     468               
     469
     470        OUTPUT:
     471        vector [v,w]
     472
     473
     474        EXAMPLES::
     475
     476            sage: var('u,v')
     477            sage: assume(cos(v)>0)
     478            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     479            sage: sphere.lie_bracket([u,v],[-v,u])
     480            (0, 0)
     481
     482            sage: EE_int = sphere.orthonormal_frame(coordinates='int')
     483            sage: sphere.lie_bracket(EE_int[1],EE_int[2])
     484            (sin(v)/cos(v)^2, 0)
     485        """
     486        vv = vector([xx for xx in v])
     487        ww = vector([xx for xx in w])
     488        uuu = [self.variables[1],self.variables[2]]
     489        Dvv = matrix([ [ diff(xxx,uu).simplify_full() for uu in uuu ] for xxx in vv])   
     490        Dww = matrix([ [ diff(xxx,uu).simplify_full() for uu in uuu ] for xxx in ww])   
     491        return simplify_vector(vector(Dvv*ww - Dww*vv))
     492
     493
     494
     495    def frame_structure_functions(self,e1,e2):
     496        """
     497        Gives the structure functions $c^k_{ij}$ for a frame field $e_1$, $e_2$, where
     498        $[e_i,e_j] = c^k_{ij}e_k$
     499         
     500       
     501        INPUT:
     502        Frame vectors $e_1$, $e_2$
     503               
     504
     505        OUTPUT:
     506        The dictionary of $c^k_{ij}$.
     507        Warning: the tuple (i,j,k) corresponds to $c^k_{ij}$
     508
     509
     510        EXAMPLES::
     511
     512            sage: var('u,v')
     513            sage: assume(cos(v)>0)
     514            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     515            sage: sphere.frame_structure_functions([u,v],[-v,u])
     516            {(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}
     517            sage: structfun = sphere.frame_structure_functions(EE_int[1],EE_int[2])
     518            sage: structfun
     519            {(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}
     520            sage: sphere.lie_bracket(EE_int[1],EE_int[2]) - CC[(1,2,1)]*EE_int[1] - CC[(1,2,2)]*EE_int[2]
     521            (0, 0)
     522            """
     523        ee1 = vector([xx for xx in e1])
     524        ee2 = vector([xx for xx in e2])
     525       
     526        lb = simplify_vector(self.lie_bracket(ee1,ee2))
     527        trans_matrix = matrix([[ee1[0],ee2[0]],[ee1[1],ee2[1]]])
     528        zz = simplify_vector(trans_matrix.inverse()*lb)
     529        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}
     530
     531
     532    @cached_method
     533    def _compute_second_order_frame_element(self, index):
     534        variables = [self.variables[i] for i in index]
     535        ddr_element = vector([diff(f, variables).simplify_full() for f in self.equation])
     536       
     537        return ddr_element
     538
     539    @cached_method
     540    def second_order_natural_frame(self):
     541       
     542        """
     543        Gives the second derivatives of the equation $\\vec r = \\vec r(u^1,u^2)$ of parametrized surface
     544         
     545       
     546        INPUT:
     547        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
     548               
     549
     550        OUTPUT:
     551        With empty argument the output is the dictionary of $\partial_{ij}\\vec r(u^1,u^2)$.
     552
     553        With given indices the output is the corresponding second derivative of the surface parametric equation,
     554        For example, index1 = 1, index2 = 2 gives $\partial_{12}\\vec r(u^1,u^2)$.
     555
     556
     557        EXAMPLES::
     558
     559            sage: var('u,v')
     560            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     561            sage: sphere.second_order_natural_frame()
     562            {(1, 2): (sin(u)*sin(v), -sin(v)*cos(u), 0), (1, 1): (-cos(u)*cos(v),
     563            -sin(u)*cos(v), 0), (2, 1): (sin(u)*sin(v), -sin(v)*cos(u), 0), (2, 2):
     564            (-cos(u)*cos(v), -sin(u)*cos(v), -sin(v))}
     565           
     566            sage: sphere.second_order_natural_frame(1,1)
     567            (-cos(u)*cos(v), -sin(u)*cos(v), 0)
     568            sage: sphere.second_order_natural_frame(1,2)
     569            (sin(u)*sin(v), -sin(v)*cos(u), 0)
     570            sage: sphere.second_order_natural_frame(2,2)
     571            (-cos(u)*cos(v), -sin(u)*cos(v), -sin(v)
     572            """
     573
     574        vectors = {}
     575        for index in self.index_list:
     576            sorted_index = tuple(sorted(index))
     577            vectors[index] = \
     578                self._compute_second_order_frame_element(sorted_index)
     579        return vectors
     580
     581    def second_order_natural_frame_element(self, index):
     582        index = tuple(sorted(index))
     583        if index not in self.index_list:
     584            raise ValueError, "Index %s out of bounds." % str(index)
     585        return self._compute_second_order_frame_element(index)
     586
     587
     588    @cached_method
     589    def _compute_second_fundamental_form_coefficient(self, index):
     590        NN = self.normal_vector(normalized=True)
     591        v  = self.second_order_natural_frame_element_new(index)
     592        return (v*NN).simplify_full()
     593
     594       
     595    def second_fundamental_form_coefficient(self, index):
     596        index = tuple(index)
     597        if index not in self.index_list:
     598            raise ValueError, "Index %s out of bounds." % str(index)
     599        return self._compute_second_fundamental_form_coefficient(index)
     600
     601
     602    @cached_method
     603    def second_fundamental_form_coefficients(self):
     604        """
     605        This function gives the coefficients $h_{ij}$ of the second fundamental form.
     606       
     607        INPUT:
     608        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
     609
     610        OUTPUT:
     611        with empty argument the output is the dictionary of coefficients of the second fundamental form
     612
     613        with given indices the output is the corresponding coefficient of the second fundamental form,
     614        for example index1 = 1, index2 = 2 gives the coefficient $h_{12}$
     615
     616        EXAMPLES::
     617
     618            sage: var('u,v')
     619            sage: assume(cos(v)>0)
     620            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     621
     622            sage: sphere.second_fundamental_form_coefficients()
     623            {(1, 2): 0, (1, 1): -cos(v)^2, (2, 1): 0, (2, 2): -1}
     624            sage: sphere.second_fundamental_form_coefficients(1,1)
     625            -cos(v)^2
     626            sage: sphere.second_fundamental_form_coefficients(2,1)
     627            0
     628
     629            sage: sphere.second_fundamental_form_coefficients(3,2)
     630            'The argument is not appropriate. Read doc'
     631        """
     632
     633        coefficients = {}
     634        for index in self.index_list:
     635            coefficients[index] = \
     636                self._compute_second_fundamental_form_coefficient(index)
     637        return coefficients
     638
     639
     640    def second_fundamental_form(self,vector1,vector2):
     641        """
     642        Finds the value of second fundamental form on two vectors
     643
     644        INPUT:
     645        Two vectors $v=(v^1,v^2)$ and $w=(w^1,w^2)$
     646
     647        OUTPUT:
     648        $h_{11} v^1 w^1 + h_{12}(v^1 w^2 + v^2 w^1) + h_{22} v^2 w^2$
     649
     650        EXAMPLES::
     651
     652           sage: var('u,v')
     653           sage: var ('v1,v2,w1,w2')
     654           sage: assume(cos(v) > 0)
     655           sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     656           sage: sphere.second_fundamental_form(vector([v1,v2]),vector([w1,w2]))
     657           -v1*w1*cos(v)^2 - v2*w2
     658           sage: vv = vector([1,2])
     659           sage: sphere.second_fundamental_form(vv,vv)
     660           -cos(v)^2 - 4
     661           sage: sphere.second_fundamental_form([1,1],[2,1])
     662           -2*cos(v)^2 - 1
     663                   
     664        """
     665        hh = self.second_fundamental_form_coefficients()
     666        return sum(hh[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list)
     667
     668
     669    @cached_method
     670    def gauss_curvature(self):
     671        """
     672        Finds the gaussian curvature $K = \\frac{h_{11}h_{22} - h_{12}^2}{g_{11}g_{22} - g_{12}^2}$.
     673
     674        INPUT:
     675        No arguments
     676
     677        OUTPUT:
     678        $K = \\frac{h_{11}h_{22} - h_{12}^2}{g_{11}g_{22} - g_{12}^2}$
     679
     680        EXAMPLES::
     681
     682           sage: var('R')
     683           sage: assume(R>0)
     684           sage: var('u,v')
     685           sage: assume(cos(v)>0)
     686           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     687           sage: sphere.gauss_curvature()
     688           R^(-2)
     689                   
     690        """
     691        hh = self.second_fundamental_form_coefficients()
     692        return ((hh[(1,1)]*hh[(2,2)]-hh[(1,2)]**2)/self.area_form_squared()).simplify_full()
     693
     694
     695    @cached_method
     696    def mean_curvature(self):
     697        """
     698        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}$.
     699
     700        INPUT:
     701        No arguments
     702
     703        OUTPUT:
     704        $H = \\frac{1}{2}\\frac{g_{22}h_{11} - 2g_{12}h_{12} + g_{11}h_{22}}{g_{11}g_{22} - g_{12}^2}$.
     705
     706        EXAMPLES::
     707           
     708           sage: var('R')
     709           sage: assume(R>0)
     710           sage: var('u,v')
     711           sage: assume(cos(v)>0)
     712           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     713           sage: sphere.mean_curvature()
     714           -1/R
     715
     716        """
     717        gg = self.first_fundamental_form_coefficients()
     718        hh = self.second_fundamental_form_coefficients()
     719        denom = 2*self.area_form_squared()
     720        enum =(gg[(2,2)]*hh[(1,1)]-2*gg[(1,2)]*hh[(1,2)]+gg[(1,1)]*hh[(2,2)]).simplify_full()
     721        return (enum/denom).simplify_full()
     722
     723
     724    @cached_method
     725    def principal_curvatures(self):
     726        """
     727        Finds the principal curvatures of the surface
     728
     729        INPUT:
     730        No arguments
     731
     732        OUTPUT:
     733        The dictionary of principal curvatures         
     734
     735        EXAMPLES::
     736                     
     737           sage: var('R')
     738           sage: assume(R>0)
     739           sage: var('u,v')
     740           sage: assume(cos(v)>0)
     741           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     742           sage: sphere.principal_curvatures()
     743           {1: -1/R, 2: -1/R}
     744
     745           sage: var('u,v')
     746           sage: var('R,r')
     747           sage: assume(R>r,r>0)
     748           sage: torus = parametrized_surface3d([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
     749           sage: torus.principal_curvatures()
     750           {1: -cos(v)/(r*cos(v) + R), 2: -1/r}
     751
     752        """
     753
     754        from sage.symbolic.assumptions import assume
     755        from sage.symbolic.relation import solve
     756        from sage.calculus.var import var
     757
     758        KK = self.gauss_curvature()
     759        HH = self.mean_curvature()
     760
     761        # jvkersch: when this assumption is uncommented, Sage raises an error stating that the assumption
     762        # is redundant... Can we safely omit this, based on some geometric reasoning?
     763       
     764        # assume(HH**2-KK>=0)
     765        # mikarm: This is a problem, I had a lot of trobles here. Of course, this inequality always hold true.
     766        # I insert this assumption because sage sometimes, in simplification, uses the complex numbers, though the roots
     767        # are, for sure, real. 
     768        # This, in turn, causes problems when we substitute coordinates into the expression of principal curvatures.
     769        # Unfortunately, I did not manage to tell Sage that they are real (to declare the variables as real).
     770        # Moreover, in a neighborhood of an umbilic point we even cannot "smoothly" order the set of principal curvatures.
     771        # So, in general, at present this method is far from the final form.
     772         
     773
     774 
     775        x = var('x')
     776        sol = solve(x**2 -2*HH*x + KK == 0,x)
     777       
     778        #k1=var('k1')
     779        #k2=var('k2')
     780
     781        # 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.
     782
     783        k1 = (x.substitute(sol[0])).simplify_full()
     784        if len(sol) == 1:
     785            k2 = k1
     786        else:
     787            k2 = (x.substitute(sol[1])).simplify_full()
     788       
     789        return {1:k1,2: k2}
     790
     791    @cached_method
     792    def principal_directions(self):
     793        """
     794        Finds the principal curvatures and principal directions of the surface
     795
     796        INPUT:
     797        No arguments
     798
     799        OUTPUT:
     800        The dictionary of lists [a principal curvature, the corresponding principal direction]         
     801
     802        If principal curvatures coincide, gives the warning that the surface is a sphere.
     803
     804        EXAMPLES::
     805
     806           sage: var('u,v')
     807           sage: var('R,r')
     808           sage: assume(R>r,r>0)
     809           sage: torus = parametrized_surface3d([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
     810           sage: pdd = torus.principal_directions()
     811           sage: pdd[1]
     812           [-cos(v)/(r*cos(v) + R), (1, 0)]
     813           sage: pdd[2]
     814           [-1/r, (0, -(R*r*cos(v) + R^2)/r)]
     815
     816           sage: var('RR')
     817           sage: assume(RR>0)
     818           sage: var('u,v')
     819           sage: assume(cos(v)>0)
     820           sage: sphere = parametrized_surface3d([RR*cos(u)*cos(v),RR*sin(u)*cos(v),RR*sin(v)],[u,v],'sphere')
     821           sage: sphere.principal_directions()
     822           'This is a sphere, so any direction is principal'
     823
     824           sage: var('aa')
     825           sage: assume(aa>0)
     826           sage: catenoid = parametrized_surface3d([aa*cosh(v)*cos(u),aa*cosh(v)*sin(u),v],[u,v],'catenoid')
     827           sage: pd = catenoid.principal_directions()
     828           sage: pd[1][1]
     829           (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))
     830           sage: pd[2][1]
     831           (1, 0)
     832           sage: pd[1][1]*pd[2][1]
     833           0
     834        """
     835        gg = self.first_fundamental_form_coefficients()
     836        hh = self.second_fundamental_form_coefficients()
     837        kk = self.principal_curvatures()
     838        if kk[1]==kk[2]:
     839            return "This is a sphere, so any direction is principal"
     840        pd1 = simplify_vector([hh[(1,2)]-kk[1]*gg[(1,2)],-hh[(1,1)]+kk[1]*gg[(1,1)]])
     841        if pd1==vector([0,0]):
     842           pd1 = vector([1,0])
     843        pd2 = simplify_vector([hh[(1,2)]-kk[2]*gg[(1,2)],-hh[(1,1)]+kk[2]*gg[(1,1)]])
     844        if pd2==vector([0,0]):
     845           pd2 = vector([1,0])
     846        return {1:[kk[1],pd1],2:[kk[2],pd2]}
     847
     848
     849    @cached_method
     850    def connection_coefficients(self):
     851        """
     852        Finds the connection coefficients of the surface
     853
     854        INPUT:
     855        No arguments
     856
     857        OUTPUT:
     858        The dictionary of connection coefficients.
     859
     860        Warning: the triple $(i,j,k)$ corresponds to $\Gamma^k_{ij}$. 
     861
     862        EXAMPLES::
     863                                 
     864           sage: var('r')
     865           sage: assume(r > 0)
     866           sage: var('u,v')
     867           sage: assume(cos(v)>0)
     868           sage: sphere = parametrized_surface3d([r*cos(u)*cos(v),r*sin(u)*cos(v),r*sin(v)],[u,v],'sphere')
     869           sage: sphere.connection_coefficients()
     870           {(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}
     871           
     872        """
     873        x = self.variables
     874        gg = self.first_fundamental_form_coefficients()
     875        gi = self.first_fundamental_form_inverse_coefficients()
     876
     877        dg = {}
     878        for kkk in self.index_list_3:
     879            dg[kkk]=gg[(kkk[1],kkk[2])].differentiate(x[kkk[0]]).simplify_full()
     880        structfun={}
     881
     882        for kkk in self.index_list_3:
     883            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()
     884        return structfun
     885
     886
     887    # jvkersch: this private method creates an ode_solver object, which can be used
     888    # to integrate the geodesic equations numerically.
     889   
     890    @cached_method
     891    def _create_geodesic_ode_system(self):
     892        from sage.ext.fast_eval import fast_float
     893        from sage.calculus.var import var
     894        from sage.gsl.ode import ode_solver
     895
     896        u1 = self.variables[1]
     897        u2 = self.variables[2]
     898        v1, v2 = var('v1, v2')
     899
     900        C = self.connection_coefficients()
     901       
     902        dv1 = - C[(1,1,1)]*v1**2 - 2*C[(1,2,1)]*v1*v2 - C[(2,2,1)]*v2**2
     903        dv2 = - C[(1,1,2)]*v1**2 - 2*C[(1,2,2)]*v1*v2 - C[(2,2,2)]*v2**2
     904        fun1 = fast_float(dv1, str(u1), str(u2), str(v1), str(v2))
     905        fun2 = fast_float(dv2, str(u1), str(u2), str(v1), str(v2))
     906
     907        geodesic_ode = ode_solver()
     908        geodesic_ode.function = \
     909                              lambda t, (u1, u2, v1, v2) : \
     910                              [v1, v2, fun1(u1, u2, v1, v2), fun2(u1, u2, v1, v2)]
     911        return geodesic_ode
     912
     913
     914    # jvkersch: integrate the geodesic equations numerically
     915    # mikarm: Very good. I cheked it restarting each time the worksheet, it works much faster.
     916       
     917    def geodesics_numerical(self, p0, v0, tinterval):
     918        """
     919        This method gives the numerical solution for the geodesic equations 
     920
     921        INPUT:
     922        p0 is the list of the coordinates of the initial point
     923       
     924        v0 is the  list of the coordinates of the initial vector
     925
     926        tinterval is the list [a,b,M], where (a,b) is the domain of the geodesic, M is the number of division points
     927
     928        OUTPUT:
     929        The list consisting of lists [t, [u1(t),u2(t)], [v1(t),v2(t)], [x1(t),x2(t),x3(t)]], where
     930
     931        t is a subdivision point
     932
     933        [u1(t),u2(t)] is the list of coordinates of the geodesic point
     934         
     935        [v1(t),v2(t)] is the list of coordinates of the vector tanget to the geodesic 
     936
     937        [x1(t),x2(t),x3(t)] is the list of coordinates of the geodesic point in the three-dimensional space
     938
     939
     940        EXAMPLES::
     941
     942           sage: var('p,q')
     943           sage: v = [p,q]
     944           sage: assume(cos(q)>0)
     945           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
     946           sage: gg_array = sphere.geodesics_numerical([0,0],[1,1],[0,2*pi,5])
     947           sage: gg_array[0]
     948           [0.0, [0.0, 0.0], [1.0, 1.0], (1, 0, 0)]
     949           sage: gg_array[1]
     950           [1.2566370614359172, [0.76440104189216407, 1.8586223516062499], [-0.2838683571264714, 1.9194187087799912], (-0.204895333443, 0.692104654602, 0.692104796553)]
     951                                           
     952         
     953        """
     954
     955        u1 = self.variables[1]
     956        u2 = self.variables[2]
     957
     958        solver = self._create_geodesic_ode_system()
     959           
     960        t_interval, n = tinterval[0:2], tinterval[2]
     961        solver.y_0 = [p0[0], p0[1], v0[0], v0[1]]
     962        solver.ode_solve(t_span=t_interval, num_points=n)
     963
     964        parsed_solution = \
     965          [[vec[0], vec[1][0:2], vec[1][2:], self.eq_callable(vec[1][0], vec[1][1])] \
     966           for vec in solver.solution]
     967
     968        return parsed_solution
     969
     970
     971    @cached_method
     972    def _create_pt_ode_system(self, curve, t):
     973        from sage.ext.fast_eval import fast_float
     974        from sage.calculus.var import var
     975        from sage.gsl.ode import ode_solver
     976
     977        u1 = self.variables[1]
     978        u2 = self.variables[2]
     979        v1, v2 = var('v1, v2')
     980
     981        du1 = diff(curve[0], t)
     982        du2 = diff(curve[1], t)
     983
     984        C = self.connection_coefficients()
     985        for coef in C:
     986            C[coef] = C[coef].subs({u1: curve[0], u2: curve[1]})
     987       
     988        dv1 = - C[(1,1,1)]*v1*du1 - C[(1,2,1)]*(du1*v2 + du2*v1) - C[(2,2,1)]*du2*v2
     989        dv2 = - C[(1,1,2)]*v1*du1 - C[(1,2,2)]*(du1*v2 + du2*v1) - C[(2,2,2)]*du2*v2
     990        fun1 = fast_float(dv1, str(t), str(v1), str(v2))
     991        fun2 = fast_float(dv2, str(t), str(v1), str(v2))
     992
     993        pt_ode = ode_solver()
     994        pt_ode.function = lambda t, (v1, v2): [fun1(t, v1, v2), fun2(t, v1, v2)]
     995        return pt_ode
     996
     997
     998    # mikarm: We should rewrite it like you did for geodesics.
     999    # jvkersch: OK, see this and the above
     1000    def parallel_translation_numerical_new(self,curve,t,v0,tinterval):
     1001        """
     1002        This method gives the numerical solution to the equation of parallel translation of a vector
     1003
     1004        INPUT:
     1005        curve equation = list of functions which determine the curve wrt the local coordinate
     1006
     1007        t - curve parameter
     1008
     1009        v0 - initial vector
     1010
     1011        tinterval = [a,b,N], (a,b) is the domain of the curve, N is the number of subdivision points
     1012
     1013        OUTPUT:
     1014        The list consisting of lists [t, [v1(t),v2(t)]], where
     1015
     1016        t is a subdivision point
     1017
     1018        [v1(t),v2(t)] is the list of coordinates of the vector translated parallely along the curve
     1019
     1020        EXAMPLES::
     1021
     1022           sage: var('p,q')
     1023           sage: v = [p,q]
     1024           sage: assume(cos(q)>0)
     1025           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
     1026           sage: var('ss')
     1027           sage: vv_array = sphere.parallel_translation_numerical([ss,ss],ss,[1,1],[0,pi/4,20])
     1028           sage: vv_array[0]
     1029           [0.0, [1.0, 1.0]]
     1030           sage: vv_array[5]
     1031           [0.19634954084936207, [0.98060182522638917, 1.0389930474425113]]
     1032
     1033        """
     1034
     1035        u1 = self.variables[1]
     1036        u2 = self.variables[2]
     1037
     1038        solver = self._create_pt_ode_system(tuple(curve), t)
     1039           
     1040        t_interval, n = tinterval[0:2], tinterval[2]
     1041        solver.y_0 = v0
     1042        solver.ode_solve(t_span=t_interval, num_points=n)
     1043
     1044        return solver.solution
     1045
     1046
     1047
     1048
     1049
     1050    def parallel_translation_numerical(self,curve,t,v0,tinterval):
     1051        """
     1052        This method gives the numerical solution to the equation of parallel translation of a vector
     1053
     1054        INPUT:
     1055        curve equation = list of functions which determine the curve wrt the local coordinate
     1056
     1057        t - curve parameter
     1058
     1059        v0 - initial vector
     1060
     1061        tinterval = [a,b,N], (a,b) is the domain of the curve, N is the number of subdivision points
     1062
     1063        OUTPUT:
     1064        The list consisting of lists [t, [v1(t),v2(t)]], where
     1065
     1066        t is a subdivision point
     1067
     1068        [v1(t),v2(t)] is the list of coordinates of the vector translated parallely along the curve
     1069
     1070        EXAMPLES::
     1071
     1072           sage: var('p,q')
     1073           sage: v = [p,q]
     1074           sage: assume(cos(q)>0)
     1075           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
     1076           sage: var('ss')
     1077           sage: vv_array = sphere.parallel_translation_numerical([ss,ss],ss,[1,1],[0,pi/4,20])
     1078           sage: vv_array[0]
     1079           [0.0, [1.0, 1.0]]
     1080           sage: vv_array[5]
     1081           [0.19634954084936207, [0.98060182522638917, 1.0389930474425113]]
     1082
     1083        """
     1084        u1 = self.variables[1]
     1085        u2 = self.variables[2]
     1086        #var('uu1,uu2')
     1087       
     1088        #f(u1,u2) = self.equation
     1089
     1090        def f(ux, vx):
     1091            return [self.equation[0].subs({u1: ux, u2: vx}),
     1092                    self.equation[1].subs({u1: ux, u2: vx}),
     1093                    self.equation[2].subs({u1: ux, u2: vx})]
     1094
     1095        #from sage.calculus.var import var
     1096        #t = var('t')
     1097        tt = t
     1098
     1099        C111 = self.connection_coefficients()[(1,1,1)]
     1100        C121 = self.connection_coefficients()[(1,2,1)]
     1101        C221 = self.connection_coefficients()[(2,2,1)]
     1102        C112 = self.connection_coefficients()[(1,1,2)]
     1103        C122 = self.connection_coefficients()[(1,2,2)]
     1104        C222 = self.connection_coefficients()[(2,2,2)]
     1105
     1106        du1=diff(curve[0],tt)
     1107        du2=diff(curve[1],tt)
     1108        uu1=curve[0]
     1109        uu2=curve[1]
     1110
     1111
     1112        def ode_system(unknown_functions,t1):
     1113            du1p = du1.subs(tt==t1)
     1114            du2p = du2.subs(tt==t1)
     1115            uu1p = uu1.subs(tt==t1)
     1116            uu2p = uu2.subs(tt==t1)
     1117
     1118            c111 = C111.subs({u1: uu1p, u2: uu2p})
     1119            c121 = C121.subs({u1: uu1p, u2: uu2p})
     1120            c221 = C221.subs({u1: uu1p, u2: uu2p})
     1121            c112 = C112.subs({u1: uu1p, u2: uu2p})
     1122            c122 = C122.subs({u1: uu1p, u2: uu2p})
     1123            c222 = C222.subs({u1: uu1p, u2: uu2p})
     1124
     1125            #print c111, c121, c221, c112, c122, c222
     1126           
     1127            v1,v2 = unknown_functions
     1128            f1 = - c111*du1p*v1 - c121*(du1p*v2 + du2p*v1) - c221*du2p*v2
     1129            f2 = - c112*du1p*v1 - c122*(du1p*v2 + du2p*v1) - c222*du2p*v2
     1130            dv1 = f1.subs(tt==t1)
     1131            dv2 = f2.subs(tt==t1)
     1132
     1133            #print dv1, dv2
     1134           
     1135            return float(dv1), float(dv2)
     1136
     1137        step = float((tinterval[1]-tinterval[0])/tinterval[2])
     1138        ttt =  float(tinterval[0])
     1139        tarray = [ ttt ]
     1140        for counter in range(1,tinterval[2]+1):
     1141           ttt = ttt + step
     1142           tarray = tarray + [ ttt ]
     1143
     1144
     1145        # import ode solving routine
     1146        import scipy.integrate
     1147
     1148        initial_data = (v0[0],v0[1])
     1149        solution = scipy.integrate.odeint(ode_system,initial_data,tarray)
     1150        return [[ tarray[counter],[solution[counter,0], solution[counter,1]]]   for counter in range(0,len(tarray))]
     1151
  • new file sage/riemann/vector_functions.py

    diff -r 120c07be6358 -r d7b6542c48b8 sage/riemann/vector_functions.py
    - +  
     1### TODO: Some of these functions have Sage equivalents -- need to use these instead.
     2
     3from sage.modules.free_module_element import vector
     4from sage.matrix.constructor import matrix
     5
     6# vector functions
     7#def scalar_product(a,b):
     8#    return(a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
     9#def norm(a):
     10#    return(sqrt(scalar_product(a,a)))
     11def vector_product(a,b):
     12  return(vector([a[1]*b[2]-a[2]*b[1],a[2]*b[0]-
     13  a[0]*b[2],a[0]*b[1]-a[1]*b[0]]))
     14def mixed_product(a,b,c):
     15    AA=matrix([[a[0],a[1],a[2]],[b[0],b[1],b[2]],[c[0],c[1],c[2]]])
     16    return AA.determinant()
     17
     18# standard vector functions
     19def vector_function_e(u):
     20    return vector([cos(u),sin(u),0])
     21def vector_function_g(u):
     22    return vector([-sin(u),cos(u),0])
     23vector_k=vector([0,0,1])
     24
     25   
     26# simplify list of functions
     27def simplify_list(a):
     28    return([i.simplify_full() for i in a])
     29   
     30def simplify_vector(a):
     31    return(vector([i.simplify_full() for i in a]))
     32
     33def simplify_matrix(A):
     34    return matrix([[ xx.simplify_full() for xx in row ] for row in A])
     35
     36
  • setup.py

    diff -r 120c07be6358 -r d7b6542c48b8 setup.py
    a b  
    985985                     'sage.quadratic_forms',
    986986                     'sage.quadratic_forms.genera',
    987987
     988                     'sage.riemann',
     989
    988990                     'sage.rings',
    989991                     'sage.rings.finite_rings',
    990992                     'sage.rings.number_field',