• Sign Up
• Login
• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# geodesic_on_sphere.py

a guest Oct 25th, 2015 317 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import sympy as sym
2. import numpy as np
3. import scipy.integrate as sc
4.
5. def Christoffel2nd_sphere():
6.     from sympy.abc import u,v
7.     from sympy import tan, cos ,sin
8.     # found using sympy.diffgeom
9.     return sym.Matrix([[(0, -2*tan(v)),  (0, 0)],[(sin(v)*cos(v), 0), (0, 0)]])
10.
11. def display_sphere_with_geodesic(u0,s0,s1,ds):
12.     C = Christoffel2nd_sphere()
13.     X = solve(C,u0,s0,s1,ds)
14.     import matplotlib.pylab as plt
15.     from mpl_toolkits.mplot3d import Axes3D
16.     N = X[:,0].shape[0]
17.     u,v = plt.meshgrid(np.linspace(0,2*np.pi,N),np.linspace(0,2*np.pi,N))
18.     x = np.cos(u)*np.cos(v)
19.     y = np.sin(u)*np.cos(v)
20.     z = np.sin(v)
21.
22.     fig = plt.figure()
23.     ax = fig.add_subplot(111, projection='3d')
24.     ax.view_init(elev=50, azim=-120)
25.     # use transparent colormap
26.     import matplotlib.cm as cm
27.     theCM = cm.get_cmap()
28.     theCM._init()
29.     alphas = -.5*np.ones(theCM.N)
30.     theCM._lut[:-3,-1] = alphas
31.     ax.plot_surface(x,y,z,linewidth=0,cmap=theCM)
32.     plt.hold('on')
33.
34.     # plot the parametrized data on to the sphere
35.     u,v = X[:,0], X[:,2]
36.     x = np.cos(u)*np.cos(v)
37.     y = np.sin(u)*np.cos(v)
38.     z = np.sin(v)
39.
40.     ax.plot(x,y,z,'--r')
41.     from math import pi
42.     s1_ = s1/pi
43.     fig.suptitle('$s\in[%.1f\, , \,%2.1f\pi]$'%(s0,s1_))
44.     plt.show()
45.
46.
47. def f(y,s,C,u,v):
48.     y0 = y[0] # u
49.     y1 = y[1] # u'
50.     y2 = y[2] # v
51.     y3 = y[3] # v'
52.     dy = np.zeros_like(y)
53.     dy[0] = y1
54.     dy[2] = y3
55.
56.     C = C.subs({u:y0,v:y2})
57.
58.     dy[1] = -C[0,0][0]*dy[0]**2 -\
59.            2*C[0,0][1]*dy[0]*dy[2] -\
60.              C[0,1][1]*dy[2]**2
61.     dy[3] = -C[1,0][0]*dy[0]**2 -\
62.            2*C[1,0][1]*dy[0]*dy[2] -\
63.              C[1,1][1]*dy[2]**2
64.     return dy
65.
66. def solve(C,u0,s0,s1,ds):
67.     s = np.arange(s0,s1+ds,ds)
68.     from sympy.abc import u,v
69.     return sc.odeint(f,u0,s,args=(C,u,v))
70.
71. from math import pi
72. u0 = [0,.1,0,.1]
73. s0 = 0
74. s1 = 18*pi
75. ds = 0.15
76. display_sphere_with_geodesic(u0,s0,s1,ds)
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top