Ticket #10132: trac_10132_revised_sage-5.13.beta2.patch

File trac_10132_revised_sage-5.13.beta2.patch, 82.4 KB (added by jvkersch, 7 years ago)
  • doc/en/reference/index.rst

    # HG changeset patch
    # User Joris Vankerschaver <joris.vankerschaver@gmail.com>
    # Date 1372872684 -3600
    #      Wed Jul 03 18:31:24 2013 +0100
    # Node ID e4a94a35818d833959d22aeeee04a49d85de6fe8
    # Parent  8377b5d138ca1db74ab7a674efc6f22a6542436f
    Trac #10132: Parametrization of (metric) surfaces in 3D Euclidean space
    
    diff --git a/doc/en/reference/index.rst b/doc/en/reference/index.rst
    a b  
    9191* :doc:`Combinatorial Geometry <geometry/index>`
    9292* :doc:`Cell Complexes and their Homology <homology/index>`
    9393* :doc:`Differential Forms <tensor/index>`
     94* :doc:`Parametrized Surfaces <riemannian_geometry/index>`
    9495
    9596Number Theory, Algebraic Geometry
    9697---------------------------------
  • new file doc/en/reference/riemannian_geometry/conf.py

    diff --git a/doc/en/reference/riemannian_geometry/conf.py b/doc/en/reference/riemannian_geometry/conf.py
    new file mode 100644
    - +  
     1# -*- coding: utf-8 -*-
     2#
     3# Sage documentation build configuration file, created by
     4# sphinx-quickstart on Thu Aug 21 20:15:55 2008.
     5#
     6# This file is execfile()d with the current directory set to its containing dir.
     7#
     8# The contents of this file are pickled, so don't put values in the namespace
     9# that aren't pickleable (module imports are okay, they're removed automatically).
     10#
     11# All configuration values have a default; values that are commented out
     12# serve to show the default.
     13
     14import sys, os
     15sys.path.append(os.environ['SAGE_DOC'])
     16from common.conf import *
     17
     18# settings for the intersphinx extension:
     19
     20ref_src = os.path.join(SAGE_DOC, 'en', 'reference')
     21ref_out = os.path.join(SAGE_DOC, 'output', 'html', 'en', 'reference')
     22intersphinx_mapping[ref_out] = None
     23
     24for doc in os.listdir(ref_src):
     25    if os.path.exists(os.path.join(ref_src, doc, 'index.rst')):
     26        intersphinx_mapping[os.path.join(ref_out, doc)] = None
     27
     28# We use the main document's title, if we can find it.
     29rst_file = open('index.rst', 'r')
     30rst_lines = rst_file.read().splitlines()
     31rst_file.close()
     32
     33title = u''
     34for i in xrange(len(rst_lines)):
     35    if rst_lines[i].startswith('==') and i > 0:
     36        title = rst_lines[i-1].strip()
     37        break
     38
     39# Otherwise, we use this directory's name.
     40name = os.path.basename(os.path.abspath('.'))
     41if not title:
     42    title = name.capitalize()
     43title = title.replace(u'`', u'$')
     44
     45# General information about the project.
     46project = u'Sage Reference Manual: ' + title
     47
     48# The name for this set of Sphinx documents.  If None, it defaults to
     49# "<project> v<release> documentation".
     50html_title = u'Sage Reference Manual v' + release + ': ' + title
     51
     52# A shorter title for the navigation bar.  Default is the same as html_title.
     53html_short_title = title
     54
     55# HTML theme (e.g., 'default', 'sphinxdoc').  The pages for the
     56# reference manual use a custom theme, a slight variant on the 'sage'
     57# theme, to set the links in the top line.
     58html_theme = 'sageref'
     59
     60# Output file base name for HTML help builder.
     61htmlhelp_basename = name
     62
     63# Grouping the document tree into LaTeX files. List of tuples (source
     64# start file, target name, title, author, document class
     65# [howto/manual]).
     66latex_documents = [
     67('index', name + '.tex', project, u'The Sage Development Team', 'manual')
     68]
     69
     70#Ignore all .rst in the _sage subdirectory
     71exclude_trees = exclude_trees + ['_sage']
     72
     73multidocs_is_master = False
  • new file doc/en/reference/riemannian_geometry/index.rst

    diff --git a/doc/en/reference/riemannian_geometry/index.rst b/doc/en/reference/riemannian_geometry/index.rst
    new file mode 100644
    - +  
     1Differential Geometry of Curves and Surfaces
     2============================================
     3
     4.. toctree::
     5   :maxdepth: 2
     6
     7   sage/geometry/riemannian_manifolds/parametrized_surface3d
     8   sage/geometry/riemannian_manifolds/surface3d_generators
     9
     10.. include:: ../footer.txt
  • sage/all.py

    diff --git a/sage/all.py b/sage/all.py
    a b  
    114114
    115115from sage.geometry.all   import *
    116116from sage.geometry.triangulation.all   import *
     117from sage.geometry.riemannian_manifolds.all   import *
    117118
    118119from sage.dynamics.all   import *
    119120
  • new file sage/geometry/riemannian_manifolds/__init__.py

    diff --git a/sage/geometry/riemannian_manifolds/__init__.py b/sage/geometry/riemannian_manifolds/__init__.py
    new file mode 100644
    - +  
     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 --git a/sage/geometry/riemannian_manifolds/all.py b/sage/geometry/riemannian_manifolds/all.py
    new file mode 100644
    - +  
     1from sage.misc.lazy_import import lazy_import
     2lazy_import('sage.geometry.riemannian_manifolds.parametrized_surface3d',
     3            'ParametrizedSurface3D')
     4lazy_import('sage.geometry.riemannian_manifolds.surface3d_generators',
     5            'surfaces')
     6
  • new file sage/geometry/riemannian_manifolds/parametrized_surface3d.py

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

    diff --git a/sage/geometry/riemannian_manifolds/surface3d_generators.py b/sage/geometry/riemannian_manifolds/surface3d_generators.py
    new file mode 100644
    - +  
     1r"""
     2Common parametrized surfaces in 3D.
     3
     4AUTHORS::
     5
     6- Joris Vankerschaver (2012-06-16)
     7
     8"""
     9
     10#*****************************************************************************
     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.symbolic.constants import pi
     19from sage.functions.log import log
     20from sage.functions.trig import sin, cos, tan
     21from sage.functions.hyperbolic import cosh, sinh, tanh
     22from sage.symbolic.ring import SR, var
     23from sage.geometry.riemannian_manifolds.parametrized_surface3d import \
     24    ParametrizedSurface3D
     25
     26
     27class SurfaceGenerators():
     28    """
     29    A class consisting of generators for several common parametrized surfaces
     30    in 3D.
     31
     32    """
     33    @staticmethod
     34    def Catenoid(c=1, name="Catenoid"):
     35        r"""
     36        Returns a catenoid surface, with parametric representation
     37
     38        .. MATH::
     39
     40            \begin{aligned}
     41              x(u, v) & = c \cosh(v/c) \cos(u); \\
     42              y(u, v) & = c \cosh(v/c) \sin(u); \\
     43              z(u, v) & = v.
     44            \end{aligned}
     45
     46        INPUT:
     47
     48        - ``c`` -- surface parameter.
     49
     50        - ``name`` -- string. Name of the surface.
     51
     52
     53        EXAMPLES::
     54
     55            sage: cat = surfaces.Catenoid(); cat
     56            Parametrized surface ('Catenoid') with equation (cos(u)*cosh(v), cosh(v)*sin(u), v)
     57            sage: cat.plot()
     58
     59        """
     60        u, v = var('u, v')
     61        catenoid_eq = [c*cosh(v/c)*cos(u), c*cosh(v/c)*sin(u), v]
     62        coords = ((u, 0, 2*pi), (v, -1, 1))
     63
     64        return ParametrizedSurface3D(catenoid_eq, coords, name)
     65
     66    @staticmethod
     67    def Crosscap(r=1, name="Crosscap"):
     68        r"""
     69        Returns a crosscap surface, with parametrization
     70
     71        .. MATH::
     72
     73            \begin{aligned}
     74              x(u, v) & = r(1 + \cos(v)) \cos(u); \\
     75              y(u, v) & = r(1 + \cos(v)) \sin(u); \\
     76              z(u, v) & = - r\tanh(u - \pi) \sin(v).
     77            \end{aligned}
     78
     79        INPUT:
     80
     81        - ``r`` -- surface parameter.
     82
     83        - ``name`` -- string. Name of the surface.
     84
     85        EXAMPLES::
     86
     87            sage: crosscap = surfaces.Crosscap(); crosscap
     88            Parametrized surface ('Crosscap') with equation ((cos(v) + 1)*cos(u), (cos(v) + 1)*sin(u), -sin(v)*tanh(-pi + u))
     89            sage: crosscap.plot()
     90
     91        """
     92
     93        u, v = var('u, v')
     94        crosscap_eq = [r*(1+cos(v))*cos(u), r*(1+cos(v))*sin(u),
     95                       -tanh(u-pi)*r*sin(v)]
     96        coords = ((u, 0, 2*pi), (v, 0, 2*pi))
     97
     98        return ParametrizedSurface3D(crosscap_eq, coords, name)
     99
     100    @staticmethod
     101    def Dini(a=1, b=1, name="Dini's surface"):
     102        r"""
     103        Returns Dini's surface, with parametrization
     104
     105        .. MATH::
     106
     107            \begin{aligned}
     108              x(u, v) & = a \cos(u)\sin(v); \\
     109              y(u, v) & = a \sin(u)\sin(v); \\
     110              z(u, v) & = u + \log(\tan(v/2)) + \cos(v).
     111            \end{aligned}
     112
     113        INPUT:
     114
     115        - ``a, b`` -- surface parameters.
     116
     117        - ``name`` -- string. Name of the surface.
     118
     119        EXAMPLES::
     120
     121            sage: dini = surfaces.Dini(a=3, b=4); dini
     122            Parametrized surface ('Dini's surface') with equation (3*cos(u)*sin(v), 3*sin(u)*sin(v), 4*u + 3*cos(v) + 3*log(tan(1/2*v)))
     123            sage: dini.plot()  # not tested -- known bug (see #10132)
     124
     125        """
     126
     127        u, v = var('u, v')
     128        dini_eq = [a*cos(u)*sin(v), a*sin(u)*sin(v),
     129                   a*(cos(v) + log(tan(v/2))) + b*u]
     130        coords = ((u, 0, 2*pi), (v, 0, 2*pi))
     131
     132
     133        return ParametrizedSurface3D(dini_eq, coords, name)
     134
     135    @staticmethod
     136    def Ellipsoid(center=(0,0,0), axes=(1,1,1), name="Ellipsoid"):
     137        r"""
     138        Returns an ellipsoid centered at ``center`` whose semi-principal axes
     139        have lengths given by the components of ``axes``. The
     140        parametrization of the ellipsoid is given by
     141
     142        .. MATH::
     143
     144            \begin{aligned}
     145              x(u, v) & = x_0 + a \cos(u) \cos(v); \\
     146              y(u, v) & = y_0 + b \sin(u) \cos(v); \\
     147              z(u, v) & = z_0 + c \sin(v).
     148            \end{aligned}
     149
     150        INPUT:
     151
     152        - ``center`` -- 3-tuple. Coordinates of the center of the ellipsoid.
     153
     154        - ``axes`` -- 3-tuple. Lengths of the semi-principal axes.
     155
     156        - ``name`` -- string. Name of the ellipsoid.
     157
     158        EXAMPLES::
     159
     160            sage: ell = surfaces.Ellipsoid(axes=(1, 2, 3)); ell
     161            Parametrized surface ('Ellipsoid') with equation (cos(u)*cos(v), 2*cos(v)*sin(u), 3*sin(v))
     162            sage: ell.plot()
     163
     164        """
     165
     166        u, v = var ('u, v')
     167        x, y, z = center
     168        a, b, c = axes
     169        ellipsoid_parametric_eq = [x + a*cos(u)*cos(v),
     170                                   y + b*sin(u)*cos(v),
     171                                   z + c*sin(v)]
     172        coords = ((u, 0, 2*pi), (v, -pi/2, pi/2))
     173
     174        return ParametrizedSurface3D(ellipsoid_parametric_eq, coords, name)
     175
     176    @staticmethod
     177    def Enneper(name="Enneper's surface"):
     178        r"""
     179        Returns Enneper's surface, with parametrization
     180
     181        .. MATH::
     182
     183            \begin{aligned}
     184              x(u, v) & = u(1 - u^2/3 + v^2)/3; \\
     185              y(u, v) & = -v(1 - v^2/3 + u^2)/3; \\
     186              z(u, v) & = (u^2 - v^2)/3.
     187            \end{aligned}
     188
     189        INPUT:
     190
     191        - ``name`` -- string. Name of the surface.
     192
     193        EXAMPLES::
     194
     195            sage: enn = surfaces.Enneper(); enn
     196            Parametrized surface ('Enneper's surface') with equation (-1/9*(u^2 - 3*v^2 - 3)*u, -1/9*(3*u^2 - v^2 + 3)*v, 1/3*u^2 - 1/3*v^2)
     197            sage: enn.plot()
     198
     199        """
     200
     201        u, v = var('u, v')
     202        enneper_eq = [u*(1-u**2/3+v**2)/3, -v*(1-v**2/3+u**2)/3, (u**2-v**2)/3]
     203        coords = ((u, -3, 3), (v, -3, 3))
     204
     205        return ParametrizedSurface3D(enneper_eq, coords, name)
     206
     207    @staticmethod
     208    def Helicoid(h=1, name="Helicoid"):
     209        r"""
     210        Returns a helicoid surface, with parametrization
     211
     212        .. MATH::
     213
     214            \begin{aligned}
     215              x(\rho, \theta) & = \rho \cos(\theta); \\
     216              y(\rho, \theta) & = \rho \sin(\theta); \\
     217              z(\rho, \theta) & = h\theta/(2\pi).
     218            \end{aligned}
     219
     220        INPUT:
     221
     222        - ``h`` -- distance along the z-axis between two
     223          successive turns of the helicoid.
     224
     225        - ``name`` -- string. Name of the surface.
     226
     227        EXAMPLES::
     228
     229            sage: helicoid = surfaces.Helicoid(h=2); helicoid
     230            Parametrized surface ('Helicoid') with equation (rho*cos(theta), rho*sin(theta), theta/pi)
     231            sage: helicoid.plot()
     232
     233        """
     234
     235        rho, theta = var('rho, theta')
     236        helicoid_eq = [rho*cos(theta), rho*sin(theta), h*theta/(2*pi)]
     237        coords = ((rho, -2, 2), (theta, 0, 2*pi))
     238
     239        return ParametrizedSurface3D(helicoid_eq, [rho, theta], name)
     240
     241    @staticmethod
     242    def Klein(r=1, name="Klein bottle"):
     243        r"""
     244        Returns the Klein bottle, in the figure-8 parametrization given by
     245
     246        .. MATH::
     247
     248            \begin{aligned}
     249              x(u, v) & = (r + \cos(u/2)\cos(v) - \sin(u/2)\sin(2v)) \cos(u); \\
     250              y(u, v) & = (r + \cos(u/2)\cos(v) - \sin(u/2)\sin(2v)) \sin(u); \\
     251              z(u, v) & = \sin(u/2)\cos(v) + \cos(u/2)\sin(2v).
     252            \end{aligned}
     253
     254        INPUT:
     255
     256        - ``r`` -- radius of the "figure-8" circle.
     257
     258        - ``name`` -- string. Name of the surface.
     259
     260        EXAMPLES::
     261
     262            sage: klein = surfaces.Klein(); klein
     263            Parametrized surface ('Klein bottle') with equation (-(sin(1/2*u)*sin(2*v) - cos(1/2*u)*sin(v) - 1)*cos(u), -(sin(1/2*u)*sin(2*v) - cos(1/2*u)*sin(v) - 1)*sin(u), cos(1/2*u)*sin(2*v) + sin(1/2*u)*sin(v))
     264            sage: klein.plot()
     265
     266        """
     267
     268        u, v = var('u, v')
     269        x = (r + cos(u/2)*sin(v) - sin(u/2)*sin(2*v))*cos(u)
     270        y = (r + cos(u/2)*sin(v) - sin(u/2)*sin(2*v))*sin(u)
     271        z = sin(u/2)*sin(v) + cos(u/2)*sin(2*v)
     272        klein_eq = [x, y, z]
     273        coords = ((u, 0, 2*pi), (v, 0, 2*pi))
     274
     275        return ParametrizedSurface3D(klein_eq, coords, name)
     276
     277    @staticmethod
     278    def MonkeySaddle(name="Monkey saddle"):
     279        r"""
     280        Returns a monkey saddle surface, with equation
     281
     282        .. MATH::
     283
     284            z = x^3 - 3xy^2.
     285
     286        INPUT:
     287
     288        - ``name`` -- string. Name of the surface.
     289
     290        EXAMPLES::
     291
     292            sage: saddle = surfaces.MonkeySaddle(); saddle
     293            Parametrized surface ('Monkey saddle') with equation (u, v, u^3 - 3*u*v^2)
     294            sage: saddle.plot()
     295
     296        """
     297
     298        u, v = var('u, v')
     299        monkey_eq = [u, v, u**3 - 3*u*v**2]
     300        coords = ((u, -2, 2), (v, -2, 2))
     301
     302        return ParametrizedSurface3D(monkey_eq, coords, name)
     303
     304    @staticmethod
     305    def Paraboloid(a=1, b=1, c=1, elliptic=True, name=None):
     306        r"""
     307        Returns a paraboloid with equation
     308
     309        .. MATH::
     310
     311            \frac{z}{c} = \pm \frac{x^2}{a^2} + \frac{y^2}{b^2}
     312
     313        When the plus sign is selected, the paraboloid is elliptic. Otherwise
     314        the surface is a hyperbolic paraboloid.
     315
     316        INPUT:
     317
     318        - ``a``, ``b``, ``c`` -- Surface parameters.
     319
     320        - ``elliptic`` (default: True) -- whether to create an elliptic or
     321          hyperbolic paraboloid.
     322
     323        - ``name`` -- string. Name of the surface.
     324
     325        EXAMPLES::
     326
     327            sage: epar = surfaces.Paraboloid(1, 3, 2); epar
     328            Parametrized surface ('Elliptic paraboloid') with equation (u, v, 2*u^2 + 2/9*v^2)
     329            sage: epar.plot()
     330
     331            sage: hpar = surfaces.Paraboloid(2, 3, 1, elliptic=False); hpar
     332            Parametrized surface ('Hyperbolic paraboloid') with equation (u, v, -1/4*u^2 + 1/9*v^2)
     333            sage: hpar.plot()
     334
     335        """
     336
     337        u, v = var('u, v')
     338        x = u; y = v
     339        if elliptic:
     340            z = c*(v**2/b**2 + u**2/a**2)
     341        else:
     342            z = c*(v**2/b**2 - u**2/a**2)
     343        paraboloid_eq = [x, y, z]
     344        coords = ((u, -3, 3), (v, -3, 3))
     345
     346        if name is None:
     347            if elliptic:
     348                name = "Elliptic paraboloid"
     349            else:
     350                name = "Hyperbolic paraboloid"
     351
     352        return ParametrizedSurface3D(paraboloid_eq, coords, name)
     353
     354    @staticmethod
     355    def Sphere(center=(0,0,0), R=1, name="Sphere"):
     356        r"""
     357        Returns a sphere of radius ``R`` centered at ``center``.
     358
     359        INPUT:
     360
     361        - ``center`` -- 3-tuple, center of the sphere.
     362
     363        - ``R`` -- Radius of the sphere.
     364
     365        - ``name`` -- string. Name of the surface.
     366
     367        EXAMPLES::
     368
     369            sage: sphere = surfaces.Sphere(center=(0, 1, -1), R=2); sphere
     370            Parametrized surface ('Sphere') with equation (2*cos(u)*cos(v), 2*cos(v)*sin(u) + 1, 2*sin(v) - 1)
     371            sage: sphere.plot()
     372
     373        Note that the radius of the sphere can be negative. The surface thus
     374        obtained is equal to the sphere (or part thereof) with positive radius,
     375        whose coordinate functions have been multiplied by -1. Compare for
     376        instant the first octant of the unit sphere with positive radius::
     377
     378            sage: octant1 = surfaces.Sphere(R=1); octant1
     379            Parametrized surface ('Sphere') with equation (cos(u)*cos(v), cos(v)*sin(u), sin(v))
     380            sage: octant1.plot((0, pi/2), (0, pi/2))
     381
     382        with the first octant of the unit sphere with negative radius::
     383
     384            sage: octant2 = surfaces.Sphere(R=-1); octant2
     385            Parametrized surface ('Sphere') with equation (-cos(u)*cos(v), -cos(v)*sin(u), -sin(v))
     386            sage: octant2.plot((0, pi/2), (0, pi/2))
     387
     388        """
     389
     390        return SurfaceGenerators.Ellipsoid(center, (R, R, R), name)
     391
     392    @staticmethod
     393    def Torus(r=2, R=3, name="Torus"):
     394        r"""
     395        Returns a torus obtained by revolving a circle of radius ``r`` around
     396        a coplanar axis ``R`` units away from the center of the circle. The
     397        parametrization used is
     398
     399        .. MATH::
     400
     401            \begin{aligned}
     402              x(u, v) & = (R + r \cos(v)) \cos(u); \\
     403              y(u, v) & = (R + r \cos(v)) \sin(u); \\
     404              z(u, v) & = r \sin(v).
     405            \end{aligned}
     406
     407        INPUT:
     408
     409        - ``r``, ``R`` -- Minor and major radius of the torus.
     410
     411        - ``name`` -- string. Name of the surface.
     412
     413        EXAMPLES::
     414
     415            sage: torus = surfaces.Torus(); torus
     416            Parametrized surface ('Torus') with equation ((2*cos(v) + 3)*cos(u), (2*cos(v) + 3)*sin(u), 2*sin(v))
     417            sage: torus.plot()
     418
     419        """
     420
     421        u, v = var('u, v')
     422        torus_eq = [(R+r*cos(v))*cos(u), (R+r*cos(v))*sin(u), r*sin(v)]
     423        coords = ((u, 0, 2*pi), (v, 0, 2*pi))
     424
     425        return ParametrizedSurface3D(torus_eq, coords, name)
     426
     427    @staticmethod
     428    def WhitneyUmbrella(name="Whitney's umbrella"):
     429        r"""
     430        Returns Whitney's umbrella, with parametric representation
     431
     432        .. MATH::
     433
     434            x(u, v) = uv, \quad y(u, v) = u, \quad z(u, v) = v^2.
     435
     436        INPUT:
     437
     438        - ``name`` -- string. Name of the surface.
     439
     440        EXAMPLES::
     441
     442            sage: whitney = surfaces.WhitneyUmbrella(); whitney
     443            Parametrized surface ('Whitney's umbrella') with equation (u*v, u, v^2)
     444            sage: whitney.plot()
     445
     446        """
     447
     448        u, v = var('u, v')
     449        whitney_eq = [u*v, u, v**2]
     450        coords = ((u, -1, 1), (v, -1, 1))
     451
     452        return ParametrizedSurface3D(whitney_eq, coords, name)
     453
     454
     455# Easy access to the surface generators
     456surfaces = SurfaceGenerators()
     457