Ticket #10132: trac_10132_diff_geom_sage_doctests.patch

File trac_10132_diff_geom_sage_doctests.patch, 55.9 KB (added by jvkersch, 10 years ago)

Latest version of the patch, with 100% coverage and doctests in place.

  • doc/en/reference/index.rst

    # HG changeset patch
    # User Joris Vankerschaver <joris.vankerschaver@gmail.com>
    # Date 1290619376 28800
    # Node ID b834b692db4c5e7fd57c369e54f8763224d2c58a
    # Parent  120c07be6358d93bcff503363d379c26b8342f2b
    #10132: Differential Geometry via Sage
    
    diff -r 120c07be6358 -r b834b692db4c doc/en/reference/index.rst
    a b  
    9292   modabvar
    9393   modmisc
    9494   tensor
     95   riemann
    9596
    9697   history_and_license
    9798
  • new file doc/en/reference/riemann.rst

    diff -r 120c07be6358 -r b834b692db4c doc/en/reference/riemann.rst
    - +  
     1Differential Geometry of Curves and Surfaces
     2============================================
     3
     4.. toctree::
     5   :maxdepth: 2
     6
     7   sage/riemann/parametrized_surface3d
     8
  • sage/all.py

    diff -r 120c07be6358 -r b834b692db4c 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 b834b692db4c 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 b834b692db4c sage/riemann/all.py
    - +  
     1from parametrized_surface3d import ParametrizedSurface3D
  • new file sage/riemann/parametrized_surface3d.py

    diff -r 120c07be6358 -r b834b692db4c sage/riemann/parametrized_surface3d.py
    - +  
     1"""
     2Differential Geometry of Parametrized Surfaces
     3
     4jvkersch: we need to have some more examples of how to use this class here.
     5
     6AUTHORS:
     7        - Mikhail Malakhaltsev (2010-09-25): initial version
     8
     9"""   
     10#*****************************************************************************
     11#       Copyright (C) 2010  Mikhail Malakhaltsev <mikarm@gmail.com>
     12#
     13#  Distributed under the terms of the GNU General Public License (GPL)
     14#                  http://www.gnu.org/licenses/
     15#*****************************************************************************
     16
     17
     18# The following is maybe too specific ...
     19from vector_functions import *
     20
     21# mikarm: This is for the moment, at the final version will be replaced by import
     22#attach('/home/prince/MY_SAGE_MODULES/vector_functions.py')
     23
     24from sage.structure.sage_object import SageObject
     25from sage.calculus.functional import diff
     26from sage.functions.other import sqrt
     27
     28from sage.misc.cachefunc import cached_method
     29
     30
     31
     32class GenericManifold(SageObject):
     33    """
     34    Placeholder class representing a smooth manifold.  Currently this class
     35    serves as a base class for ``ParametrizedSurface3D``, but does not have any
     36    other functionality.
     37   
     38    TODO:
     39   
     40     - Make this class aware of ``CoordinatePatch`` functionality.
     41   
     42    """
     43    pass
     44
     45
     46class ParametrizedSurface3D(GenericManifold):
     47    r"""
     48    This is a class for finding basic invariants of surfaces.
     49
     50    The orientation of the surface is determined by the parametrization, that is, the
     51    natural frame with positive orientation is given by $\partial_1 \vec r$, $\partial_2 \vec r$.         
     52
     53    jvkersch: This description, as well as the doctests, should be expanded.
     54
     55    INPUT:
     56
     57     - ``surface_equation`` - list/vector specifying a parametric representation of the surface
     58
     59     - ``variables`` - coordinates on the surface
     60
     61     - ``name_of_surface`` - string with the name of the surface (optional).
     62       
     63    EXAMPLES::
     64
     65      sage: u, v = var('u,v')
     66      sage: paraboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'paraboloid'); paraboloid
     67      Parametrized surface ('paraboloid') with equation [u, v, u^2 + v^2]
     68
     69    """
     70
     71    def __init__(self,equation,variables,*name):
     72        """
     73        See ``ParametrizedSurface3D`` for full documentation.
     74
     75        EXAMPLES::
     76
     77            sage: u, v = var('u,v')
     78            sage: eq = [3*u + 3*u*v^2 - u^3, 3*v + 3*u^2*v - v^3, 3*(u^2-v^2)]
     79            sage: enneper = ParametrizedSurface3D(eq,[u,v],'enneper_surface')
     80           
     81        """
     82        self.equation = equation
     83        self.variables_list = variables
     84        self.variables = {1:self.variables_list[0],2:self.variables_list[1]}
     85        self.name = name
     86
     87        # Callable version of the underlying equation
     88        def eq_callable(u, v):
     89            u1, u2 = self.variables_list
     90            return [fun.subs({u1: u, u2: v}) for fun in self.equation]
     91           
     92        self.eq_callable = eq_callable
     93
     94        # Various index tuples       
     95        self.index_list   = [(i,j) for i in (1,2) for j in (1,2)]
     96        self.index_list_3 = [(i,j,k) for i in (1,2) for j in (1,2) for k in (1,2)]
     97
     98
     99    def _latex_(self):
     100        r"""
     101        Returns the LaTeX representation of this parametrized surface.
     102
     103        EXAMPLES::
     104
     105            sage: u, v = var('u, v')
     106            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     107            sage: latex(sphere)
     108            \left[\cos\left(u\right) \cos\left(v\right), \sin\left(u\right) \cos\left(v\right), \sin\left(v\right)\right]
     109            sage: sphere._latex_()
     110            \left[\cos\left(u\right) \cos\left(v\right), \sin\left(u\right) \cos\left(v\right), \sin\left(v\right)\right]
     111
     112        """       
     113        from sage.misc.latex import latex
     114        return latex(self.equation)
     115
     116
     117    def _repr_(self):
     118        r"""
     119        Returns the string representation of this parametrized surface.
     120
     121        EXAMPLES::
     122
     123            sage: u, v = var('u, v')
     124            sage: eq = [3*u + 3*u*v^2 - u^3, 3*v + 3*u^2*v - v^3, 3*(u^2-v^2)]
     125            sage: enneper = ParametrizedSurface3D(eq,[u,v],'enneper_surface')
     126            sage: print enneper
     127            Parametrized surface ('enneper_surface') with equation [-u^3 + 3*u*v^2 + 3*u, 3*u^2*v - v^3 + 3*v, 3*u^2 - 3*v^2]
     128            sage: enneper._repr_()
     129            "Parametrized surface ('enneper_surface') with equation [-u^3 + 3*u*v^2 + 3*u, 3*u^2*v - v^3 + 3*v, 3*u^2 - 3*v^2]"
     130
     131        """
     132        name = 'Parametrized surface'
     133        if self.name is not None:
     134            name += " ('%s')" % self.name
     135        return '%(designation)s with equation %(eq)s' % \
     136               {'designation': name, 'eq': str(self.equation)}
     137
     138
     139    def plot3d(self, urange, vrange):
     140        """
     141        Enable easy plotting directly from the surface class.
     142
     143        INPUT:
     144
     145         - ``urange`` - 3-tuple ``(u, u_min, u_max)``
     146         - ``vrange`` - 3-tuple ``(v, v_min, v_max)``
     147         
     148        EXAMPLES::
     149       
     150            sage: u, v = var('u, v')
     151            sage: eq = [3*u + 3*u*v^2 - u^3, 3*v + 3*u^2*v - v^3, 3*(u^2-v^2)]
     152            sage: enneper = ParametrizedSurface3D(eq,[u,v],'enneper_surface')
     153            sage: enneper.plot3d((u, -5, 5), (v, -5, 5))
     154        """
     155        from sage.plot.plot3d.parametric_plot3d import parametric_plot3d
     156        P = parametric_plot3d(self.equation, urange, vrange)
     157        return P
     158
     159
     160    @cached_method           
     161    def natural_frame(self):
     162        """
     163        Returns the natural tangent frame on the parametrized surface.  The vectors of
     164        this frame are tangent to the coordinate lines on the surface.
     165       
     166        OUTPUT:
     167
     168        - The natural frame as a dictionary.
     169       
     170        EXAMPLES:: 
     171
     172            sage: u, v = var('u, v')
     173            sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
     174            sage: eparaboloid.natural_frame()
     175            {1: (1, 0, 2*u), 2: (0, 1, 2*v)}
     176        """
     177
     178        dr1 = vector([diff(f,self.variables[1]).simplify_full() for f in self.equation])
     179        dr2 = vector([diff(f,self.variables[2]).simplify_full() for f in self.equation])
     180
     181        return {1:dr1, 2:dr2}
     182
     183
     184    @cached_method
     185    def normal_vector(self, normalized=False):
     186        """
     187        Returns the normal vector field of the parametrized surface.
     188       
     189        INPUT:
     190       
     191          - ``normalized`` - default ``False`` - specifies whether the normal vector should be normalized.
     192                 
     193        OUTPUT:
     194
     195         - Normal vector field.
     196       
     197        EXAMPLES::
     198
     199            sage: u, v = var('u, v')
     200            sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
     201            sage: eparaboloid.normal_vector(normalized=False)
     202            (-2*u, -2*v, 1)
     203            sage: eparaboloid.normal_vector(normalized=True)
     204            (-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))
     205
     206        """
     207
     208        dr = self.natural_frame()
     209        normal = dr[1].cross_product(dr[2])
     210
     211        if normalized:
     212            normal /= normal.norm()
     213        return simplify_vector(normal)
     214       
     215 
     216    @cached_method
     217    def _compute_first_fundamental_form_coefficient(self, index):
     218        """
     219        Helper function to compute coefficients of the first fundamental form.
     220
     221        Do not call this method directly; instead use ``first_fundamental_form_coefficient``.
     222        This method expects its argument to be a list, and hence can be cached.
     223
     224        EXAMPLES::
     225
     226            sage: u, v = var('u, v')
     227            sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v])
     228            sage: eparaboloid._compute_first_fundamental_form_coefficient((1,2))
     229            4*u*v
     230
     231        """
     232        dr = self.natural_frame()
     233        return (dr[index[0]]*dr[index[1]]).simplify_full()       
     234
     235   
     236    def first_fundamental_form_coefficient(self, index):
     237        r"""
     238        Compute a single component $g_{ij}$ of the first fundamental form.  If the
     239        parametric representation of the surface is given by the vector function
     240        $\vec r(u^i)$, where $u^i$, $i = 1, 2$ are curvilinear coordinates, then
     241
     242         .. math::
     243
     244           g_{ij} = \frac{\partial \vec r}{\partial u^i} \cdot
     245             \frac{\partial \vec r}{\partial u^j}.
     246
     247        INPUT:
     248
     249         - ``index`` - length-2 index ``(i, j)`` of the component
     250
     251        OUTPUT:
     252
     253         - Component ``g_ij`` of the first fundamental form
     254
     255        EXAMPLES::
     256
     257            sage: u, v = var('u, v')
     258            sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v])
     259            sage: eparaboloid.first_fundamental_form_coefficient((1,2))
     260            4*u*v
     261
     262        When the index is invalid, an error is raised::
     263
     264            sage: u, v = var('u, v')
     265            sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v])
     266            sage: eparaboloid.first_fundamental_form_coefficient((1,5))
     267            Traceback (most recent call last):
     268            ...
     269            ValueError: Index (1, 5) out of bounds.
     270           
     271        """
     272        index = tuple(sorted(index))
     273        if index not in self.index_list:
     274            raise ValueError, "Index %s out of bounds." % str(index)
     275        return self._compute_first_fundamental_form_coefficient(index)
     276
     277
     278    @cached_method
     279    def first_fundamental_form_coefficients(self):
     280        """
     281        Returns the coefficients of the first fundamental form as a dictionary. The keys
     282        are tuples $(i, j)$, where $i$ and $j$ range over $1, 2$, while the values are
     283        the corresponding coefficients $g_{ij}$.
     284       
     285        OUTPUT:
     286
     287         - Dictionary of first fundamental form coefficients.
     288       
     289        EXAMPLES::
     290
     291            sage: u, v = var('u,v')
     292            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     293            sage: sphere.first_fundamental_form_coefficients()
     294            {(1, 2): 0, (1, 1): cos(v)^2, (2, 1): 0, (2, 2): 1}
     295
     296        """       
     297        coefficients = {}
     298        for index in self.index_list:
     299            sorted_index = list(sorted(index))
     300            coefficients[index] = \
     301                self._compute_first_fundamental_form_coefficient(index)
     302        return coefficients
     303
     304
     305    def first_fundamental_form(self, vector1, vector2):
     306        r"""
     307        Evaluate the first fundamental form on two vectors expressed with
     308        respect to the natural coordinate frame on the surface. In other words,
     309        if the vectors are $v = (v^1, v^2)$ and $w = (w^1, w^2)$, calculate
     310        $g_{11} v^1 w^1 + g_{12}(v^1 w^2 + v^2 w^1) + g_{22} v^2 w^2$, with
     311        $g_{ij}$ the coefficients of the first fundamental form.
     312       
     313        INPUT:
     314       
     315         - ``vector1``, ``vector2`` - vectors on the surface.
     316
     317        OUTPUT:
     318
     319         - First fundamental form evaluated on the input vectors.
     320
     321        EXAMPLES::
     322         
     323            sage: u, v = var('u, v')
     324            sage: v1, v2, w1, w2 = var('v1, v2, w1, w2')
     325            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     326            sage: sphere.first_fundamental_form(vector([v1,v2]),vector([w1,w2]))
     327            v1*w1*cos(v)^2 + v2*w2
     328
     329            sage: vv = vector([1,2])
     330            sage: sphere.first_fundamental_form(vv,vv)
     331            cos(v)^2 + 4
     332
     333            sage: sphere.first_fundamental_form([1,1],[2,1])
     334            2*cos(v)^2 + 1
     335        """
     336        gg = self.first_fundamental_form_coefficients()
     337        return sum(gg[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list)
     338
     339
     340    @cached_method
     341    def area_form_squared(self):
     342        """
     343        Returns the square of the coefficient of the area form on the surface.
     344        In terms of the coefficients $g_{ij}$ (where $i, j = 1, 2$) of the first fundamental
     345        form, this invariant is given by
     346
     347        .. math::
     348
     349           A^2 = g_{11}g_{22} - g_{12}^2.
     350
     351        See also ``area_form``.
     352       
     353        OUTPUT:
     354       
     355         - Square of the area form
     356
     357        EXAMPLES::
     358
     359            sage: u, v = var('u, v')
     360            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     361            sage: sphere.area_form_squared()
     362            cos(v)^2
     363
     364        """
     365        gg = self.first_fundamental_form_coefficients()
     366        return  (gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2).simplify_full()
     367
     368
     369    @cached_method
     370    def area_form(self):
     371        r"""
     372        Returns the coefficient of the area form on the surface.  In terms of the
     373        coefficients $g_{ij}$ (where $i, j = 1, 2$) of the first fundamental form, the
     374        coefficient of the area form is given by
     375
     376        .. math::
     377
     378           A = \sqrt{g_{11}g_{22} - g_{12}^2}.
     379
     380        See also ``area_form_squared``.
     381       
     382        OUTPUT: 
     383       
     384         - Coefficient of the area form
     385
     386        EXAMPLES::
     387
     388            sage: u, v = var('u,v')
     389            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     390            sage: sphere.area_form()
     391            abs(cos(v))
     392
     393        """
     394        return sqrt(self.area_form_squared()).simplify_full()
     395
     396
     397    @cached_method
     398    def first_fundamental_form_inverse_coefficients(self):
     399        r"""
     400        Returns the coefficients $g^{ij}$ of the inverse of the fundamental form,
     401        as a dictionary.  The inverse coefficients are defined by
     402
     403         .. math::
     404
     405               \sum_j g^{ij} g_{jk} = \delta^i_k
     406
     407        with $\delta^i_k$ the Kronecker delta.
     408
     409        OUTPUT:
     410 
     411         - Dictionary of the inverse coefficients.
     412
     413        EXAMPLES::
     414
     415            sage: u, v = var('u, v')
     416            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     417            sage: sphere.first_fundamental_form_inverse_coefficients()
     418            {(1, 2): 0, (1, 1): cos(v)^(-2), (2, 1): 0, (2, 2): 1}
     419           
     420        """
     421
     422        gg = self.first_fundamental_form_coefficients()
     423        DD = gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2
     424       
     425        gi11 = (gg[(2,2)]/DD).simplify_full()
     426        gi12 = (-gg[(1,2)]/DD).simplify_full()
     427        gi21 = gi12
     428        gi22 = (gg[(1,1)]/DD).simplify_full()
     429
     430        return {(1,1): gi11, (1,2): gi12, (2,1): gi21, (2,2): gi22}
     431
     432
     433    def first_fundamental_form_inverse_coefficient(self, index):
     434        r"""
     435        Returns a specific component $g^{ij}$ of the inverse of the fundamental form.
     436
     437        INPUT:
     438
     439         - ``index`` - index ``(i, j)`` of the component ``g^{ij}``.
     440
     441        OUTPUT:
     442
     443         - Component of the inverse of the fundamental form.
     444
     445        EXAMPLES::
     446
     447            sage: u, v = var('u, v')
     448            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     449            sage: sphere.first_fundamental_form_inverse_coefficient((1, 2))
     450            0
     451            sage: sphere.first_fundamental_form_inverse_coefficient((1, 1))
     452            cos(v)^(-2)
     453
     454        """
     455
     456        index = tuple(sorted(index))
     457        if index not in self.index_list:
     458            raise ValueError, "Index %s out of bounds." % str(index)
     459        return self.first_fundamental_form_inverse_coefficients()[index]
     460
     461
     462    @cached_method
     463    def rotation(self,theta):
     464        r"""
     465        Gives the matrix of the rotation operator over a given angle $\theta$ with respect to the natural frame.
     466       
     467        INPUT:
     468       
     469         - ``theta`` - rotation angle
     470
     471        OUTPUT: 
     472
     473         - Rotation matrix with respect to the natural frame.
     474
     475        ALGORITHM:
     476
     477        The operator of rotation over $\pi/2$ is $J^i_j = g^{ik}\omega_{jk}$, where $\omega$ is the area form.  The operator of rotation over an angle $\theta$ is $\cos(\theta) I + sin(\theta) J$
     478
     479        EXAMPLES::
     480
     481            sage: u, v = var('u, v')
     482            sage: assume(cos(v)>0)
     483            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     484
     485        We first compute the matrix of rotation over pi/3::
     486
     487            sage: rotation = sphere.rotation(pi/3); rotation
     488            [                1/2 -1/2*sqrt(3)/cos(v)]
     489            [ 1/2*sqrt(3)*cos(v)                 1/2]
     490
     491        We verify that three succesive rotations over pi/3 yield minus the identity::
     492
     493            sage: rotation^3
     494            [-1  0]
     495            [ 0 -1]
     496           
     497        """
     498
     499        from sage.functions.trig import sin, cos
     500
     501        gi = self.first_fundamental_form_inverse_coefficients()
     502        w12 = self.area_form()
     503        RR11 = (cos(theta) + sin(theta)*gi[1,2]*w12).simplify_full()
     504        RR12 = (- sin(theta)*gi[1,1]*w12).simplify_full()
     505        RR21 = (sin(theta)*gi[2,2]*w12).simplify_full()
     506        RR22 = (cos(theta) - sin(theta)*gi[2,1]*w12).simplify_full()       
     507        rotation = matrix([[RR11,RR12],[RR21,RR22]])
     508        return rotation
     509
     510
     511
     512    @cached_method
     513    def orthonormal_frame(self, coordinates='ext'):
     514        r"""
     515        Returns the orthonormal frame field on the surface, expressed either in exterior
     516        coordinates (i.e. expressed as vector fields in the ambient space $\mathbb{R}^3$,
     517        the default) or interior coordinates (with respect to the natural frame)
     518       
     519        INPUT:
     520
     521         - ``coordinates`` - either ``ext`` (default) or ``int``.
     522
     523        OUTPUT:
     524       
     525         - Orthogonal frame field as a dictionary.
     526
     527        ALGORITHM:
     528
     529        We normalize the first vector $\vec e_1$ of the natural frame and then get the second frame vector as $\vec e_2 = [\vec n, \vec e_1]$, where $\vec n$ is the unit normal to the surface.
     530
     531        EXAMPLES::
     532
     533            sage: u, v = var('u,v')
     534            sage: assume(cos(v)>0)
     535            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v), sin(u)*cos(v), sin(v)], [u, v],'sphere')
     536            sage: frame = sphere.orthonormal_frame(); frame
     537            {1: (-sin(u), cos(u), 0), 2: (-sin(v)*cos(u), -sin(u)*sin(v), cos(v))}
     538            sage: (frame[1]*frame[1]).simplify_full()
     539            1
     540            sage: (frame[1]*frame[2]).simplify_full()
     541            0
     542            sage: frame[1] == sphere.orthonormal_frame_vector(1)
     543            True
     544
     545        We compute the orthonormal frame with respect to the natural frame on the surface::
     546       
     547            sage: frame_int = sphere.orthonormal_frame(coordinates='int'); frame_int
     548            {1: (1/cos(v), 0), 2: (0, 1)}
     549            sage: sphere.first_fundamental_form(frame_int[1], frame_int[1])
     550            1
     551            sage: sphere.first_fundamental_form(frame_int[1], frame_int[2])
     552            0
     553            sage: sphere.first_fundamental_form(frame_int[2], frame_int[2])
     554            1
     555
     556        """
     557
     558
     559        from sage.symbolic.constants import pi
     560
     561        if coordinates not in ['ext', 'int']:
     562            raise ValueError, \
     563                  r"Coordinate system must be exterior ('ext') or interior ('int')."
     564
     565        c  = self.first_fundamental_form_coefficient([1,1])
     566        if coordinates == 'ext':
     567            f1 = self.natural_frame()[1]
     568
     569            E1 = simplify_vector(f1/sqrt(c))
     570            E2 = simplify_vector(self.normal_vector(normalized=True).cross_product(E1))
     571        else:
     572            E1 =  vector([(1/sqrt(c)).simplify_full(), 0])
     573            E2 = simplify_vector(self.rotation(pi/2)*E1)
     574        return  {1:E1, 2:E2}
     575
     576
     577    @cached_method
     578    def orthonormal_frame_vector(self, index, coordinates='ext'):
     579        r"""
     580        Returns a specific basis vector field of the orthonormal frame field on the surface, expressed in exterior or interior coordinates.  See ``orthogonal_frame`` for more details.
     581       
     582        INPUT:
     583
     584         - ``index`` - index ``(i, j)`` of the basis vector;
     585         - ``coordinates`` - either ``ext`` (default) or ``int``.
     586
     587        OUTPUT:
     588
     589         - Orthonormal frame vector field.
     590
     591        EXAMPLES::
     592
     593            sage: u, v = var('u, v')
     594            sage: assume(cos(v)>0)
     595            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     596            sage: V1 = sphere.orthonormal_frame_vector(1); V1
     597            (-sin(u), cos(u), 0)
     598            sage: V2 = sphere.orthonormal_frame_vector(2); V2
     599            (-sin(v)*cos(u), -sin(u)*sin(v), cos(v))
     600            sage: (V1*V1).simplify_full()
     601            1
     602            sage: (V1*V2).simplify_full()
     603            0
     604       
     605            sage: from sage.riemann.vector_functions import simplify_vector
     606            sage: n = sphere.normal_vector(normalized=True)
     607            sage: simplify_vector(V1.cross_product(V2) - n)
     608            (0, 0, 0)
     609        """
     610       
     611        return self.orthonormal_frame(coordinates)[index]
     612
     613       
     614    def lie_bracket(self, v, w):
     615        """
     616        Returns the Lie bracket of two vector fields that are tangent to the surface.
     617        The vector fields should be given in coordinates with respect to the natural frame.
     618         
     619        INPUT:
     620
     621         - ``v`` and ``w`` - vector fields on the surface, expressed as tuples of functions.
     622               
     623        OUTPUT:
     624
     625         - The Lie bracket ``[v, w]``
     626
     627
     628        EXAMPLES::
     629
     630            sage: u, v = var('u, v')
     631            sage: assume(cos(v)>0)
     632            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     633            sage: sphere.lie_bracket([u,v],[-v,u])
     634            (0, 0)
     635
     636            sage: EE_int = sphere.orthonormal_frame(coordinates='int')
     637            sage: sphere.lie_bracket(EE_int[1],EE_int[2])
     638            (sin(v)/cos(v)^2, 0)
     639        """
     640        vv = vector([xx for xx in v])
     641        ww = vector([xx for xx in w])
     642        uuu = self.variables_list
     643        Dvv = matrix([ [ diff(xxx,uu).simplify_full() for uu in uuu ] for xxx in vv])   
     644        Dww = matrix([ [ diff(xxx,uu).simplify_full() for uu in uuu ] for xxx in ww])   
     645        return simplify_vector(vector(Dvv*ww - Dww*vv))
     646
     647
     648    def frame_structure_functions(self, e1, e2):
     649        r"""
     650        Returns the structure functions $c^k_{ij}$ for a frame field $e_1, e_2$, i.e.
     651        a pair of vector fields on the surface which are linearly independent at
     652        each point.  The structure functions are defined using the Lie bracket by
     653
     654         .. math::
     655
     656              [e_i,e_j] = c^k_{ij}e_k.
     657         
     658        INPUT:
     659
     660         - ``e1``, ``e2`` - vector fields on the surface, expressed as pairs of functions.
     661
     662        OUTPUT:
     663
     664        Dictionary of structure functions, where the key (i, j, k) refers to the function
     665        $c^k_{ij}$.
     666
     667        EXAMPLES::
     668
     669            sage: u, v = var('u, v')
     670            sage: assume(cos(v) > 0)
     671            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     672            sage: sphere.frame_structure_functions([u,v],[-v,u])
     673            {(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}
     674
     675        We construct the structure functions of the orthonormal frame on the surface::
     676           
     677            sage: EE_int = sphere.orthonormal_frame(coordinates='int')
     678            sage: CC = sphere.frame_structure_functions(EE_int[1],EE_int[2]); CC
     679            {(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}
     680            sage: sphere.lie_bracket(EE_int[1],EE_int[2]) - CC[(1,2,1)]*EE_int[1] - CC[(1,2,2)]*EE_int[2]
     681            (0, 0)
     682            """
     683        ee1 = vector([xx for xx in e1])
     684        ee2 = vector([xx for xx in e2])
     685       
     686        lb = simplify_vector(self.lie_bracket(ee1,ee2))
     687        trans_matrix = matrix([[ee1[0],ee2[0]],[ee1[1],ee2[1]]])
     688        zz = simplify_vector(trans_matrix.inverse()*lb)
     689        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}
     690
     691
     692    @cached_method
     693    def _compute_second_order_frame_element(self, index):
     694        """
     695        Compute an element of the second order frame of the surface.  See
     696        ``second_order_natural_frame`` for more details.
     697
     698        This method expects its arguments in tuple form for caching.  As it does no input
     699        checking, it should not be called directly.
     700
     701        EXAMPLES::
     702
     703            sage: u, v = var('u, v')
     704            sage: paraboloid = ParametrizedSurface3D([u, v, u^2 + v^2], [u,v], 'paraboloid')
     705            sage: paraboloid._compute_second_order_frame_element((1, 2))
     706            (0, 0, 0)
     707            sage: paraboloid._compute_second_order_frame_element((2, 2))
     708            (0, 0, 2)
     709
     710        """
     711        variables = [self.variables[i] for i in index]
     712        ddr_element = vector([diff(f, variables).simplify_full() for f in self.equation])
     713       
     714        return ddr_element
     715
     716
     717    @cached_method
     718    def second_order_natural_frame(self):
     719        r"""
     720        Returns the second-order frame of the surface, i.e. computes the second-order derivatives
     721        (with respect to the parameters on the surface) of the parametric expression
     722        $\vec r = \vec r(u^1,u^2)$ of the surface.
     723         
     724        OUTPUT:
     725
     726         - Dictionary where the keys are 2-tuples ``(i, j)`` and the values are the corresponding derivatives ``r_ij``.
     727
     728        EXAMPLES:
     729
     730        We compute the second-order natural frame of the sphere::
     731
     732            sage: u, v = var('u, v')
     733            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     734            sage: sphere.second_order_natural_frame()
     735            {(1, 2): (sin(u)*sin(v), -sin(v)*cos(u), 0), (1, 1): (-cos(u)*cos(v),
     736            -sin(u)*cos(v), 0), (2, 1): (sin(u)*sin(v), -sin(v)*cos(u), 0), (2, 2):
     737            (-cos(u)*cos(v), -sin(u)*cos(v), -sin(v))}
     738
     739        """
     740
     741        vectors = {}
     742        for index in self.index_list:
     743            sorted_index = tuple(sorted(index))
     744            vectors[index] = \
     745                self._compute_second_order_frame_element(sorted_index)
     746        return vectors
     747
     748
     749    def second_order_natural_frame_element(self, index):
     750        r"""
     751        Returns a vector in the second-order frame of the surface, i.e. computes the
     752        second-order derivatives of the parametric expression $\vec{r}$ of the
     753        surface with respect to the parameters listed in the argument.
     754
     755        INPUT:
     756
     757         - ``index`` - a 2-tuple ``(i, j)`` specifying the element of the second-order frame.
     758         
     759        OUTPUT:
     760
     761         - The second-order derivative ``r_ij``.
     762
     763        EXAMPLES::
     764
     765            sage: u, v = var('u, v')
     766            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     767            sage: sphere.second_order_natural_frame_element((1, 2))
     768            (sin(u)*sin(v), -sin(v)*cos(u), 0)
     769           
     770        """
     771       
     772        index = tuple(sorted(index))
     773        if index not in self.index_list:
     774            raise ValueError, "Index %s out of bounds." % str(index)
     775        return self._compute_second_order_frame_element(index)
     776
     777
     778    @cached_method
     779    def _compute_second_fundamental_form_coefficient(self, index):
     780        """
     781        Compute a coefficient of the second fundamental form of the surface.  See
     782        ``second_fundamental_form_coefficient`` for more details.
     783
     784        This method expects its arguments in tuple form for caching.  As it does no input
     785        checking, it should not be called directly.
     786
     787        EXAMPLES::
     788
     789            sage: u, v = var('u,v')
     790            sage: paraboloid = ParametrizedSurface3D([u, v, u^2+v^2], [u, v], 'paraboloid')
     791            sage: paraboloid._compute_second_fundamental_form_coefficient((1,1))
     792            2/sqrt(4*u^2 + 4*v^2 + 1)
     793
     794        """
     795        NN = self.normal_vector(normalized=True)
     796        v  = self.second_order_natural_frame_element(index)
     797        return (v*NN).simplify_full()
     798
     799       
     800    def second_fundamental_form_coefficient(self, index):
     801        r"""
     802        Returns the coefficient $h_{ij}$ of the second fundamental form corresponding to the
     803        index $(i, j)$.  If the equation of the surface is $\vec{r}(u^1, u^2)$, then
     804
     805         .. math::
     806
     807            h_{ij} = \vec{r}_{u^i u^j} \cdot \vec{n},
     808
     809        where $\vec{n}$ is the unit normal.
     810
     811        INPUT:
     812
     813         - ``index`` - a 2-tuple ``(i, j)``
     814
     815        OUTPUT:
     816
     817         - Component ``h_{ij}`` of the second fundamental form.
     818
     819        EXAMPLES::
     820       
     821            sage: u, v = var('u,v')
     822            sage: assume(cos(v)>0)
     823            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     824            sage: sphere.second_fundamental_form_coefficient((1, 1))
     825            -cos(v)^2
     826            sage: sphere.second_fundamental_form_coefficient((2, 1))
     827            0
     828
     829        """
     830        index = tuple(index)
     831        if index not in self.index_list:
     832            raise ValueError, "Index %s out of bounds." % str(index)
     833        return self._compute_second_fundamental_form_coefficient(index)
     834
     835
     836    @cached_method
     837    def second_fundamental_form_coefficients(self):
     838        """
     839        Returns the coefficients $h_{ij}$ of the second fundamental form as a dictionary, where the
     840        keys are the indices $(i, j)$ and the values are the corresponding components $h_{ij}$.
     841
     842        When only one component is needed, consider the function ``second_fundamental_form_coefficient``.
     843       
     844        OUTPUT:
     845
     846        Dictionary of second fundamental form coefficients.
     847           
     848        EXAMPLES::
     849
     850            sage: u, v = var('u, v')
     851            sage: assume(cos(v)>0)
     852            sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     853            sage: sphere.second_fundamental_form_coefficients()
     854            {(1, 2): 0, (1, 1): -cos(v)^2, (2, 1): 0, (2, 2): -1}
     855
     856        """
     857
     858        coefficients = {}
     859        for index in self.index_list:
     860            coefficients[index] = \
     861                self._compute_second_fundamental_form_coefficient(index)
     862        return coefficients
     863
     864
     865    def second_fundamental_form(self,vector1,vector2):
     866        r"""
     867        Evaluates the second fundamental form on two vectors on the surface. If the
     868        vectors are given by $v=(v^1,v^2)$ and $w=(w^1,w^2)$, the result of this function
     869        is $h_{11} v^1 w^1 + h_{12}(v^1 w^2 + v^2 w^1) + h_{22} v^2 w^2$.
     870
     871        INPUT:
     872
     873         - ``vector1``, ``vector2`` - 2-tuples representing the input vectors.
     874
     875        OUTPUT:
     876
     877         - Value of the second fundamental form evaluated on the given vectors.
     878
     879        EXAMPLES:
     880
     881        We evaluate the second fundamental form on two symbolic vectors::
     882
     883           sage: u, v = var('u, v')
     884           sage: v1, v2, w1, w2 = var('v1, v2, w1, w2')
     885           sage: assume(cos(v) > 0)
     886           sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     887           sage: sphere.second_fundamental_form(vector([v1, v2]), vector([w1, w2]))
     888           -v1*w1*cos(v)^2 - v2*w2
     889
     890        We evaluate the second fundamental form on vectors with numerical components::
     891           
     892           sage: vect = vector([1,2])
     893           sage: sphere.second_fundamental_form(vect, vect)
     894           -cos(v)^2 - 4
     895           sage: sphere.second_fundamental_form([1,1], [2,1])
     896           -2*cos(v)^2 - 1
     897                   
     898        """
     899        hh = self.second_fundamental_form_coefficients()
     900        return sum(hh[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list)
     901
     902
     903    @cached_method
     904    def gauss_curvature(self):
     905        r"""
     906        Finds the gaussian curvature of the surface, given by $K = \frac{h_{11}h_{22} - h_{12}^2}{g_{11}g_{22} - g_{12}^2}$, where $g_{ij}$ and $h_{ij}$ are the coefficients of the first and second fundamental form, respectively.
     907
     908        OUTPUT:
     909
     910         - Gaussian curvature of the surface.
     911
     912        EXAMPLES::
     913
     914           sage: R = var('R')
     915           sage: assume(R>0)
     916           sage: u, v = var('u,v')
     917           sage: assume(cos(v)>0)
     918           sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     919           sage: sphere.gauss_curvature()
     920           R^(-2)
     921                   
     922        """
     923        hh = self.second_fundamental_form_coefficients()
     924        return ((hh[(1,1)]*hh[(2,2)]-hh[(1,2)]**2)/self.area_form_squared()).simplify_full()
     925
     926
     927    @cached_method
     928    def mean_curvature(self):
     929        r"""
     930        Finds the mean curvature of the surface, given by
     931        $H = \frac{1}{2}\frac{g_{22}h_{11} - 2g_{12}h_{12} + g_{11}h_{22}}{g_{11}g_{22} - g_{12}^2}$,
     932        where $g_{ij}$ and $h_{ij}$ are the components of the first and second fundamental forms,
     933        respectively.
     934
     935        OUTPUT:
     936
     937         - Mean curvature of the surface
     938
     939        EXAMPLES::
     940           
     941           sage: R = var('R')
     942           sage: assume(R>0)
     943           sage: u, v = var('u,v')
     944           sage: assume(cos(v)>0)
     945           sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     946           sage: sphere.mean_curvature()
     947           -1/R
     948
     949        """
     950        gg = self.first_fundamental_form_coefficients()
     951        hh = self.second_fundamental_form_coefficients()
     952        denom = 2*self.area_form_squared()
     953        enum =(gg[(2,2)]*hh[(1,1)]-2*gg[(1,2)]*hh[(1,2)]+gg[(1,1)]*hh[(2,2)]).simplify_full()
     954        return (enum/denom).simplify_full()
     955
     956
     957    @cached_method
     958    def principal_curvatures(self):
     959        """
     960        Finds the principal curvatures of the surface
     961
     962        OUTPUT:
     963        The two principal curvatures, given as a dictionary.
     964
     965        EXAMPLES::
     966                     
     967           sage: R = var('R')
     968           sage: assume(R>0)
     969           sage: u, v = var('u,v')
     970           sage: assume(cos(v)>0)
     971           sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     972           sage: sphere.principal_curvatures()
     973           {1: -1/R, 2: -1/R}
     974
     975           sage: u, v = var('u,v')
     976           sage: R, r = var('R, r')
     977           sage: assume(R>r,r>0)
     978           sage: torus = ParametrizedSurface3D([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
     979           sage: torus.principal_curvatures()
     980           {1: -cos(v)/(r*cos(v) + R), 2: -1/r}
     981
     982        """
     983
     984        from sage.symbolic.assumptions import assume
     985        from sage.symbolic.relation import solve
     986        from sage.calculus.var import var
     987
     988        KK = self.gauss_curvature()
     989        HH = self.mean_curvature()
     990
     991        # jvkersch: when this assumption is uncommented, Sage raises an error stating that the assumption
     992        # is redundant... Can we safely omit this, based on some geometric reasoning?
     993       
     994        # assume(HH**2-KK>=0)
     995        # mikarm: This is a problem, I had a lot of trobles here. Of course, this inequality always hold true.
     996        # I insert this assumption because sage sometimes, in simplification, uses the complex numbers, though the roots
     997        # are, for sure, real. 
     998        # This, in turn, causes problems when we substitute coordinates into the expression of principal curvatures.
     999        # Unfortunately, I did not manage to tell Sage that they are real (to declare the variables as real).
     1000        # Moreover, in a neighborhood of an umbilic point we even cannot "smoothly" order the set of principal curvatures.
     1001        # So, in general, at present this method is far from the final form.
     1002         
     1003
     1004 
     1005        x = var('x')
     1006        sol = solve(x**2 -2*HH*x + KK == 0,x)
     1007       
     1008        #k1=var('k1')
     1009        #k2=var('k2')
     1010
     1011        # jvkersch: when I tried to run the previous version of the code, I ran into the problem that if the equation for the principal curvatures had a double root (as in the case of the sphere example in the worksheet), solve returned only one root.  Maybe this is a difference due to having different versions of sage.
     1012
     1013        k1 = (x.substitute(sol[0])).simplify_full()
     1014        if len(sol) == 1:
     1015            k2 = k1
     1016        else:
     1017            k2 = (x.substitute(sol[1])).simplify_full()
     1018       
     1019        return {1:k1,2: k2}
     1020
     1021
     1022    @cached_method
     1023    def shape_operator_coefficients(self):
     1024        r"""
     1025        Returns the components of the shape operator of the surface as a dictionary.
     1026        See ``shape_operator`` for more information.
     1027
     1028        OUTPUT:
     1029
     1030         - Dictionary where the keys are two-tuples ``(i, j)``, with values the
     1031           corresponding component of the shape operator.
     1032
     1033        EXAMPLES::
     1034
     1035           sage: R = var('R')
     1036           sage: u, v = var('u,v')
     1037           sage: assume(cos(v)>0)
     1038           sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     1039           sage: sphere.shape_operator_coefficients()
     1040           {(1, 2): 0, (1, 1): -1/R, (2, 1): 0, (2, 2): -1/R}
     1041
     1042        """
     1043
     1044        gi = self.first_fundamental_form_inverse_coefficients()
     1045        hh = self.second_fundamental_form_coefficients()
     1046
     1047        shop11 = (gi[(1,1)]*hh[(1,1)] + gi[(1,2)]*hh[(1,2)]).simplify_full()
     1048        shop12 = (gi[(1,1)]*hh[(2,1)] + gi[(1,2)]*hh[(2,2)]).simplify_full()
     1049        shop21 = (gi[(2,1)]*hh[(1,1)] + gi[(2,2)]*hh[(1,2)]).simplify_full()
     1050        shop22 = (gi[(2,1)]*hh[(2,1)] + gi[(2,2)]*hh[(2,2)]).simplify_full()
     1051
     1052        return {(1,1): shop11, (1,2): shop12, (2,1): shop21, (2,2): shop22}
     1053
     1054
     1055
     1056    # jvkersch: I've changed the definition of the shape operator to return the matrix
     1057    # of the shape operator, rather than a matrix times a vector.  This will make it easier
     1058    # to compute eigenvalues and eigenvectors of the shape operator.
     1059   
     1060    def shape_operator(self):
     1061        r"""
     1062        Returns the shape operator of the surface as a matrix.  The shape operator is
     1063        defined as the derivative of the Gauss map, and is computed here in terms of the
     1064        first and second fundamental form by means of the Weingarten equations.
     1065
     1066        OUTPUT:
     1067
     1068         - Matrix of the shape operator
     1069
     1070        EXAMPLES::
     1071       
     1072           sage: R = var('R')
     1073           sage: assume(R>0)
     1074           sage: u, v = var('u,v')
     1075           sage: assume(cos(v)>0)
     1076           sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     1077           sage: S = sphere.shape_operator(); S
     1078           [-1/R    0]
     1079           [   0 -1/R]
     1080
     1081        The eigenvalues of the shape operator are the principal curvatures of the surface::
     1082       
     1083            sage: u, v = var('u,v')
     1084            sage: paraboloid = ParametrizedSurface3D([u, v, u^2+v^2], [u, v], 'paraboloid')
     1085            sage: S = paraboloid.shape_operator(); S
     1086            [2*(4*v^2 + 1)/(4*u^2 + 4*v^2 + 1)^(3/2)        -8*u*v/(4*u^2 + 4*v^2 + 1)^(3/2)]
     1087            [       -8*u*v/(4*u^2 + 4*v^2 + 1)^(3/2) 2*(4*u^2 + 1)/(4*u^2 + 4*v^2 + 1)^(3/2)]
     1088            sage: S.eigenvalues()
     1089            [2*sqrt(4*u^2 + 4*v^2 + 1)/(16*u^4 + 8*u^2 + 16*v^4 + 8*(4*u^2 + 1)*v^2 + 1), 2/sqrt(4*u^2 + 4*v^2 + 1)]
     1090            sage: paraboloid.principal_curvatures()
     1091            {1: 2/(4*u^2 + 4*v^2 + 1)^(3/2), 2: 2/sqrt(4*u^2 + 4*v^2 + 1)}
     1092
     1093        """
     1094       
     1095        #vv = vector([xx for xx in v])
     1096        shop = self.shape_operator_coefficients()
     1097        shop_matrix=matrix([[shop[(1,1)],shop[(1,2)]],[shop[(2,1)],shop[(2,2)]]])
     1098        return simplify_matrix(shop_matrix)
     1099
     1100
     1101    @cached_method
     1102    def principal_directions(self):
     1103        """
     1104        Finds the principal curvatures and principal directions of the surface
     1105
     1106        OUTPUT:
     1107
     1108         - The dictionary of lists [principal curvature, corresponding principal direction]         
     1109
     1110        If principal curvatures coincide, gives the warning that the surface is a sphere.
     1111
     1112        EXAMPLES::
     1113
     1114           sage: u, v = var('u, v')
     1115           sage: R, r = var('R,r')
     1116           sage: assume(R>r,r>0)
     1117           sage: torus = ParametrizedSurface3D([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
     1118           sage: pdd = torus.principal_directions()
     1119           sage: pdd[1]
     1120           [-cos(v)/(r*cos(v) + R), (1, 0)]
     1121           sage: pdd[2]
     1122           [-1/r, (0, -(R*r*cos(v) + R^2)/r)]
     1123
     1124           sage: R = var('R')
     1125           sage: assume(R>0)
     1126           sage: assume(cos(v)>0)
     1127           sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     1128           sage: sphere.principal_directions()
     1129           'This is a sphere, so any direction is principal'
     1130
     1131           sage: catenoid = ParametrizedSurface3D([R*cosh(v)*cos(u), R*cosh(v)*sin(u),v],[u,v],'catenoid')
     1132           sage: pd = catenoid.principal_directions()
     1133           sage: pd[1][1]
     1134           (1, 0)
     1135           sage: pd[2][1]
     1136           (0, 1/2*(sqrt(R^2*sinh(v)^2 + 1)*sqrt((4*sinh(v)^2*cosh(v)^2 + 1)*R^4 + 2*(2*cosh(v)^2 - 1)*R^2 + 1)*R*cosh(v) + sqrt(R^2*sinh(v)^2 + 1)*((2*sinh(v)^2*cosh(v) + cosh(v))*R^3 + R*cosh(v)))/(R^4*sinh(v)^4 + 2*R^2*sinh(v)^2 + 1))
     1137           sage: pd[1][1]*pd[2][1]
     1138           0
     1139        """
     1140        gg = self.first_fundamental_form_coefficients()
     1141        hh = self.second_fundamental_form_coefficients()
     1142        kk = self.principal_curvatures()
     1143        if kk[1]==kk[2]:
     1144            return "This is a sphere, so any direction is principal"
     1145        pd1 = simplify_vector([hh[(1,2)]-kk[1]*gg[(1,2)],-hh[(1,1)]+kk[1]*gg[(1,1)]])
     1146        if pd1==vector([0,0]):
     1147           pd1 = vector([1,0])
     1148        pd2 = simplify_vector([hh[(1,2)]-kk[2]*gg[(1,2)],-hh[(1,1)]+kk[2]*gg[(1,1)]])
     1149        if pd2==vector([0,0]):
     1150           pd2 = vector([1,0])
     1151        return {1:[kk[1],pd1],2:[kk[2],pd2]}
     1152
     1153
     1154    ### Jvkersch: doctests checked up to this point.
     1155
     1156    # jvkersch: alternative definition of principal_directions.  The behavior of this
     1157    # function is consistent with the eigenvector  function in Sage (since it is just
     1158    # looking for eigenvectors of the shape operator).  In particular, for a spherical
     1159    # surface it just returns an arbitrary pair of vectors.
     1160    def principal_directions_new(self):
     1161        r"""
     1162        Finds the principal curvatures and principal directions of the surface
     1163
     1164        OUTPUT:
     1165
     1166        For each principal curvature, returns a list of the form ``(rho, V, n)``,
     1167        where ``rho`` is the principal curvature, ``V`` is the corresponding
     1168        principal direction, and ``n`` is the multiplicity.
     1169
     1170        EXAMPLES::
     1171
     1172           sage: u, v = var('u, v')
     1173           sage: R, r = var('R,r')
     1174           sage: assume(R>r,r>0)
     1175           sage: torus = ParametrizedSurface3D([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
     1176           sage: torus.principal_directions_new()
     1177           [(-cos(v)/(r*cos(v) + R), [(1, 0)], 1), (-1/r, [(0, 1)], 1)]
     1178
     1179        """
     1180        return self.shape_operator().eigenvectors_left()
     1181   
     1182
     1183    @cached_method
     1184    def connection_coefficients(self):
     1185        r"""
     1186        Computes the connection coefficients or Christoffel symbols $\Gamma^k_{ij}$ of the surface.
     1187        If the coefficients of the first fundamental form are given by $g_{ij}$
     1188        (where $i, j = 1, 2$), then
     1189       
     1190         .. math::
     1191
     1192           \Gamma^k_{ij} = \frac{1}{2} g^{kl} \left( \frac{\partial g_{li}}{\partial x^j}
     1193             - \frac{\partial g_{ij}}{\partial x^l}
     1194             + \frac{\partial g_{lj}}{\partial x^i} \right).
     1195
     1196        Here, $(g^{kl})$ is the inverse of the matrix $(g_{ij})$, with $i, j = 1, 2$.
     1197
     1198        OUTPUT:
     1199       
     1200        Dictionary of connection coefficients, where the keys are 3-tuples $(i,j,k)$ and the values are the corresponding coefficients $\Gamma^k_{ij}$. 
     1201
     1202        EXAMPLES::
     1203                                 
     1204           sage: r = var('r')
     1205           sage: assume(r > 0)
     1206           sage: u, v = var('u,v')
     1207           sage: assume(cos(v)>0)
     1208           sage: sphere = ParametrizedSurface3D([r*cos(u)*cos(v),r*sin(u)*cos(v),r*sin(v)],[u,v],'sphere')
     1209           sage: sphere.connection_coefficients()
     1210           {(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}
     1211           
     1212        """
     1213        x = self.variables
     1214        gg = self.first_fundamental_form_coefficients()
     1215        gi = self.first_fundamental_form_inverse_coefficients()
     1216
     1217        dg = {}
     1218        for i,j,k in self.index_list_3:
     1219            dg[(i,j,k)]=gg[(j,k)].differentiate(x[i]).simplify_full()
     1220           
     1221        structfun={}
     1222        for i,j,k in self.index_list_3:
     1223            structfun[(i,j,k)]=sum(gi[(k,s)]*(dg[(i,j,s)]+dg[(j,i,s)]-dg[(s,i,j)])/2 for s in (1,2)).full_simplify()
     1224        return structfun
     1225
     1226   
     1227    @cached_method
     1228    def _create_geodesic_ode_system(self):
     1229        r"""
     1230        Helper method to create a fast floating-point version of the geodesic equations,
     1231        used by ``geodesics_numerical``.
     1232
     1233        EXAMPLES::
     1234           sage: p, q = var('p,q')
     1235           sage: sphere = ParametrizedSurface3D([cos(q)*cos(p),sin(q)*cos(p),sin(p)],[p,q],'sphere')
     1236           sage: ode = sphere._create_geodesic_ode_system()
     1237           sage: ode.function(0.0, (1.0, 0.0, 1.0, 1.0))
     1238           [1.00000000000000, 1.00000000000000, -0.45464871341284091, 3.1148154493098041]
     1239           
     1240        """
     1241        from sage.ext.fast_eval import fast_float
     1242        from sage.calculus.var import var
     1243        from sage.gsl.ode import ode_solver
     1244
     1245        u1 = self.variables[1]
     1246        u2 = self.variables[2]
     1247        v1, v2 = var('v1, v2')
     1248
     1249        C = self.connection_coefficients()
     1250       
     1251        dv1 = - C[(1,1,1)]*v1**2 - 2*C[(1,2,1)]*v1*v2 - C[(2,2,1)]*v2**2
     1252        dv2 = - C[(1,1,2)]*v1**2 - 2*C[(1,2,2)]*v1*v2 - C[(2,2,2)]*v2**2
     1253        fun1 = fast_float(dv1, str(u1), str(u2), str(v1), str(v2))
     1254        fun2 = fast_float(dv2, str(u1), str(u2), str(v1), str(v2))
     1255
     1256        geodesic_ode = ode_solver()
     1257        geodesic_ode.function = \
     1258                              lambda t, (u1, u2, v1, v2) : \
     1259                              [v1, v2, fun1(u1, u2, v1, v2), fun2(u1, u2, v1, v2)]
     1260        return geodesic_ode
     1261
     1262       
     1263    def geodesics_numerical(self, p0, v0, tinterval):
     1264        r"""
     1265        Numerical integration of the geodesic equations.  Explicitly, the
     1266        geodesic equations are given by
     1267
     1268         .. math::
     1269
     1270           \frac{d^2 u^i}{dt^2} + \Gamma^i_{jk} \frac{d u^j}{dt} \frac{d u^k}{dt} = 0.
     1271
     1272        Solving these equations gives the coordinates $(u^1, u^2)$ on the surface
     1273        of the geodesic.  The coordinates in space can then be found by substituting
     1274        $(u^1, u^2)$ into the vector representating $\vec{r}(u^1, u^2)$ of the surface.
     1275       
     1276        ALGORITHM:
     1277
     1278        The geodesic equations are integrated forward in time using the
     1279        ode solvers from ``sage.gsl.ode``.  See ``_create_geodesic_ode_system``
     1280        for more details.
     1281
     1282        INPUT:
     1283
     1284         - ``p0`` - 2-tuple with coordinates of the initial point.
     1285
     1286         - ``v0`` - 2-tuple with components of the initial tangent vector to the geodesic.
     1287
     1288         - ``tinterval`` - List ``[a,b,M]``, where ``(a,b)`` is the domain of the geodesic and ``M`` is the number of division points
     1289
     1290        OUTPUT:
     1291
     1292        List ``[t, [u1(t),u2(t)], [v1(t),v2(t)], [x1(t),x2(t),x3(t)]]``, where
     1293
     1294         - ``t`` is a subdivision point;
     1295
     1296         - ``[u1(t),u2(t)]`` is the list of coordinates of the geodesic point;
     1297         
     1298         - ``[v1(t),v2(t)]`` is the list of coordinates of the vector tanget to the geodesic; 
     1299
     1300         - ``[x1(t),x2(t),x3(t)]`` is the list of coordinates of the geodesic point in the three-dimensional space.
     1301
     1302        EXAMPLES::
     1303
     1304           sage: p, q = var('p,q')
     1305           sage: v = [p,q]
     1306           sage: assume(cos(q)>0)
     1307           sage: sphere = ParametrizedSurface3D([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
     1308           sage: gg_array = sphere.geodesics_numerical([0.0,0.0],[1.0,1.0],[0,2*pi,5])
     1309           sage: gg_array[0]
     1310           [0, [0.000000000000000, 0.000000000000000], [1.00000000000000, 1.00000000000000], [1, 0, 0]]
     1311           sage: gg_array[1]
     1312           [1.2566370614359172, [0.76440092815484362, 1.8586224292405702], [-0.28386842687533731, 1.9194187166389509], [-0.204895409519981, 0.692104714174527, 0.692104714458033]]
     1313        """
     1314
     1315        u1 = self.variables[1]
     1316        u2 = self.variables[2]
     1317
     1318        solver = self._create_geodesic_ode_system()
     1319           
     1320        t_interval, n = tinterval[0:2], tinterval[2]
     1321        solver.y_0 = [p0[0], p0[1], v0[0], v0[1]]
     1322        solver.ode_solve(t_span=t_interval, num_points=n)
     1323
     1324        parsed_solution = \
     1325          [[vec[0], vec[1][0:2], vec[1][2:], self.eq_callable(vec[1][0], vec[1][1])] \
     1326           for vec in solver.solution]
     1327
     1328        return parsed_solution
     1329
     1330
     1331    @cached_method
     1332    def _create_pt_ode_system(self, curve, t):
     1333        """
     1334        Helper method to create a fast floating-point version of the parallel transport
     1335        equations, used by ``parallel_translation_numerical``.
     1336
     1337        INPUT:
     1338
     1339         - ``curve`` - curve in intrinsic coordinates along which to do parallel transport.
     1340         - ``t`` - curve parameter
     1341
     1342        EXAMPLES::
     1343           sage: p, q = var('p,q')
     1344           sage: sphere = ParametrizedSurface3D([cos(q)*cos(p),sin(q)*cos(p),sin(p)],[p,q],'sphere')
     1345           sage: s = var('s')
     1346           sage: ode = sphere._create_pt_ode_system((s, s), s)
     1347           sage: ode.function(0.0, (1.0, 1.0))
     1348           [-0.0, 0.0]
     1349                     
     1350        """
     1351
     1352        from sage.ext.fast_eval import fast_float
     1353        from sage.calculus.var import var
     1354        from sage.gsl.ode import ode_solver
     1355
     1356        u1 = self.variables[1]
     1357        u2 = self.variables[2]
     1358        v1, v2 = var('v1, v2')
     1359
     1360        du1 = diff(curve[0], t)
     1361        du2 = diff(curve[1], t)
     1362
     1363        C = self.connection_coefficients()
     1364        for coef in C:
     1365            C[coef] = C[coef].subs({u1: curve[0], u2: curve[1]})
     1366       
     1367        dv1 = - C[(1,1,1)]*v1*du1 - C[(1,2,1)]*(du1*v2 + du2*v1) - C[(2,2,1)]*du2*v2
     1368        dv2 = - C[(1,1,2)]*v1*du1 - C[(1,2,2)]*(du1*v2 + du2*v1) - C[(2,2,2)]*du2*v2
     1369        fun1 = fast_float(dv1, str(t), str(v1), str(v2))
     1370        fun2 = fast_float(dv2, str(t), str(v1), str(v2))
     1371
     1372        pt_ode = ode_solver()
     1373        pt_ode.function = lambda t, (v1, v2): [fun1(t, v1, v2), fun2(t, v1, v2)]
     1374        return pt_ode
     1375
     1376
     1377    def parallel_translation_numerical(self,curve,t,v0,tinterval):
     1378        r"""
     1379        Numerically solves the equations for parallel translation of a vector along a curve
     1380        on the surface.  Explicitly, the equations for parallel translation are given by
     1381
     1382         .. math::
     1383
     1384           \frac{d u^i}{dt} + u^j \frac{d c^k}{dt} \Gamma^i_{jk} = 0,
     1385
     1386        where $\Gamma^i_{jk}$ are the connection coefficients of the surface, the vector to
     1387        be transported has components $u^j$ and the curve along which to transport has
     1388        components $c^k$.
     1389
     1390        ALGORITHM:
     1391
     1392        The parallel transport equations are integrated forward in time using the
     1393        ode solvers from ``sage.gsl.ode``.  See ``_create_pt_ode_system``
     1394        for more details.
     1395
     1396        INPUT:
     1397
     1398         - ``curve`` - 2-tuple of functions which determine the curve with respect to
     1399           the local coordinate;
     1400
     1401         - ``t`` - symbolic variable denoting the curve parameter;
     1402
     1403         - ``v0`` - 2-tuple representing the initial vector;
     1404
     1405         - ``tinterval`` - list ``[a,b,N]``, where ``(a,b)`` is the domain of the curve
     1406           and ``N`` is the number of subdivision points.
     1407
     1408        OUTPUT:
     1409        The list consisting of lists ``[t, [v1(t),v2(t)]]``, where
     1410
     1411         - ``t`` is a subdivision point;
     1412
     1413         - ``[v1(t),v2(t)]`` is the list of coordinates of the vector parallel translated
     1414           along the curve.
     1415
     1416        EXAMPLES::
     1417
     1418           sage: p, q = var('p,q')
     1419           sage: v = [p,q]
     1420           sage: assume(cos(q)>0)
     1421           sage: sphere = ParametrizedSurface3D([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
     1422           sage: ss = var('ss')
     1423           sage: vv_array = sphere.parallel_translation_numerical([ss,ss],ss,[1.0,1.0],[0.0,pi/4,20])
     1424           sage: vv_array[0]
     1425           (0.000000000000000, [1.00000000000000, 1.00000000000000])
     1426           sage: vv_array[5]
     1427           (0.19634954084936207, [0.98060186565184926, 1.0389928973994338])
     1428
     1429        """
     1430
     1431        u1 = self.variables[1]
     1432        u2 = self.variables[2]
     1433
     1434        solver = self._create_pt_ode_system(tuple(curve), t)
     1435           
     1436        t_interval, n = tinterval[0:2], tinterval[2]
     1437        solver.y_0 = v0
     1438        solver.ode_solve(t_span=t_interval, num_points=n)
     1439
     1440        return solver.solution
     1441
     1442   
     1443
     1444
  • new file sage/riemann/vector_functions.py

    diff -r 120c07be6358 -r b834b692db4c 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 b834b692db4c 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',