# Ticket #10132: parametrized_surface3d.2.py

File parametrized_surface3d.2.py, 40.1 KB (added by mikarm, 11 years ago)

parametrized_surface3d.py with small changes in parallel_translation_numeric_new

Line
1"""
2Basic invariants of parametrized surfaces.
3
4
5AUTHORS:
6        - Mikhail Malakhaltsev (2010-09-25): initial version
7
8"""
9#*****************************************************************************
10#       Copyright (C) 2010  Mikhail Malakhaltsev <mikarm@gmail.com>
11#
14#*****************************************************************************
15
16
17# The following is maybe too specific ...
18#from vector_functions import *
19
20# mikarm: This is for the moment, at the final version will be replaced by import
21attach('/home/prince/MY_SAGE_MODULES/vector_functions.py')
22
23from sage.structure.sage_object import SageObject
24from sage.calculus.functional import diff
25from sage.functions.other import sqrt
26
27from sage.misc.cachefunc import cached_method
28
29
30
31# TODO (distant future): update class...
32class GenericManifold(SageObject):
33    pass
34
35
36class ParametrizedSurface3D(GenericManifold):
37    """
38    This is a class for finding basic invariants of surfaces.
39
40    USAGE:
41    parametrized_surface3d(surface_equation,variables,name_of_surface)
42    where
43      surface_equation is a vector or list of three components
44      variables is a list of two variables
45      name_of_surface is a string (optional)
46
47    Warning: The orientation on the surface is given by the parametrization, that is the positive frame is
48    $\partial_1 \\vec r$, $\partial_2 \\vec r$.
49
50    This class includes the following methods:
51      natural_frame
52      normal_vector
53      first_fundamental_form_coefficients
54      first_fundamental_form
55      first_fundamental_form_inverse_coefficients
56      first_fundamental_form_inverse_coefficients
57      area_form_squared
58      area_form
59      rotation
60      orthonormal_frame
61      lie_bracket
62      frame_structure_functions
63      second_order_natural_frame
64      second_fundamental_form_coefficients
65      second_fundamental_form
66      gauss_curvature
67      mean_curvature
68      principal curvatures
69      principal directions
70      connection_coefficients
71      geodesics_numerical
72      parallel_translation_numerical
73
74
75    EXAMPLES::
76
77      sage: var('u,v')
78      sage: paraboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'paraboloid')
79      sage: paraboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v])
80
81    """
82
83    def __init__(self,equation,variables,*name):
84
85        self.equation = equation
86        self.variables_list = variables
87        self.variables = {1:self.variables_list[0],2:self.variables_list[1]}
88        self.name = name
89
90        # Callable version of the underlying equation
91        def eq_callable(u, v):
92            u1, u2 = self.variables_list
93            return [fun.subs({u1: u, u2: v}) for fun in self.equation]
94
95        self.eq_callable = eq_callable
96
97        # Various index tuples
98        self.index_list = [(i,j) for i in (1,2) for j in (1,2)]
99        self.index_list_3=[(i,j,k) for i in (1,2) for j in (1,2) for k in (1,2)]
100
101
102    @cached_method
103    def natural_frame(self,vector_number=0):
104        """
105        This functions find the natural frame of the parametrized surface
106        INPUT:
107        the argument can be empty, or equal to 1, or to 2
108
109        OUTPUT:
110        if argument is equal to 1 or 2, the output is the corresponding vector of the natural frame
111        if argument is empty, the output is the dictionary of the vector fields of the natural frame
112
113        EXAMPLES::
114
115            sage:  eparaboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
116            sage:  eparaboloid.natural_frame()
117            {1: (1, 0, 2*u), 2: (0, 1, 2*v)}
118            sage:  eparaboloid.natural_frame(1)
119            (1, 0, 2*u)
120            sage:  eparaboloid.natural_frame(2)
121            (0, 1, 2*v)
122        """
123
124        dr1 = vector([diff(f,self.variables[1]).simplify_full() for f in self.equation])
125        dr2 = vector([diff(f,self.variables[2]).simplify_full() for f in self.equation])
126
127        return {1:dr1, 2:dr2}
128
129
130    @cached_method
131    def normal_vector(self, normalized=False):
132        """
133        This functions find the normal vector of the parametrized surface
134
135        INPUT:
136        empty argument, of a nonzero real number l
137
138        OUTPUT:
139        without arguments gives the normal vector which is the vector product of the
140        natural frame vectors;
141
142        with the argument l gives the normal vector of length l (the orientation determined
143        by the sign of l)
144
145        EXAMPLES::
146
147          sage: eparaboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
148          sage: eparaboloid.normal_vector()
149          (-2*u,-2*v,1)
150          sage: eparaboloid.normal_vector(1)
151          (-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))
152
153        """
154
155        dr = self.natural_frame()
156        normal = dr[1].cross_product(dr[2])
157
158        if normalized:
159            normal /= normal.norm()
160        return simplify_vector(normal)
161
162
163
164    # The following can become part of the docstring at some point...
165
166    # I've defined a cached (private) method which does not check its arguments and
167    # computes one component at a time.  Note that caching is just a matter of
168    # prepending "@cached_method" to the method signature.  For this to work, the
169    # arguments of the method must be hashable (e.g tuple, string, ...) and the method
170    # must not change the internal state of the object.  For all of the methods in this
171    # class, this is the case.
172    #
173    # I then added a method that does simple argument checking, converting the arguments to
174    # tuples, and invokes the private method to do the work.  Note that this method should
175    # not be cached, as the user might specify the index of the component as a list, which is
176    # not hashable.
177    #
178    # A second method returns a dictionary of all the components, again invoking the hidden
179    # helper method.
180
181    @cached_method
182    def _compute_first_fundamental_form_coefficient(self, index):
183        dr = self.natural_frame()
184        return (dr[index[0]]*dr[index[1]]).simplify_full()
185
186    def first_fundamental_form_coefficient(self, index):
187        index = tuple(sorted(index))
188        if index not in self.index_list:
189            raise ValueError, "Index %s out of bounds." % str(index)
190        return self._compute_first_fundamental_form_coefficient(index)
191
192    @cached_method
193    def first_fundamental_form_coefficients(self):
194
195        """
196        This function gives the coefficients of the first fundamental form
197
198        INPUT:
199        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
200
201        OUTPUT:
202        with empty argument the output is the dictionary of coefficients of the first fundamental form
203
204        with given indices the output is the corresponding coefficient of the first fundamental form,
205        for example index1 = 1, index2 = 2 gives the coefficient $g_{12}$
206
207        EXAMPLES::
208
209            sage: var('u,v')
210            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
211            sage: sphere.first_fundamental_form_coefficients()
212            {(1, 2): 0, (1, 1): cos(v)^2, (2, 1): 0, (2, 2): 1}
213
214            sage: sphere.first_fundamental_form_coefficients(1,1)
215            cos(v)^2
216
217            sage: sphere.first_fundamental_form_coefficients(2,1)
218            0
219
220            sage: sphere.first_fundamental_form_coefficients(3,2)
221            'The argument is not appropriate. Read doc'
222        """
223        coefficients = {}
224        for index in self.index_list:
225            sorted_index = list(sorted(index))
226            coefficients[index] = \
227                self._compute_first_fundamental_form_coefficient(index)
228        return coefficients
229
230
231
232    def first_fundamental_form(self,vector1,vector2):
233        """
234        Finds the value of first fundamental form on two vectors
235
236        INPUT:
237        Two vectors $v=(v^1,v^2)$ and $w=(w^1,w^2)$
238
239        OUTPUT:
240        $g_{11} v^1 w^1 + g_{12}(v^1 w^2 + v^2 w^1) + g_{22} v^2 w^2$
241
242        #EXAMPLES::
243
244            sage: var('u,v')
245            sage: var ('v1,v2,w1,w2')
246            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
247
248            sage: sphere.first_fundamental_form(vector([v1,v2]),vector([w1,w2]))
249            v1*w1*cos(v)^2 + v2*w2
250
251            sage: vv = vector([1,2])
252            sage: sphere.first_fundamental_form(vv,vv)
253            cos(v)^2 + 4
254
255            sage: sphere.first_fundamental_form([1,1],[2,1])
256            2*cos(v)^2 + 1
257        """
258        gg = self.first_fundamental_form_coefficients()
259        return sum(gg[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list)
260
261
262    @cached_method
263    def area_form_squared(self):
264        """
265        Gives $g_{11}g_{22} - g_{12}^2$, which is the square of the coefficient of the area form
266
267        INPUT:
268        No arguments
269
270        OUTPUT:
271        $g_{11}g_{22} - g_{12}^2$
272
273        EXAMPLES::
274
275            sage: var('u,v')
276            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
277            sage: sphere.area_form_squared()
278            cos(v)^2
279
280        """
281        gg = self.first_fundamental_form_coefficients()
282        return  (gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2).simplify_full()
283
284
285    @cached_method
286    def area_form(self):
287        """
288        Gives $\sqrt{g_{11}g_{22} - g_{12}^2}$, which is the coefficient of the area form
289
290        INPUT:
291        No arguments
292
293        OUTPUT:
294        $\sqrt{g_{11}g_{22} - g_{12}^2}$
295
296        EXAMPLES::
297
298            sage: var('u,v')
299            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
300            sage: sphere.area_form_squared()
301            abs(cos(v))
302
303        """
304        return sqrt(self.area_form_squared()).simplify_full()
305
306
307    @cached_method
308    def first_fundamental_form_inverse_coefficients(self):
309        """
310        Gives $g^{ij}$
311
312        INPUT:
313        No arguments
314
315        OUTPUT:
316        dictionary {(1,1):$g^{11}$, (1,2):$g^{12}$, (2,1):$g^{21}$, (2,2):$g^{22}$
317
318        EXAMPLES::
319
320            sage: var('u,v')
321            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
322            sage: sphere.first_fundamental_form_inverse_coefficients()
323            {(1, 2): 0, (1, 1): cos(v)^(-2), (2, 1): 0, (2, 2): 1}
324        """
325
326        gg = self.first_fundamental_form_coefficients()
327        DD = gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2
328
329        gi11 = (gg[(2,2)]/DD).simplify_full()
330        gi12 = (-gg[(1,2)]/DD).simplify_full()
331        gi21 = gi12
332        gi22 = (gg[(1,1)]/DD).simplify_full()
333
334        return {(1,1):gi11,(1,2):gi12,(2,1):gi21,(2,2):gi22}
335
336    def first_fundamental_form_inverse_coefficient(self, index):
337        index = tuple(sorted(index))
338        if index not in self.index_list:
339            raise ValueError, "Index %s out of bounds." % str(index)
340        return self._compute_first_fundamental_form_inverse_coefficients()[index]
341
342
343    @cached_method
344    def rotation(self,theta):
345        """
346        Gives the matrix of the operator of rotation on the given angle $\\theta$ with respect to the natural frame
347
348        INPUT:
349        given angle $\\theta$
350
351        OUTPUT:
352        matrix of the operator of rotation on $\\theta$ with respect to the natural frame
353
354        ALGORITHM:
355        The operator of rotation on $\pi/2$ is $J^i_j = g^{ik}\omega_{jk}$, where $\omega$ is the area form
356        The operator of rotation on angle $\\theta$ is $\cos(\\theta) I + sin(\\theta) J$
357
358        EXAMPLES::
359
360            sage: var('u,v')
361            sage: assume(cos(v)>0)
362            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
363            sage: rotation = sphere.rotation(pi/3)
364            sage: rotation^3
365            [-1  0]
366            [ 0 -1]
367            # it is true because the rotation at $\pi/3$ applied three times gives rotation a $\pi$
368        """
369
370        from sage.functions.trig import sin, cos
371
372        gi = self.first_fundamental_form_inverse_coefficients()
373        w12 = self.area_form()
374        RR11 = (cos(theta) + sin(theta)*gi[1,2]*w12).simplify_full()
375        RR12 = (- sin(theta)*gi[1,1]*w12).simplify_full()
376        RR21 = (sin(theta)*gi[2,2]*w12).simplify_full()
377        RR22 = (cos(theta) - sin(theta)*gi[2,1]*w12).simplify_full()
378        rotation = matrix([[RR11,RR12],[RR21,RR22]])
379        return rotation
380
381
382    # XXXXXX
383
384
385
386
387
388    @cached_method
389    def orthonormal_frame_vector(self, vector_number=0, coordinates='ext'):
390        """
391        Gives the orthonormal frame field of the surface
392
393        INPUT:
394        (vector_number,coordinates), where
395        vector number is 0 (default), 1 or 2
396        coordinates is 'ext'(default) or 'int'
397
398
399        OUTPUT:
400        If vector_number is equal to 1 or 2, the output is the corresponding vector of an orthonormal frame.
401        If argument vector_number is empty, the output is the dictionary of the vector fields of an orthonormal frame.
402
403        If coordinates='ext', output is coordinates in $\mathbb{R}^3$
404        If coordinates='int', output is coordinates with respect to the natural frame
405
406        ALGORITHM:
407        We normalize the first vector $\\vec e_1$ of the natural frame and then get the second frame vector
408        $\\vec e_2 = [\\vec n, \\vec e_1]$.
409
410        EXAMPLES::
411
412            sage: var('u,v')
413            sage: assume(cos(v)>0)
414            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
415            sage: EE = sphere.orthonormal_frame()
416            sage: (EE[1]*EE[1]).simplify_full()
417            1
418            sage: (EE[1]*EE[2]).simplify_full()
419            0
420            sage: simplify_vector(vector_product(EE[1],EE[2])-sphere.normal_vector(1))
421            [0,0,0]
422
423            sage: EE_int = sphere.orthonormal_frame(coordinates='int')
424            sage: EE_int
425            {1: (1/cos(v), 0), 2: (0, 1)}
426            sage: sphere.first_fundamental_form(EE_int[1],EE_int[1])
427            1
428            sage: sphere.first_fundamental_form(EE_int[1],EE_int[2])
429            0
430            sage: sphere.first_fundamental_form(EE_int[2],EE_int[2])
431            1
432
433        """
434
435        if vector_number not in [1, 2]:
436            raise ValueError, "Basis vector number out of range."
437        return self.orthonormal_frame(coordinates)[vector_number]
438
439
440    @cached_method
441    def orthonormal_frame(self, coordinates='ext'):
442
443        from sage.symbolic.constants import pi
444
445        if coordinates not in ['ext', 'int']:
446            raise ValueError, \
447                  r"Coordinate system must be exterior ('ext') or interior ('int')."
448
449        c  = self.first_fundamental_form_coefficient([1,1])
450        if coordinates == 'ext':
451            f1 = self.natural_frame()[1]
452
453            E1 = simplify_vector(f1/sqrt(c))
454            E2 = simplify_vector(self.normal_vector(normalized=True).cross_product(E1))
455        else:
456            E1 =  vector([(1/sqrt(c)).simplify_full(), 0])
457            E2 = simplify_vector(self.rotation(pi/2)*E1)
458        return  {1:E1, 2:E2}
459
460
461    def lie_bracket(self,v,w):
462        """
463        Gives the Lie bracket of two vector fields which are given by coordinates with respect to the natural frame
464
465
466        INPUT:
467        v and w are vectors
468
469
470        OUTPUT:
471        vector [v,w]
472
473
474        EXAMPLES::
475
476            sage: var('u,v')
477            sage: assume(cos(v)>0)
478            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
479            sage: sphere.lie_bracket([u,v],[-v,u])
480            (0, 0)
481
482            sage: EE_int = sphere.orthonormal_frame(coordinates='int')
483            sage: sphere.lie_bracket(EE_int[1],EE_int[2])
484            (sin(v)/cos(v)^2, 0)
485        """
486        vv = vector([xx for xx in v])
487        ww = vector([xx for xx in w])
488        uuu = [self.variables[1],self.variables[2]]
489        Dvv = matrix([ [ diff(xxx,uu).simplify_full() for uu in uuu ] for xxx in vv])
490        Dww = matrix([ [ diff(xxx,uu).simplify_full() for uu in uuu ] for xxx in ww])
491        return simplify_vector(vector(Dvv*ww - Dww*vv))
492
493
494
495    def frame_structure_functions(self,e1,e2):
496        """
497        Gives the structure functions $c^k_{ij}$ for a frame field $e_1$, $e_2$, where
498        $[e_i,e_j] = c^k_{ij}e_k$
499
500
501        INPUT:
502        Frame vectors $e_1$, $e_2$
503
504
505        OUTPUT:
506        The dictionary of $c^k_{ij}$.
507        Warning: the tuple (i,j,k) corresponds to $c^k_{ij}$
508
509
510        EXAMPLES::
511
512            sage: var('u,v')
513            sage: assume(cos(v)>0)
514            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
515            sage: sphere.frame_structure_functions([u,v],[-v,u])
516            {(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}
517            sage: structfun = sphere.frame_structure_functions(EE_int[1],EE_int[2])
518            sage: structfun
519            {(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}
520            sage: sphere.lie_bracket(EE_int[1],EE_int[2]) - CC[(1,2,1)]*EE_int[1] - CC[(1,2,2)]*EE_int[2]
521            (0, 0)
522            """
523        ee1 = vector([xx for xx in e1])
524        ee2 = vector([xx for xx in e2])
525
526        lb = simplify_vector(self.lie_bracket(ee1,ee2))
527        trans_matrix = matrix([[ee1[0],ee2[0]],[ee1[1],ee2[1]]])
528        zz = simplify_vector(trans_matrix.inverse()*lb)
529        return {(1,1,1):0, (1,1,2):0, (1,2,1):zz[0], (1,2,2):zz[1], (2,1,1):-zz[0], (2,1,2):-zz[1],(2,2,1):0, (2,2,2):0}
530
531
532    @cached_method
533    def _compute_second_order_frame_element(self, index):
534        variables = [self.variables[i] for i in index]
535        ddr_element = vector([diff(f, variables).simplify_full() for f in self.equation])
536
537        return ddr_element
538
539    @cached_method
540    def second_order_natural_frame(self):
541
542        """
543        Gives the second derivatives of the equation $\\vec r = \\vec r(u^1,u^2)$ of parametrized surface
544
545
546        INPUT:
547        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
548
549
550        OUTPUT:
551        With empty argument the output is the dictionary of $\partial_{ij}\\vec r(u^1,u^2)$.
552
553        With given indices the output is the corresponding second derivative of the surface parametric equation,
554        For example, index1 = 1, index2 = 2 gives $\partial_{12}\\vec r(u^1,u^2)$.
555
556
557        EXAMPLES::
558
559            sage: var('u,v')
560            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
561            sage: sphere.second_order_natural_frame()
562            {(1, 2): (sin(u)*sin(v), -sin(v)*cos(u), 0), (1, 1): (-cos(u)*cos(v),
563            -sin(u)*cos(v), 0), (2, 1): (sin(u)*sin(v), -sin(v)*cos(u), 0), (2, 2):
564            (-cos(u)*cos(v), -sin(u)*cos(v), -sin(v))}
565
566            sage: sphere.second_order_natural_frame(1,1)
567            (-cos(u)*cos(v), -sin(u)*cos(v), 0)
568            sage: sphere.second_order_natural_frame(1,2)
569            (sin(u)*sin(v), -sin(v)*cos(u), 0)
570            sage: sphere.second_order_natural_frame(2,2)
571            (-cos(u)*cos(v), -sin(u)*cos(v), -sin(v)
572            """
573
574        vectors = {}
575        for index in self.index_list:
576            sorted_index = tuple(sorted(index))
577            vectors[index] = \
578                self._compute_second_order_frame_element(sorted_index)
579        return vectors
580
581    def second_order_natural_frame_element(self, index):
582        index = tuple(sorted(index))
583        if index not in self.index_list:
584            raise ValueError, "Index %s out of bounds." % str(index)
585        return self._compute_second_order_frame_element(index)
586
587
588    @cached_method
589    def _compute_second_fundamental_form_coefficient(self, index):
590        NN = self.normal_vector(normalized=True)
591        v  = self.second_order_natural_frame_element_new(index)
592        return (v*NN).simplify_full()
593
594
595    def second_fundamental_form_coefficient(self, index):
596        index = tuple(index)
597        if index not in self.index_list:
598            raise ValueError, "Index %s out of bounds." % str(index)
599        return self._compute_second_fundamental_form_coefficient(index)
600
601
602    @cached_method
603    def second_fundamental_form_coefficients(self):
604        """
605        This function gives the coefficients $h_{ij}$ of the second fundamental form.
606
607        INPUT:
608        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
609
610        OUTPUT:
611        with empty argument the output is the dictionary of coefficients of the second fundamental form
612
613        with given indices the output is the corresponding coefficient of the second fundamental form,
614        for example index1 = 1, index2 = 2 gives the coefficient $h_{12}$
615
616        EXAMPLES::
617
618            sage: var('u,v')
619            sage: assume(cos(v)>0)
620            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
621
622            sage: sphere.second_fundamental_form_coefficients()
623            {(1, 2): 0, (1, 1): -cos(v)^2, (2, 1): 0, (2, 2): -1}
624            sage: sphere.second_fundamental_form_coefficients(1,1)
625            -cos(v)^2
626            sage: sphere.second_fundamental_form_coefficients(2,1)
627            0
628
629            sage: sphere.second_fundamental_form_coefficients(3,2)
630            'The argument is not appropriate. Read doc'
631        """
632
633        coefficients = {}
634        for index in self.index_list:
635            coefficients[index] = \
636                self._compute_second_fundamental_form_coefficient(index)
637        return coefficients
638
639
640    def second_fundamental_form(self,vector1,vector2):
641        """
642        Finds the value of second fundamental form on two vectors
643
644        INPUT:
645        Two vectors $v=(v^1,v^2)$ and $w=(w^1,w^2)$
646
647        OUTPUT:
648        $h_{11} v^1 w^1 + h_{12}(v^1 w^2 + v^2 w^1) + h_{22} v^2 w^2$
649
650        EXAMPLES::
651
652           sage: var('u,v')
653           sage: var ('v1,v2,w1,w2')
654           sage: assume(cos(v) > 0)
655           sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
656           sage: sphere.second_fundamental_form(vector([v1,v2]),vector([w1,w2]))
657           -v1*w1*cos(v)^2 - v2*w2
658           sage: vv = vector([1,2])
659           sage: sphere.second_fundamental_form(vv,vv)
660           -cos(v)^2 - 4
661           sage: sphere.second_fundamental_form([1,1],[2,1])
662           -2*cos(v)^2 - 1
663
664        """
665        hh = self.second_fundamental_form_coefficients()
666        return sum(hh[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list)
667
668
669    @cached_method
670    def gauss_curvature(self):
671        """
672        Finds the gaussian curvature $K = \\frac{h_{11}h_{22} - h_{12}^2}{g_{11}g_{22} - g_{12}^2}$.
673
674        INPUT:
675        No arguments
676
677        OUTPUT:
678        $K = \\frac{h_{11}h_{22} - h_{12}^2}{g_{11}g_{22} - g_{12}^2}$
679
680        EXAMPLES::
681
682           sage: var('R')
683           sage: assume(R>0)
684           sage: var('u,v')
685           sage: assume(cos(v)>0)
686           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
687           sage: sphere.gauss_curvature()
688           R^(-2)
689
690        """
691        hh = self.second_fundamental_form_coefficients()
692        return ((hh[(1,1)]*hh[(2,2)]-hh[(1,2)]**2)/self.area_form_squared()).simplify_full()
693
694
695    @cached_method
696    def mean_curvature(self):
697        """
698        Finds the mean curvature $H = \\frac{1}{2}\\frac{g_{22}h_{11} - 2g_{12}h_{12} + g_{11}h_{22}}{g_{11}g_{22} - g_{12}^2}$.
699
700        INPUT:
701        No arguments
702
703        OUTPUT:
704        $H = \\frac{1}{2}\\frac{g_{22}h_{11} - 2g_{12}h_{12} + g_{11}h_{22}}{g_{11}g_{22} - g_{12}^2}$.
705
706        EXAMPLES::
707
708           sage: var('R')
709           sage: assume(R>0)
710           sage: var('u,v')
711           sage: assume(cos(v)>0)
712           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
713           sage: sphere.mean_curvature()
714           -1/R
715
716        """
717        gg = self.first_fundamental_form_coefficients()
718        hh = self.second_fundamental_form_coefficients()
719        denom = 2*self.area_form_squared()
720        enum =(gg[(2,2)]*hh[(1,1)]-2*gg[(1,2)]*hh[(1,2)]+gg[(1,1)]*hh[(2,2)]).simplify_full()
721        return (enum/denom).simplify_full()
722
723
724    @cached_method
725    def principal_curvatures(self):
726        """
727        Finds the principal curvatures of the surface
728
729        INPUT:
730        No arguments
731
732        OUTPUT:
733        The dictionary of principal curvatures
734
735        EXAMPLES::
736
737           sage: var('R')
738           sage: assume(R>0)
739           sage: var('u,v')
740           sage: assume(cos(v)>0)
741           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
742           sage: sphere.principal_curvatures()
743           {1: -1/R, 2: -1/R}
744
745           sage: var('u,v')
746           sage: var('R,r')
747           sage: assume(R>r,r>0)
748           sage: torus = parametrized_surface3d([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
749           sage: torus.principal_curvatures()
750           {1: -cos(v)/(r*cos(v) + R), 2: -1/r}
751
752        """
753
754        from sage.symbolic.assumptions import assume
755        from sage.symbolic.relation import solve
756        from sage.calculus.var import var
757
758        KK = self.gauss_curvature()
759        HH = self.mean_curvature()
760
761        # jvkersch: when this assumption is uncommented, Sage raises an error stating that the assumption
762        # is redundant... Can we safely omit this, based on some geometric reasoning?
763
764        # assume(HH**2-KK>=0)
765        # mikarm: This is a problem, I had a lot of trobles here. Of course, this inequality always hold true.
766        # I insert this assumption because sage sometimes, in simplification, uses the complex numbers, though the roots
767        # are, for sure, real.
768        # This, in turn, causes problems when we substitute coordinates into the expression of principal curvatures.
769        # Unfortunately, I did not manage to tell Sage that they are real (to declare the variables as real).
770        # Moreover, in a neighborhood of an umbilic point we even cannot "smoothly" order the set of principal curvatures.
771        # So, in general, at present this method is far from the final form.
772
773
774
775        x = var('x')
776        sol = solve(x**2 -2*HH*x + KK == 0,x)
777
778        #k1=var('k1')
779        #k2=var('k2')
780
781        # jvkersch: when I tried to run the previous version of the code, I ran into the problem that if the equation for the principal curvatures had a double root (as in the case of the sphere example in the worksheet), solve returned only one root.  Maybe this is a difference due to having different versions of sage.
782
783        k1 = (x.substitute(sol[0])).simplify_full()
784        if len(sol) == 1:
785            k2 = k1
786        else:
787            k2 = (x.substitute(sol[1])).simplify_full()
788
789        return {1:k1,2: k2}
790
791    @cached_method
792    def principal_directions(self):
793        """
794        Finds the principal curvatures and principal directions of the surface
795
796        INPUT:
797        No arguments
798
799        OUTPUT:
800        The dictionary of lists [a principal curvature, the corresponding principal direction]
801
802        If principal curvatures coincide, gives the warning that the surface is a sphere.
803
804        EXAMPLES::
805
806           sage: var('u,v')
807           sage: var('R,r')
808           sage: assume(R>r,r>0)
809           sage: torus = parametrized_surface3d([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
810           sage: pdd = torus.principal_directions()
811           sage: pdd[1]
812           [-cos(v)/(r*cos(v) + R), (1, 0)]
813           sage: pdd[2]
814           [-1/r, (0, -(R*r*cos(v) + R^2)/r)]
815
816           sage: var('RR')
817           sage: assume(RR>0)
818           sage: var('u,v')
819           sage: assume(cos(v)>0)
820           sage: sphere = parametrized_surface3d([RR*cos(u)*cos(v),RR*sin(u)*cos(v),RR*sin(v)],[u,v],'sphere')
821           sage: sphere.principal_directions()
822           'This is a sphere, so any direction is principal'
823
824           sage: var('aa')
825           sage: assume(aa>0)
826           sage: catenoid = parametrized_surface3d([aa*cosh(v)*cos(u),aa*cosh(v)*sin(u),v],[u,v],'catenoid')
827           sage: pd = catenoid.principal_directions()
828           sage: pd[1][1]
829           (0, 1/2*(2*aa^3*sinh(v)^2*cosh(v) + sqrt(4*aa^4*sinh(v)^2*cosh(v)^2 + aa^4 + 4*aa^2*cosh(v)^2 - 2*aa^2 + 1)*aa*cosh(v) + (aa^3 + aa)*cosh(v))/(aa^2*sinh(v)^2 + 1)^(3/2))
830           sage: pd[2][1]
831           (1, 0)
832           sage: pd[1][1]*pd[2][1]
833           0
834        """
835        gg = self.first_fundamental_form_coefficients()
836        hh = self.second_fundamental_form_coefficients()
837        kk = self.principal_curvatures()
838        if kk[1]==kk[2]:
839            return "This is a sphere, so any direction is principal"
840        pd1 = simplify_vector([hh[(1,2)]-kk[1]*gg[(1,2)],-hh[(1,1)]+kk[1]*gg[(1,1)]])
841        if pd1==vector([0,0]):
842           pd1 = vector([1,0])
843        pd2 = simplify_vector([hh[(1,2)]-kk[2]*gg[(1,2)],-hh[(1,1)]+kk[2]*gg[(1,1)]])
844        if pd2==vector([0,0]):
845           pd2 = vector([1,0])
846        return {1:[kk[1],pd1],2:[kk[2],pd2]}
847
848
849    @cached_method
850    def connection_coefficients(self):
851        """
852        Finds the connection coefficients of the surface
853
854        INPUT:
855        No arguments
856
857        OUTPUT:
858        The dictionary of connection coefficients.
859
860        Warning: the triple $(i,j,k)$ corresponds to $\Gamma^k_{ij}$.
861
862        EXAMPLES::
863
864           sage: var('r')
865           sage: assume(r > 0)
866           sage: var('u,v')
867           sage: assume(cos(v)>0)
868           sage: sphere = parametrized_surface3d([r*cos(u)*cos(v),r*sin(u)*cos(v),r*sin(v)],[u,v],'sphere')
869           sage: sphere.connection_coefficients()
870           {(1, 2, 1): -sin(v)/cos(v), (2, 2, 2): 0, (1, 2, 2): 0, (2, 1, 1): -sin(v)/cos(v), (1, 1, 2): sin(v)*cos(v), (2, 2, 1): 0, (2, 1, 2): 0, (1, 1, 1): 0}
871
872        """
873        x = self.variables
874        gg = self.first_fundamental_form_coefficients()
875        gi = self.first_fundamental_form_inverse_coefficients()
876
877        dg = {}
878        for kkk in self.index_list_3:
879            dg[kkk]=gg[(kkk[1],kkk[2])].differentiate(x[kkk[0]]).simplify_full()
880        structfun={}
881
882        for kkk in self.index_list_3:
883            structfun[kkk]=sum(gi[(kkk[2],s)]*(dg[(kkk[0],kkk[1],s)]+dg[(kkk[1],kkk[0],s)]-dg[(s,kkk[0],kkk[1])])/2 for s in (1,2)).full_simplify()
884        return structfun
885
886
887    # jvkersch: this private method creates an ode_solver object, which can be used
888    # to integrate the geodesic equations numerically.
889
890    @cached_method
891    def _create_geodesic_ode_system(self):
892        from sage.ext.fast_eval import fast_float
893        from sage.calculus.var import var
894        from sage.gsl.ode import ode_solver
895
896        u1 = self.variables[1]
897        u2 = self.variables[2]
898        v1, v2 = var('v1, v2')
899
900        C = self.connection_coefficients()
901
902        dv1 = - C[(1,1,1)]*v1**2 - 2*C[(1,2,1)]*v1*v2 - C[(2,2,1)]*v2**2
903        dv2 = - C[(1,1,2)]*v1**2 - 2*C[(1,2,2)]*v1*v2 - C[(2,2,2)]*v2**2
904        fun1 = fast_float(dv1, str(u1), str(u2), str(v1), str(v2))
905        fun2 = fast_float(dv2, str(u1), str(u2), str(v1), str(v2))
906
907        geodesic_ode = ode_solver()
908        geodesic_ode.function = \
909                              lambda t, (u1, u2, v1, v2) : \
910                              [v1, v2, fun1(u1, u2, v1, v2), fun2(u1, u2, v1, v2)]
911        return geodesic_ode
912
913
914    # jvkersch: integrate the geodesic equations numerically
915    # mikarm: Very good. I cheked it restarting each time the worksheet, it works much faster.
916
917    def geodesics_numerical(self, p0, v0, tinterval):
918        """
919        This method gives the numerical solution for the geodesic equations
920
921        INPUT:
922        p0 is the list of the coordinates of the initial point
923
924        v0 is the  list of the coordinates of the initial vector
925
926        tinterval is the list [a,b,M], where (a,b) is the domain of the geodesic, M is the number of division points
927
928        OUTPUT:
929        The list consisting of lists [t, [u1(t),u2(t)], [v1(t),v2(t)], [x1(t),x2(t),x3(t)]], where
930
931        t is a subdivision point
932
933        [u1(t),u2(t)] is the list of coordinates of the geodesic point
934
935        [v1(t),v2(t)] is the list of coordinates of the vector tanget to the geodesic
936
937        [x1(t),x2(t),x3(t)] is the list of coordinates of the geodesic point in the three-dimensional space
938
939
940        EXAMPLES::
941
942           sage: var('p,q')
943           sage: v = [p,q]
944           sage: assume(cos(q)>0)
945           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
946           sage: gg_array = sphere.geodesics_numerical([0,0],[1,1],[0,2*pi,5])
947           sage: gg_array[0]
948           [0.0, [0.0, 0.0], [1.0, 1.0], (1, 0, 0)]
949           sage: gg_array[1]
950           [1.2566370614359172, [0.76440104189216407, 1.8586223516062499], [-0.2838683571264714, 1.9194187087799912], (-0.204895333443, 0.692104654602, 0.692104796553)]
951
952
953        """
954
955        u1 = self.variables[1]
956        u2 = self.variables[2]
957
958        solver = self._create_geodesic_ode_system()
959
960        t_interval, n = tinterval[0:2], tinterval[2]
961        solver.y_0 = [p0[0], p0[1], v0[0], v0[1]]
962        solver.ode_solve(t_span=t_interval, num_points=n)
963
964        parsed_solution = \
965          [[vec[0], vec[1][0:2], vec[1][2:], self.eq_callable(vec[1][0], vec[1][1])] \
966           for vec in solver.solution]
967
968        return parsed_solution
969
970
971    @cached_method
972    def _create_pt_ode_system(self, curve, t):
973        from sage.ext.fast_eval import fast_float
974        from sage.calculus.var import var
975        from sage.gsl.ode import ode_solver
976
977        u1 = self.variables[1]
978        u2 = self.variables[2]
979        v1, v2 = var('v1, v2')
980
981        du1 = diff(curve[0], t)
982        du2 = diff(curve[1], t)
983
984        C = self.connection_coefficients()
985        # mikarm: here is the mistake
986        #for coef in C:
987        #    C[coef] = C[coef].subs({u1: curve[0], u2: curve[1]})
988        # This affected the connection coefficients causing them being now functions of t
989        # To tell the truth I do not understand why, because setting C = self.connection.coefficients
990        # and then changing C should not change the output of self.connection.coefficients() but it did
991        # at the same time, I think that
992        # "for coef in C" is wrong here anyway
993
994        C_1 = dict([[coef,C[coef].subs({u1: curve[0], u2: curve[1]})] for coef in self.index_list_3])
995
996        dv1 = - C_1[(1,1,1)]*v1*du1 - C_1[(1,2,1)]*(du1*v2 + du2*v1) - C_1[(2,2,1)]*du2*v2
997        dv2 = - C_1[(1,1,2)]*v1*du1 - C_1[(1,2,2)]*(du1*v2 + du2*v1) - C_1[(2,2,2)]*du2*v2
998        fun1 = fast_float(dv1, str(t), str(v1), str(v2))
999        fun2 = fast_float(dv2, str(t), str(v1), str(v2))
1000        # mikarm: debug
1001        #print dv1,dv2
1002        #
1003        pt_ode = ode_solver()
1004        pt_ode.function = lambda t, (v1, v2): [fun1(t, v1, v2), fun2(t, v1, v2)]
1005        return pt_ode
1006
1007
1008    # mikarm: We should rewrite it like you did for geodesics.
1009    # jvkersch: OK, see this and the above
1010
1011    def parallel_translation_numerical_new(self,curve,t,v0,tinterval):
1012        """
1013        This method gives the numerical solution to the equation of parallel translation of a vector
1014
1015        INPUT:
1016        curve equation = list of functions which determine the curve wrt the local coordinate
1017
1018        t - curve parameter
1019
1020        v0 - initial vector
1021
1022        tinterval = [a,b,N], (a,b) is the domain of the curve, N is the number of subdivision points
1023
1024        OUTPUT:
1025        The list consisting of lists [t, [v1(t),v2(t)]], where
1026
1027        t is a subdivision point
1028
1029        [v1(t),v2(t)] is the list of coordinates of the vector translated parallely along the curve
1030
1031        EXAMPLES::
1032
1033           sage: var('p,q')
1034           sage: v = [p,q]
1035           sage: assume(cos(q)>0)
1036           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
1037           sage: var('ss')
1038           sage: vv_array = sphere.parallel_translation_numerical([ss,ss],ss,[1,1],[0,pi/4,20])
1039           sage: vv_array[0]
1040           [0.0, [1.0, 1.0]]
1041           sage: vv_array[5]
1042           [0.19634954084936207, [0.98060182522638917, 1.0389930474425113]]
1043
1044        """
1045
1046        u1 = self.variables[1]
1047        u2 = self.variables[2]
1048
1049        solver = self._create_pt_ode_system(tuple(curve), t)
1050
1051        t_interval, n = tinterval[0:2], tinterval[2]
1052        solver.y_0 = v0
1053        solver.ode_solve(t_span=t_interval, num_points=n)
1054
1055        return solver.solution
1056
1057
1058
1059
1060
1061    def parallel_translation_numerical(self,curve,t,v0,tinterval):
1062        """
1063        This method gives the numerical solution to the equation of parallel translation of a vector
1064
1065        INPUT:
1066        curve equation = list of functions which determine the curve wrt the local coordinate
1067
1068        t - curve parameter
1069
1070        v0 - initial vector
1071
1072        tinterval = [a,b,N], (a,b) is the domain of the curve, N is the number of subdivision points
1073
1074        OUTPUT:
1075        The list consisting of lists [t, [v1(t),v2(t)]], where
1076
1077        t is a subdivision point
1078
1079        [v1(t),v2(t)] is the list of coordinates of the vector translated parallely along the curve
1080
1081        EXAMPLES::
1082
1083           sage: var('p,q')
1084           sage: v = [p,q]
1085           sage: assume(cos(q)>0)
1086           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
1087           sage: var('ss')
1088           sage: vv_array = sphere.parallel_translation_numerical([ss,ss],ss,[1,1],[0,pi/4,20])
1089           sage: vv_array[0]
1090           [0.0, [1.0, 1.0]]
1091           sage: vv_array[5]
1092           [0.19634954084936207, [0.98060182522638917, 1.0389930474425113]]
1093
1094        """
1095        u1 = self.variables[1]
1096        u2 = self.variables[2]
1097        #var('uu1,uu2')
1098
1099        #f(u1,u2) = self.equation
1100
1101        def f(ux, vx):
1102            return [self.equation[0].subs({u1: ux, u2: vx}),
1103                    self.equation[1].subs({u1: ux, u2: vx}),
1104                    self.equation[2].subs({u1: ux, u2: vx})]
1105
1106        #from sage.calculus.var import var
1107        #t = var('t')
1108        tt = t
1109
1110        C111 = self.connection_coefficients()[(1,1,1)]
1111        C121 = self.connection_coefficients()[(1,2,1)]
1112        C221 = self.connection_coefficients()[(2,2,1)]
1113        C112 = self.connection_coefficients()[(1,1,2)]
1114        C122 = self.connection_coefficients()[(1,2,2)]
1115        C222 = self.connection_coefficients()[(2,2,2)]
1116
1117        du1=diff(curve[0],tt)
1118        du2=diff(curve[1],tt)
1119        uu1=curve[0]
1120        uu2=curve[1]
1121
1122
1123        def ode_system(unknown_functions,t1):
1124            du1p = du1.subs(tt==t1)
1125            du2p = du2.subs(tt==t1)
1126            uu1p = uu1.subs(tt==t1)
1127            uu2p = uu2.subs(tt==t1)
1128
1129            c111 = C111.subs({u1: uu1p, u2: uu2p})
1130            c121 = C121.subs({u1: uu1p, u2: uu2p})
1131            c221 = C221.subs({u1: uu1p, u2: uu2p})
1132            c112 = C112.subs({u1: uu1p, u2: uu2p})
1133            c122 = C122.subs({u1: uu1p, u2: uu2p})
1134            c222 = C222.subs({u1: uu1p, u2: uu2p})
1135
1136            #print c111, c121, c221, c112, c122, c222
1137
1138            v1,v2 = unknown_functions
1139            f1 = - c111*du1p*v1 - c121*(du1p*v2 + du2p*v1) - c221*du2p*v2
1140            f2 = - c112*du1p*v1 - c122*(du1p*v2 + du2p*v1) - c222*du2p*v2
1141            dv1 = f1.subs(tt==t1)
1142            dv2 = f2.subs(tt==t1)
1143
1144            #print dv1, dv2
1145
1146            return float(dv1), float(dv2)
1147
1148        step = float((tinterval[1]-tinterval[0])/tinterval[2])
1149        ttt =  float(tinterval[0])
1150        tarray = [ ttt ]
1151        for counter in range(1,tinterval[2]+1):
1152           ttt = ttt + step
1153           tarray = tarray + [ ttt ]
1154
1155
1156        # import ode solving routine
1157        import scipy.integrate
1158
1159        initial_data = (v0[0],v0[1])
1160        solution = scipy.integrate.odeint(ode_system,initial_data,tarray)
1161        return [[ tarray[counter],[solution[counter,0], solution[counter,1]]]   for counter in range(0,len(tarray))]
1162