Ticket #10132: trac_10123_final_allfiles_done.patch

File trac_10123_final_allfiles_done.patch, 58.7 KB (added by jvkersch, 10 years ago)
  • doc/en/reference/index.rst

    # HG changeset patch
    # User Joris Vankerschaver <joris.vankerschaver@gmail.com>
    # Date 1309210275 25200
    # Node ID 2c8cabc1aa8f2c59d1dd952e80fc435571637640
    # Parent  841ede959e174beb02127f6ca660efd491480cb9
    #10132: Parametrization of (metric) surfaces in 3D Euclidean space
    
    diff -r 841ede959e17 -r 2c8cabc1aa8f doc/en/reference/index.rst
    a b  
    9292   modabvar
    9393   modmisc
    9494   tensor
     95   riemannian_manifolds
    9596
    9697   history_and_license
    9798
  • new file doc/en/reference/riemannian_manifolds.rst

    diff -r 841ede959e17 -r 2c8cabc1aa8f 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 841ede959e17 -r 2c8cabc1aa8f 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 841ede959e17 -r 2c8cabc1aa8f 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 841ede959e17 -r 2c8cabc1aa8f sage/geometry/riemannian_manifolds/all.py
    - +  
     1from parametrized_surface3d import ParametrizedSurface3D
  • new file sage/geometry/riemannian_manifolds/parametrized_surface3d.py

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

    diff -r 841ede959e17 -r 2c8cabc1aa8f setup.py
    a b  
    961961
    962962                     'sage.geometry',
    963963                     'sage.geometry.triangulation',
     964                     'sage.geometry.riemannian_manifolds',
    964965
    965966                     'sage.games',
    966967