Ticket #10132: trac_10132_final_all_files_done.patch

File trac_10132_final_all_files_done.patch, 61.5 KB (added by mikarm, 10 years ago)
  • doc/en/reference/index.rst

    # HG changeset patch
    # User Joris Vankerschaver <joris.vankerschaver@gmail.com>
    # Date 1299994529 28800
    # Node ID ddfbc751857d1981b3deeddb3c5569510d999542
    # Parent  120c07be6358d93bcff503363d379c26b8342f2b
    [mq]: 10132allfiles
    
    diff -r 120c07be6358 -r ddfbc751857d 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 ddfbc751857d 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 ddfbc751857d 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 ddfbc751857d 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 ddfbc751857d sage/riemann/all.py
    - +  
     1from parametrized_surface3d import ParametrizedSurface3D
  • new file sage/riemann/parametrized_surface3d.py

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

    diff -r 120c07be6358 -r ddfbc751857d sage/riemann/vector_functions.py
    - +  
     1"""
     2Utility functions to simplify vectors and matrices over the symbolic ring.
     3
     4AUTHORS:
     5        - Mikhail Malakhaltsev (2010-09-25): initial version
     6        - Joris Vankerschaver (2011-03-06): docstrings
     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
     16from sage.modules.free_module_element import vector
     17from sage.matrix.constructor import matrix
     18
     19
     20def simplify_vector(a):
     21    """
     22    Simplify a symbolic vector.
     23
     24    EXAMPLES::
     25
     26        sage: from sage.riemann.vector_functions import simplify_vector
     27        sage: x, y = var('x, y')
     28        sage: v = vector([cos(x)**2 + sin(x)**2, 0]); v
     29        (sin(x)^2 + cos(x)^2, 0)
     30        sage: simplify_vector(v)
     31        (1, 0)
     32
     33    """
     34    return(vector([i.simplify_full() for i in a]))
     35
     36
     37def simplify_matrix(A):
     38    """
     39    Simplify a symbolic matrix.
     40
     41    EXAMPLES::
     42
     43        sage: from sage.riemann.vector_functions import simplify_matrix
     44        sage: x, y = var('x, y')
     45        sage: A = matrix([[1, x], [x, 1]])
     46        sage: B = A*A.inverse(); B
     47        [                                   1                                    0]
     48        [-(x^2/(x^2 - 1) - 1)*x + x/(x^2 - 1)          x^2/(x^2 - 1) - 1/(x^2 - 1)]
     49        sage: simplify_matrix(B)
     50        [1 0]
     51        [0 1]
     52
     53    """
     54
     55    return matrix([[ xx.simplify_full() for xx in row ] for row in A])
     56
     57
  • setup.py

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