Ticket #10132: trac_10132_final_allfiles_done_no_matrix_simplify.patch

File trac_10132_final_allfiles_done_no_matrix_simplify.patch, 60.5 KB (added by jvkersch, 9 years ago)
  • doc/en/reference/index.rst

    # HG changeset patch
    # User Joris Vankerschaver <joris.vankerschaver@gmail.com>
    # Date 1309210275 25200
    # Node ID 24d9f3cf95558ab757fae2ec5c7d9d9b06e5ed2a
    # Parent  f2d77baca5f74d1119bf657eb48524443a99cfbe
    #10132: Parametrization of (metric) surfaces in 3D Euclidean space
    
    diff -r f2d77baca5f7 -r 24d9f3cf9555 doc/en/reference/index.rst
    a b  
    9393   modabvar
    9494   modmisc
    9595   tensor
     96   riemannian_manifolds
    9697
    9798   history_and_license
    9899
  • new file doc/en/reference/riemannian_manifolds.rst

    diff -r f2d77baca5f7 -r 24d9f3cf9555 doc/en/reference/riemannian_manifolds.rst
    - +  
     1Differential Geometry of Curves and Surfaces
     2============================================
     3
     4.. toctree::
     5   :maxdepth: 2
     6
     7   sage/geometry/riemannian_manifolds/parametrized_surface3d
  • sage/all.py

    diff -r f2d77baca5f7 -r 24d9f3cf9555 sage/all.py
    a b  
    119119
    120120from sage.geometry.all   import *
    121121from sage.geometry.triangulation.all   import *
     122from sage.geometry.riemannian_manifolds.all   import *
    122123
    123124from sage.homology.all   import *
    124125
  • new file sage/geometry/riemannian_manifolds/__init__.py

    diff -r f2d77baca5f7 -r 24d9f3cf9555 sage/geometry/riemannian_manifolds/__init__.py
    - +  
     1# This comment is here so the file is non-empty (so Mercurial will check it in).
  • new file sage/geometry/riemannian_manifolds/all.py

    diff -r f2d77baca5f7 -r 24d9f3cf9555 sage/geometry/riemannian_manifolds/all.py
    - +  
     1from parametrized_surface3d import ParametrizedSurface3D
  • new file sage/geometry/riemannian_manifolds/parametrized_surface3d.py

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

    diff -r f2d77baca5f7 -r 24d9f3cf9555 setup.py
    a b  
    961961
    962962                     'sage.geometry',
    963963                     'sage.geometry.triangulation',
     964                     'sage.geometry.riemannian_manifolds',
    964965
    965966                     'sage.games',
    966967