In [9]:
import numpy as np
import pylab as pl
import os
from mpl_toolkits.mplot3d import Axes3D
path = '/data/'
data = []
for f in os.listdir(os.getcwd()+path):
r = np.loadtxt(os.getcwd()+path+f)
data.append(r)
In our data set the first column corresponds to each iteration of the trayectory computed. The second column corresponds to $t$ which is the relativistic time, the second, third and fourth columns corresponds to the spatial coordinates $x(t)$, $y(t)$ and $z(t)$ respectively.
Let's plot $t$ vs $x(t)$ of the third file which corresponds to the Innermost Stable Circular Orbit (ISCO).
In [26]:
# create fig
pl.rc('font',size=16)
pl.figure(figsize=(16,8))
#set plot limits
pl.xlim([0,1000])
pl.ylim([-6,6])
#set labels
pl.title('Oscilation of a particle in the ISCO')
pl.xlabel('$t$')
pl.ylabel('$X(t)$')
pl.legend(loc=4,fancybox=True, shadow=True)
#set t & x(t)
t = data[2][:,1]
x = data[2][:,2]
#plot
pl.plot(t,x,'r', lw=1.0, label='$n=1$')
pl.show()
Let's plot the spatial trayectories.
In [28]:
#create figure
pl.rc('font',size=10)
fig = pl.figure(figsize=(9,8),dpi=600)
ax = fig.gca(projection='3d')
#draw lines
c = ['b','g','r','c','k']
labels = ['$r=4M$','$r=5M$','$ISCO$','$r=6M$','$r=7M$']
for i, d in enumerate(data):
ax.plot(d[:,2],d[:,3],d[:,4],c[i], label=labels[i])
#create and plot a sphere
r = 2
pi = np.pi
cos = np.cos
sin = np.sin
phi,theta = np.mgrid[0.0:pi:100j,0.0:2.0*pi:100j]
x = r*sin(phi)*cos(theta)
y = r*sin(phi)*sin(theta)
z = r*cos(phi)
ax.plot_surface(x,y,z, rstride=1,cstride=1,color='k', alpha=0.6,linewidth=0)
#set plot limits
ax.set_xlim3d([-8,8])
ax.set_ylim3d([-8,8])
ax.set_zlim3d([-7,7])
#labels
pl.title('Timelike geodesics in the Schwarzschild equatorial plane')
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
ax.set_zlabel('Z axis')
pl.legend(loc=4,fancybox=True, shadow=True)
pl.show()
In order to interact with the 3D plot, you should run the file "02 - 3D Plot.py"
In [ ]: