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. OK, I Understand
 
Top