Ticket #10132: parametrized_surface3d.3.py

File parametrized_surface3d.3.py, 40.4 KB (added by mikarm, 10 years ago)

Next version of ParametrizedSurface3D class

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#
12#  Distributed under the terms of the GNU General Public License (GPL)
13#                  http://www.gnu.org/licenses/
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        v  = self.second_order_natural_frame_element(index) 
593        return (v*NN).simplify_full()
594
595       
596    def second_fundamental_form_coefficient(self, index):
597        index = tuple(index)
598        if index not in self.index_list:
599            raise ValueError, "Index %s out of bounds." % str(index)
600        return self._compute_second_fundamental_form_coefficient(index)
601
602
603    @cached_method
604    def second_fundamental_form_coefficients(self):
605        """
606        This function gives the coefficients $h_{ij}$ of the second fundamental form.
607       
608        INPUT:
609        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
610
611        OUTPUT:
612        with empty argument the output is the dictionary of coefficients of the second fundamental form
613
614        with given indices the output is the corresponding coefficient of the second fundamental form,
615        for example index1 = 1, index2 = 2 gives the coefficient $h_{12}$
616
617        EXAMPLES::
618
619            sage: var('u,v')
620            sage: assume(cos(v)>0)
621            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
622
623            sage: sphere.second_fundamental_form_coefficients()
624            {(1, 2): 0, (1, 1): -cos(v)^2, (2, 1): 0, (2, 2): -1}
625            sage: sphere.second_fundamental_form_coefficients(1,1)
626            -cos(v)^2
627            sage: sphere.second_fundamental_form_coefficients(2,1)
628            0
629
630            sage: sphere.second_fundamental_form_coefficients(3,2)
631            'The argument is not appropriate. Read doc'
632        """
633
634        coefficients = {}
635        for index in self.index_list:
636            coefficients[index] = \
637                self._compute_second_fundamental_form_coefficient(index)
638        return coefficients
639
640
641    def second_fundamental_form(self,vector1,vector2):
642        """
643        Finds the value of second fundamental form on two vectors
644
645        INPUT:
646        Two vectors $v=(v^1,v^2)$ and $w=(w^1,w^2)$
647
648        OUTPUT:
649        $h_{11} v^1 w^1 + h_{12}(v^1 w^2 + v^2 w^1) + h_{22} v^2 w^2$
650
651        EXAMPLES::
652
653           sage: var('u,v')
654           sage: var ('v1,v2,w1,w2')
655           sage: assume(cos(v) > 0)
656           sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
657           sage: sphere.second_fundamental_form(vector([v1,v2]),vector([w1,w2]))
658           -v1*w1*cos(v)^2 - v2*w2
659           sage: vv = vector([1,2])
660           sage: sphere.second_fundamental_form(vv,vv)
661           -cos(v)^2 - 4
662           sage: sphere.second_fundamental_form([1,1],[2,1])
663           -2*cos(v)^2 - 1
664                   
665        """
666        hh = self.second_fundamental_form_coefficients()
667        return sum(hh[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list)
668
669
670    @cached_method
671    def gauss_curvature(self):
672        """
673        Finds the gaussian curvature $K = \\frac{h_{11}h_{22} - h_{12}^2}{g_{11}g_{22} - g_{12}^2}$.
674
675        INPUT:
676        No arguments
677
678        OUTPUT:
679        $K = \\frac{h_{11}h_{22} - h_{12}^2}{g_{11}g_{22} - g_{12}^2}$
680
681        EXAMPLES::
682
683           sage: var('R')
684           sage: assume(R>0)
685           sage: var('u,v')
686           sage: assume(cos(v)>0)
687           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
688           sage: sphere.gauss_curvature()
689           R^(-2)
690                   
691        """
692        hh = self.second_fundamental_form_coefficients()
693        return ((hh[(1,1)]*hh[(2,2)]-hh[(1,2)]**2)/self.area_form_squared()).simplify_full()
694
695
696    @cached_method
697    def mean_curvature(self):
698        """
699        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}$.
700
701        INPUT:
702        No arguments
703
704        OUTPUT:
705        $H = \\frac{1}{2}\\frac{g_{22}h_{11} - 2g_{12}h_{12} + g_{11}h_{22}}{g_{11}g_{22} - g_{12}^2}$.
706
707        EXAMPLES::
708           
709           sage: var('R')
710           sage: assume(R>0)
711           sage: var('u,v')
712           sage: assume(cos(v)>0)
713           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
714           sage: sphere.mean_curvature()
715           -1/R
716
717        """
718        gg = self.first_fundamental_form_coefficients()
719        hh = self.second_fundamental_form_coefficients()
720        denom = 2*self.area_form_squared()
721        enum =(gg[(2,2)]*hh[(1,1)]-2*gg[(1,2)]*hh[(1,2)]+gg[(1,1)]*hh[(2,2)]).simplify_full()
722        return (enum/denom).simplify_full()
723
724
725    @cached_method
726    def principal_curvatures(self):
727        """
728        Finds the principal curvatures of the surface
729
730        INPUT:
731        No arguments
732
733        OUTPUT:
734        The dictionary of principal curvatures         
735
736        EXAMPLES::
737                     
738           sage: var('R')
739           sage: assume(R>0)
740           sage: var('u,v')
741           sage: assume(cos(v)>0)
742           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
743           sage: sphere.principal_curvatures()
744           {1: -1/R, 2: -1/R}
745
746           sage: var('u,v')
747           sage: var('R,r')
748           sage: assume(R>r,r>0)
749           sage: torus = parametrized_surface3d([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
750           sage: torus.principal_curvatures()
751           {1: -cos(v)/(r*cos(v) + R), 2: -1/r}
752
753        """
754
755        from sage.symbolic.assumptions import assume
756        from sage.symbolic.relation import solve
757        from sage.calculus.var import var
758
759        KK = self.gauss_curvature()
760        HH = self.mean_curvature()
761
762        # jvkersch: when this assumption is uncommented, Sage raises an error stating that the assumption
763        # is redundant... Can we safely omit this, based on some geometric reasoning?
764       
765        # assume(HH**2-KK>=0)
766        # mikarm: This is a problem, I had a lot of trobles here. Of course, this inequality always hold true.
767        # I insert this assumption because sage sometimes, in simplification, uses the complex numbers, though the roots
768        # are, for sure, real. 
769        # This, in turn, causes problems when we substitute coordinates into the expression of principal curvatures.
770        # Unfortunately, I did not manage to tell Sage that they are real (to declare the variables as real).
771        # Moreover, in a neighborhood of an umbilic point we even cannot "smoothly" order the set of principal curvatures.
772        # So, in general, at present this method is far from the final form.
773         
774
775 
776        x = var('x')
777        sol = solve(x**2 -2*HH*x + KK == 0,x)
778       
779        #k1=var('k1')
780        #k2=var('k2')
781
782        # 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.
783
784        k1 = (x.substitute(sol[0])).simplify_full()
785        if len(sol) == 1:
786            k2 = k1
787        else:
788            k2 = (x.substitute(sol[1])).simplify_full()
789       
790        return {1:k1,2: k2}
791
792        #mikarm: add method for finding shape operator
793    @cached_method
794    def shape_operator_coefficients(self):
795
796        gi = self.first_fundamental_form_inverse_coefficients()
797        hh = self.second_fundamental_form_coefficients()
798
799        shop11 = (gi[(1,1)]*hh[(1,1)] + gi[(1,2)]*hh[(1,2)]).simplify_full() 
800        shop12 = (gi[(1,1)]*hh[(2,1)] + gi[(1,2)]*hh[(2,2)]).simplify_full() 
801        shop21 = (gi[(2,1)]*hh[(1,1)] + gi[(2,2)]*hh[(1,2)]).simplify_full() 
802        shop22 = (gi[(2,1)]*hh[(2,1)] + gi[(2,2)]*hh[(2,2)]).simplify_full() 
803
804        return {(1,1):shop11,(1,2):shop12,(2,1):shop21,(2,2):shop22}
805
806    def shape_operator(self,v):
807       vv = vector([xx for xx in v])
808       shop = self.shape_operator_coefficients()
809       shop_matrix=matrix([[shop[(1,1)],shop[(1,2)]],[shop[(2,1)],shop[(2,2)]]])
810       return simplify_vector(shop_matrix*vv) 
811
812    @cached_method
813    def principal_directions(self):
814        """
815        Finds the principal curvatures and principal directions of the surface
816
817        INPUT:
818        No arguments
819
820        OUTPUT:
821        The dictionary of lists [a principal curvature, the corresponding principal direction]         
822
823        If principal curvatures coincide, gives the warning that the surface is a sphere.
824
825        EXAMPLES::
826
827           sage: var('u,v')
828           sage: var('R,r')
829           sage: assume(R>r,r>0)
830           sage: torus = parametrized_surface3d([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
831           sage: pdd = torus.principal_directions()
832           sage: pdd[1]
833           [-cos(v)/(r*cos(v) + R), (1, 0)]
834           sage: pdd[2]
835           [-1/r, (0, -(R*r*cos(v) + R^2)/r)]
836
837           sage: var('RR')
838           sage: assume(RR>0)
839           sage: var('u,v')
840           sage: assume(cos(v)>0)
841           sage: sphere = parametrized_surface3d([RR*cos(u)*cos(v),RR*sin(u)*cos(v),RR*sin(v)],[u,v],'sphere')
842           sage: sphere.principal_directions()
843           'This is a sphere, so any direction is principal'
844
845           sage: var('aa')
846           sage: assume(aa>0)
847           sage: catenoid = parametrized_surface3d([aa*cosh(v)*cos(u),aa*cosh(v)*sin(u),v],[u,v],'catenoid')
848           sage: pd = catenoid.principal_directions()
849           sage: pd[1][1]
850           (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))
851           sage: pd[2][1]
852           (1, 0)
853           sage: pd[1][1]*pd[2][1]
854           0
855        """
856        gg = self.first_fundamental_form_coefficients()
857        hh = self.second_fundamental_form_coefficients()
858        kk = self.principal_curvatures()
859        if kk[1]==kk[2]:
860            return "This is a sphere, so any direction is principal"
861        pd1 = simplify_vector([hh[(1,2)]-kk[1]*gg[(1,2)],-hh[(1,1)]+kk[1]*gg[(1,1)]])
862        if pd1==vector([0,0]): 
863           pd1 = vector([1,0]) 
864        pd2 = simplify_vector([hh[(1,2)]-kk[2]*gg[(1,2)],-hh[(1,1)]+kk[2]*gg[(1,1)]])
865        if pd2==vector([0,0]):
866           pd2 = vector([1,0]) 
867        return {1:[kk[1],pd1],2:[kk[2],pd2]}
868
869
870    @cached_method
871    def connection_coefficients(self):
872        """
873        Finds the connection coefficients of the surface
874
875        INPUT:
876        No arguments
877
878        OUTPUT:
879        The dictionary of connection coefficients.
880
881        Warning: the triple $(i,j,k)$ corresponds to $\Gamma^k_{ij}$. 
882
883        EXAMPLES::
884                                 
885           sage: var('r')
886           sage: assume(r > 0)
887           sage: var('u,v')
888           sage: assume(cos(v)>0)
889           sage: sphere = parametrized_surface3d([r*cos(u)*cos(v),r*sin(u)*cos(v),r*sin(v)],[u,v],'sphere')
890           sage: sphere.connection_coefficients()
891           {(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}
892           
893        """
894        x = self.variables
895        gg = self.first_fundamental_form_coefficients()
896        gi = self.first_fundamental_form_inverse_coefficients()
897
898        dg = {}
899        for kkk in self.index_list_3:
900            dg[kkk]=gg[(kkk[1],kkk[2])].differentiate(x[kkk[0]]).simplify_full()
901        structfun={}
902
903        for kkk in self.index_list_3:
904            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()
905        return structfun
906
907
908    # jvkersch: this private method creates an ode_solver object, which can be used
909    # to integrate the geodesic equations numerically.
910   
911    @cached_method
912    def _create_geodesic_ode_system(self):
913        from sage.ext.fast_eval import fast_float
914        from sage.calculus.var import var
915        from sage.gsl.ode import ode_solver
916
917        u1 = self.variables[1]
918        u2 = self.variables[2]
919        v1, v2 = var('v1, v2')
920
921        C = self.connection_coefficients()
922       
923        dv1 = - C[(1,1,1)]*v1**2 - 2*C[(1,2,1)]*v1*v2 - C[(2,2,1)]*v2**2
924        dv2 = - C[(1,1,2)]*v1**2 - 2*C[(1,2,2)]*v1*v2 - C[(2,2,2)]*v2**2
925        fun1 = fast_float(dv1, str(u1), str(u2), str(v1), str(v2))
926        fun2 = fast_float(dv2, str(u1), str(u2), str(v1), str(v2))
927
928        geodesic_ode = ode_solver()
929        geodesic_ode.function = \
930                              lambda t, (u1, u2, v1, v2) : \
931                              [v1, v2, fun1(u1, u2, v1, v2), fun2(u1, u2, v1, v2)]
932        return geodesic_ode
933
934
935    # jvkersch: integrate the geodesic equations numerically
936    # mikarm: Very good. I cheked it restarting each time the worksheet, it works much faster.
937       
938    def geodesics_numerical(self, p0, v0, tinterval):
939        """
940        This method gives the numerical solution for the geodesic equations 
941
942        INPUT:
943        p0 is the list of the coordinates of the initial point
944       
945        v0 is the  list of the coordinates of the initial vector
946
947        tinterval is the list [a,b,M], where (a,b) is the domain of the geodesic, M is the number of division points
948
949        OUTPUT:
950        The list consisting of lists [t, [u1(t),u2(t)], [v1(t),v2(t)], [x1(t),x2(t),x3(t)]], where
951
952        t is a subdivision point
953
954        [u1(t),u2(t)] is the list of coordinates of the geodesic point
955         
956        [v1(t),v2(t)] is the list of coordinates of the vector tanget to the geodesic 
957
958        [x1(t),x2(t),x3(t)] is the list of coordinates of the geodesic point in the three-dimensional space
959
960
961        EXAMPLES::
962
963           sage: var('p,q')
964           sage: v = [p,q]
965           sage: assume(cos(q)>0)
966           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
967           sage: gg_array = sphere.geodesics_numerical([0,0],[1,1],[0,2*pi,5])
968           sage: gg_array[0]
969           [0.0, [0.0, 0.0], [1.0, 1.0], (1, 0, 0)]
970           sage: gg_array[1]
971           [1.2566370614359172, [0.76440104189216407, 1.8586223516062499], [-0.2838683571264714, 1.9194187087799912], (-0.204895333443, 0.692104654602, 0.692104796553)]
972                                           
973         
974        """ 
975
976        u1 = self.variables[1]
977        u2 = self.variables[2]
978
979        solver = self._create_geodesic_ode_system()
980           
981        t_interval, n = tinterval[0:2], tinterval[2]
982        solver.y_0 = [p0[0], p0[1], v0[0], v0[1]]
983        solver.ode_solve(t_span=t_interval, num_points=n)
984
985        parsed_solution = \
986          [[vec[0], vec[1][0:2], vec[1][2:], self.eq_callable(vec[1][0], vec[1][1])] \
987           for vec in solver.solution]
988
989        return parsed_solution
990
991
992    @cached_method
993    def _create_pt_ode_system(self, curve, t):
994        from sage.ext.fast_eval import fast_float
995        from sage.calculus.var import var
996        from sage.gsl.ode import ode_solver
997
998        u1 = self.variables[1]
999        u2 = self.variables[2]
1000        v1, v2 = var('v1, v2')
1001
1002        du1 = diff(curve[0], t)
1003        du2 = diff(curve[1], t)
1004
1005        C = self.connection_coefficients()
1006        for coef in C:
1007            C[coef] = C[coef].subs({u1: curve[0], u2: curve[1]})
1008       
1009        dv1 = - C[(1,1,1)]*v1*du1 - C[(1,2,1)]*(du1*v2 + du2*v1) - C[(2,2,1)]*du2*v2
1010        dv2 = - C[(1,1,2)]*v1*du1 - C[(1,2,2)]*(du1*v2 + du2*v1) - C[(2,2,2)]*du2*v2
1011        fun1 = fast_float(dv1, str(t), str(v1), str(v2))
1012        fun2 = fast_float(dv2, str(t), str(v1), str(v2))
1013
1014        pt_ode = ode_solver()
1015        pt_ode.function = lambda t, (v1, v2): [fun1(t, v1, v2), fun2(t, v1, v2)]
1016        return pt_ode
1017
1018
1019    # mikarm: We should rewrite it like you did for geodesics.
1020    # jvkersch: OK, see this and the above
1021    def parallel_translation_numerical_new(self,curve,t,v0,tinterval): 
1022        """
1023        This method gives the numerical solution to the equation of parallel translation of a vector
1024
1025        INPUT:
1026        curve equation = list of functions which determine the curve wrt the local coordinate
1027
1028        t - curve parameter
1029
1030        v0 - initial vector
1031
1032        tinterval = [a,b,N], (a,b) is the domain of the curve, N is the number of subdivision points
1033
1034        OUTPUT:
1035        The list consisting of lists [t, [v1(t),v2(t)]], where
1036
1037        t is a subdivision point
1038
1039        [v1(t),v2(t)] is the list of coordinates of the vector translated parallely along the curve
1040
1041        EXAMPLES::
1042
1043           sage: var('p,q')
1044           sage: v = [p,q]
1045           sage: assume(cos(q)>0)
1046           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
1047           sage: var('ss')
1048           sage: vv_array = sphere.parallel_translation_numerical([ss,ss],ss,[1,1],[0,pi/4,20])
1049           sage: vv_array[0]
1050           [0.0, [1.0, 1.0]]
1051           sage: vv_array[5]
1052           [0.19634954084936207, [0.98060182522638917, 1.0389930474425113]]
1053
1054        """
1055
1056        u1 = self.variables[1]
1057        u2 = self.variables[2]
1058
1059        solver = self._create_pt_ode_system(tuple(curve), t)
1060           
1061        t_interval, n = tinterval[0:2], tinterval[2]
1062        solver.y_0 = v0
1063        solver.ode_solve(t_span=t_interval, num_points=n)
1064
1065        return solver.solution
1066
1067
1068
1069
1070
1071    def parallel_translation_numerical(self,curve,t,v0,tinterval): 
1072        """
1073        This method gives the numerical solution to the equation of parallel translation of a vector
1074
1075        INPUT:
1076        curve equation = list of functions which determine the curve wrt the local coordinate
1077
1078        t - curve parameter
1079
1080        v0 - initial vector
1081
1082        tinterval = [a,b,N], (a,b) is the domain of the curve, N is the number of subdivision points
1083
1084        OUTPUT:
1085        The list consisting of lists [t, [v1(t),v2(t)]], where
1086
1087        t is a subdivision point
1088
1089        [v1(t),v2(t)] is the list of coordinates of the vector translated parallely along the curve
1090
1091        EXAMPLES::
1092
1093           sage: var('p,q')
1094           sage: v = [p,q]
1095           sage: assume(cos(q)>0)
1096           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
1097           sage: var('ss')
1098           sage: vv_array = sphere.parallel_translation_numerical([ss,ss],ss,[1,1],[0,pi/4,20])
1099           sage: vv_array[0]
1100           [0.0, [1.0, 1.0]]
1101           sage: vv_array[5]
1102           [0.19634954084936207, [0.98060182522638917, 1.0389930474425113]]
1103
1104        """
1105        u1 = self.variables[1]
1106        u2 = self.variables[2]
1107        #var('uu1,uu2')
1108       
1109        #f(u1,u2) = self.equation
1110
1111        def f(ux, vx):
1112            return [self.equation[0].subs({u1: ux, u2: vx}),
1113                    self.equation[1].subs({u1: ux, u2: vx}),
1114                    self.equation[2].subs({u1: ux, u2: vx})]
1115
1116        #from sage.calculus.var import var
1117        #t = var('t')
1118        tt = t
1119
1120        C111 = self.connection_coefficients()[(1,1,1)] 
1121        C121 = self.connection_coefficients()[(1,2,1)] 
1122        C221 = self.connection_coefficients()[(2,2,1)] 
1123        C112 = self.connection_coefficients()[(1,1,2)] 
1124        C122 = self.connection_coefficients()[(1,2,2)] 
1125        C222 = self.connection_coefficients()[(2,2,2)] 
1126
1127        du1=diff(curve[0],tt)
1128        du2=diff(curve[1],tt)
1129        uu1=curve[0]
1130        uu2=curve[1]
1131
1132
1133        def ode_system(unknown_functions,t1):
1134            du1p = du1.subs(tt==t1)
1135            du2p = du2.subs(tt==t1)
1136            uu1p = uu1.subs(tt==t1)
1137            uu2p = uu2.subs(tt==t1)
1138
1139            c111 = C111.subs({u1: uu1p, u2: uu2p})
1140            c121 = C121.subs({u1: uu1p, u2: uu2p})
1141            c221 = C221.subs({u1: uu1p, u2: uu2p})
1142            c112 = C112.subs({u1: uu1p, u2: uu2p})
1143            c122 = C122.subs({u1: uu1p, u2: uu2p})
1144            c222 = C222.subs({u1: uu1p, u2: uu2p})
1145
1146            #print c111, c121, c221, c112, c122, c222
1147           
1148            v1,v2 = unknown_functions
1149            f1 = - c111*du1p*v1 - c121*(du1p*v2 + du2p*v1) - c221*du2p*v2
1150            f2 = - c112*du1p*v1 - c122*(du1p*v2 + du2p*v1) - c222*du2p*v2
1151            dv1 = f1.subs(tt==t1)
1152            dv2 = f2.subs(tt==t1)
1153
1154            #print dv1, dv2
1155           
1156            return float(dv1), float(dv2) 
1157
1158        step = float((tinterval[1]-tinterval[0])/tinterval[2])
1159        ttt =  float(tinterval[0]) 
1160        tarray = [ ttt ]
1161        for counter in range(1,tinterval[2]+1):
1162           ttt = ttt + step
1163           tarray = tarray + [ ttt ] 
1164
1165
1166        # import ode solving routine
1167        import scipy.integrate
1168
1169        initial_data = (v0[0],v0[1])
1170        solution = scipy.integrate.odeint(ode_system,initial_data,tarray)
1171        return [[ tarray[counter],[solution[counter,0], solution[counter,1]]]   for counter in range(0,len(tarray))]
1172