This workbook describes the: Trajectory of Light Signal using Jacobi’s Elliptic Function. The light deflection arround the sun is determined.
The theory behind the method is found in paper: 'Equations of Orbits, Deflection of Light, Mercury’s Perihelion Shift, Period of Revolution and Time Delay of Signals in Schwarzchild’s Geometry. Some closed-form solutions using Weinberg’s approach. New Version (19/08/2016) Solomon M. Antoniou'. Researchgate
This project was started based upon a topic on: Wetenschapsforum.
Here follows a summary of the result found on forum Wetenschapsforum. The mean set of formulas required are found in article. Where $r_s$ is Schwarzschild radius and $r_0$ radius lightbeam from center mass and $\varphi$ is the angle following the light signal. The radius $r$ with angle $\varphi$ are computed with $(6)$. Finally the startcondition $\sigma$ and $\varphi$ range are computed:
$$e_1 = \frac{r_0 - r_s + \sqrt{(r_0 - r_s)(r_0 + 3r_s)}}{2 r_s r_0} \tag{1}$$$$e_2 = \frac{1}{r_0} \tag{2}$$$$e_3 = \frac{r_0 - r_s - \sqrt{(r_0 - r_s)(r_0 + 3r_s)}}{2 r_s r_0} \tag{3}$$$$\tau = \sqrt{ \frac{ r_s (e_1 - e_3)}{ 4 }} \tag{4}$$$$h = \sqrt{ \frac{ e_2 - e_3}{ e_1 - e_3 }} \tag{5}$$$$\boxed{ r = \frac{1}{ e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \varphi + \sigma ; h^2)}} \tag{6}$$Start conditions for $r=r_{0}$ at $\varphi=\pi/2$ determine $\sigma$ substitude $(2)$:
$$r_{0} = \frac{1}{ e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \frac{\pi}{2} + \sigma ; h^2)}=\frac{1}{e_2} \tag{7}$$After simplification:
$$\mathrm{sn}( \tau \frac{\pi}{2} + \sigma ; h^2) = 1 \tag{8}$$Property elliptic integrals:
$$\mathrm{sn}(K,h^2)=\sin \left( \frac{\pi}{2} \right) \tag{10} $$Where $\tau$ and $h^2$, are given and the startcondition $\sigma$ can be solved with standard form complete elliptic integral Wiki. In the simulation the following integral solver was used from scipy: [special.ellipk]. Definition of complete Jacobi elliptic integral:
$$K(h^2)=\int_{0}^{\frac{\pi}{2}} \frac{d t}{\sqrt{1-h^2 \sin^2(t) }} =\tau \frac{\pi}{2} + \sigma \tag{9} $$Not all solutions are valid, the root of $(6)$ van have zeros. The $x$-range is found by searching the roots of $(6)$.
$$e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \varphi + \sigma ; h^2) = 0 \tag{12}$$$$ \mathrm{sn}( \tau \varphi + \sigma ; h^2) =\sqrt{ -\frac{e_3}{e_2 - e_3}} \tag{13}$$$$\mathrm{sn}(K|h^2)=\sin \left( \phi \right) \tag{14} $$$$\phi=\arcsin \left( \sqrt{ -\frac{e_3}{e_2 - e_3}} \right) \tag{15}$$Where $\tau$, $h^2$, $\phi$ and $\sigma$ are given and can be solved with incomplete elliptic integral [scipy: special.ellipkinc]. The maximum angle $\varphi$ can be calulated with:
$$K(\phi|h^2)=\int_{0}^{\phi} \frac{d t}{\sqrt{1-h^2 \sin^2(t) }} =\tau \varphi + \sigma \tag{16} $$#Version 25.0 Schwarzschild, Jacobi Elliptic
#https://www.wetenschapsforum.nl/viewtopic.php?f=85&t=212427
#Open pyplot in separate interactive window
#from IPython import get_ipython
#get_ipython().run_line_magic('matplotlib', 'qt5')
import numpy as np
import matplotlib.pyplot as plt
#import matplotlib.gridspec as gridspec
import scipy.special as sps
from scipy import optimize
from matplotlib.widgets import Slider, CheckButtons
import matplotlib.ticker as ticker
%matplotlib widget
fig= plt.figure(figsize=(9, 5))
ax1 = plt.subplot2grid((3, 3), (0, 0),rowspan=2)
ax2 = plt.subplot2grid((3, 3), (0, 1),rowspan=2)
ax3 = plt.subplot2grid((3, 3), (0, 2),rowspan=2)
ax4 = plt.subplot2grid((3, 3), (2, 0),colspan=3)
ax4.axis('off')
plt.suptitle('Schwarzschild, Jacobi Elliptic', fontsize=20)
#fig.set_tight_layout(True)
plt.tight_layout(pad=1.5, w_pad=1, h_pad=1)
# Draw slider alpha
axcolor = 'lightgoldenrodyellow'
axalpha = plt.axes([0.25, 0.01,0.4, 0.04], facecolor=axcolor)
s_zoom = Slider(axalpha, 'zoom: ', 0, 1, valinit=0.9)
axsm = plt.axes([0.25, 0.06,0.4, 0.04], facecolor=axcolor)
s_sm = Slider(axsm, 'solar masses: ', 1,157280, valinit=1)
radius = plt.axes([0.25, 0.11,0.4, 0.04], facecolor=axcolor)
s_sr = Slider(radius, 'radius: ', 0.1, 5, valinit=1)
rax = plt.axes([0.75, 0.01, 0.075, 0.14])
check = CheckButtons(rax, ('x=y', 'sun'), (True, True))
#Set constants note G is divided by c^2
global M
global G
global Ro
M=1.989e30
G=6.67408e-11/((3e8)**2)
Ro=696340000
def e1(sm,sr):
esqrt=np.sqrt((sr*Ro-sm*2*M*G)*(sr*Ro+sm*6*M*G))
e1=(sr*Ro-sm*2*M*G+esqrt)/(sm*4*M*G*sr*Ro)
return e1
def e2(sr):
e2=1/(sr*Ro)
return e2
def e3(sm,sr):
esqrt=np.sqrt((sr*Ro-sm*2*M*G)*(sr*Ro+sm*6*M*G))
e3=(sr*Ro-sm*2*M*G-esqrt)/(sm*4*M*G*sr*Ro)
return e3
def r(phi):
sm = s_sm.val
sr = s_sr.val
e1v=e1(sm,sr)
e2v=e2(sr)
e3v=e3(sm,sr)
tau=np.sqrt(sm*M*G*(e1v-e3v)/2)
h = (e2v-e3v)/(e1v-e3v)
sigma = -tau*np.pi/2 - sps.ellipk(h)
sn,_,_,_=sps.ellipj(tau*phi+sigma, h)
r=e3v+(e2v-e3v)*sn**2
return r
def alpha():
#Determine the root of solution
sm = s_sm.val
sr = s_sr.val
e1v=e1(sm,sr)
e2v=e2(sr)
e3v=e3(sm,sr)
h = (e2v-e3v)/(e1v-e3v)
phi=np.arcsin(np.sqrt(-e3v/(e2v-e3v)))
angle= sps.ellipkinc(phi, h)
tau=np.sqrt(sm*M*G*(e1v-e3v)/2)
sigma = -tau*np.pi/2 - sps.ellipk(h)
alphav = (angle-sigma) /tau-2*np.pi
return alphav
#The following function updateplot runs whenever slider changes value
def updateplot(val):
#Solar masses/radius and magnification alpha
zoom = s_zoom.val
sm = s_sm.val
sr = s_sr.val
#Find x range for all valid solutions
alphamax=alpha()
#Set angle range and calculate radius
phi= np.linspace(alphamax/2,np.pi-alphamax/2,500)
phi=phi[int(0.5*(1-zoom)*500):int(0.5*(1+zoom)*500)]
#phi= np.linspace(0,np.pi,500)
rv=1/r(phi)
clearsetplot()
#Check if results a biggen than photon sphere
if sr*Ro>sm*3*G*M:
x=rv*np.cos(phi)
y=rv*np.sin(phi)
#plot lightpath curve
ax1.plot(x,y)
ax1.grid()
#Angular Deflection
dx=np.diff(x)
dy=np.diff(y)
dydx=dy/dx
da=np.arctan(dydx)
dam=np.max(da)
ax2.plot(x[:-1],-da)
ax2.grid()
#Angular Distribution
ddy=np.diff(da)
ddydx=ddy/dx[:-1]
ax3.plot(x[:-2],-ddydx)
ax3.grid()
plotdatatext(sm,sr,dam)
#scaling method and draw sun
plotoptions(x,y,-da,-ddydx,sr,sm)
else:
#Plot Text No solution
plotnosolutiontext()
fig.canvas.draw_idle()
plt.plot
def plotdatatext(sm,sr,dam):
e1v=e1(sm,sr)
e2v=e2(sr)
e3v=e3(sm,sr)
tau=np.sqrt(sm*M*G*(e1v-e3v)/2)
h = (e2v-e3v)/(e1v-e3v)
text=('$r / r_s$: ' + str('{:.4e}'.format(Ro/(sm*2*G*M))) + '\n' +
str( '$r / r_p$: ' +'{:.4e}'.format(Ro/(sm*3*G*M))) )
ax1.text(0.98,0.98,text, horizontalalignment='right',verticalalignment='top', transform=ax1.transAxes, fontsize=10)
text=(str('{:.6f}'.format(dam)) + ' $[rad]$\n' +
str('{:.6f}'.format(np.degrees(dam))) + '$^{\circ}$\n' +
str('{:.6f}'.format(np.degrees(dam)*3600)) + "$''$\n")
ax2.text(0.02,0.98,text, horizontalalignment='left',verticalalignment='top', transform=ax2.transAxes, fontsize=10)
text=('$e1$: ' + str('{:.4e}'.format(e1v)) + '\n' +
str('$e2$: ' + '{:.4e}'.format(e2v)) + '\n' +
str('$e3$: ' + '{:.4e}'.format(e3v)) + '\n' +
str('$h^2$: ' + '{:.4e}'.format(h)) )
ax3.text(0.98,0.98,text, horizontalalignment='right',verticalalignment='top', transform=ax3.transAxes, fontsize=10)
def plotnosolutiontext():
text=('No solution')
ax2.annotate(text,xy=(0.5,0.5), xycoords='axes fraction',ha='center',color='red', fontsize=20)
#This function set scaling axis and drawes sun
def plotoptions(x,y,dy,ddy,sr,sm):
i=0
for r in check.get_status():
if i==0:
#Set scaling auto or equal
if r==False:
ax1.axis('auto')
#Set y scaling plot ax1
range=np.max(y)-np.min(y)
binsy=np.linspace(np.min(y)-range*2/100,np.max(y)+range*2/100, 5)
ax1.yaxis.set_ticks(binsy)
ax1.set_yticklabels(binsy, rotation=90, va="center")
ax1.yaxis.set_major_formatter(ticker.FormatStrFormatter('%.1e'))
#plt.setp(ax2.get_yticklabels()[-1], visible=False)
else:
ax1.axis('equal')
#Set y scaling plot ax1 (equal x=y and 5 bins)
a=ax1.get_yticks().tolist()
ax1.yaxis.set_major_locator(ticker.FixedLocator(a))
ax1.set_yticklabels(a,rotation=90, va="center")
ax1.yaxis.set_major_locator(ticker.MaxNLocator(5))
ax1.yaxis.set_major_formatter(ticker.FormatStrFormatter('%.1e'))
if i==1:
#Draw Sun yes/no
if r==True:
circle=plt.Circle((0,0),Ro, color='orange',label='$r_{sun}$: sun')
ax1.add_patch(circle)
circle=plt.Circle((0,0),sr*Ro,fill=False, ls='--')
ax1.add_patch(circle)
circle=plt.Circle((0,0),sm*3*G*M, color='gray',label='$r_p$: photon sphere')
ax1.add_patch(circle)
circle=plt.Circle((0,0),sm*2*G*M, color='black',label='$r_s$: Schwarzschild radius')
ax1.add_patch(circle)
ax1.legend(loc=3,fontsize='8',markerscale=0.25,facecolor='none',frameon=False,prop={'size': 8})
i=i+1
#Set x axis ticks
range=np.max(x)-np.min(x)
max=np.max(x)
binsx=np.linspace(-max-range*2/100,max+range*2/100, 5)
ax1.xaxis.set_ticks(binsx)
ax1.set_xticklabels(binsx, rotation=0, ha="left")
ax1.xaxis.set_major_formatter(ticker.FormatStrFormatter('%.1e'))
plt.setp(ax1.get_xticklabels()[-1], visible=False)
ax2.xaxis.set_ticks(binsx)
ax2.set_xticklabels(binsx, rotation=0, ha="left")
ax2.xaxis.set_major_formatter(ticker.FormatStrFormatter('%.1e'))
plt.setp(ax2.get_xticklabels()[-1], visible=False)
ax3.xaxis.set_ticks(binsx)
ax3.set_xticklabels(binsx, rotation=0, ha="left")
ax3.xaxis.set_major_formatter(ticker.FormatStrFormatter('%.1e'))
plt.setp(ax3.get_xticklabels()[-1], visible=False)
#Set y axis plots ax2 and ax3
range=np.max(dy)-np.min(dy)
binsy=np.linspace(np.min(dy)-range*2/100,np.max(dy)+range*2/100, 5)
ax2.yaxis.set_ticks(binsy)
ax2.set_yticklabels(binsy, rotation=90, va="center")
ax2.yaxis.set_major_formatter(ticker.FormatStrFormatter('%.1e'))
#plt.setp(ax2.get_yticklabels()[-1], visible=False)
range=np.max(ddy)-np.min(ddy)
binsy=np.linspace(np.min(ddy)-range*2/100,np.max(ddy)+range*2/100, 5)
ax3.yaxis.set_ticks(binsy)
ax3.set_yticklabels(binsy, rotation=90, va="center")
ax3.yaxis.set_major_formatter(ticker.FormatStrFormatter('%.1e'))
#plt.setp(ax3.get_yticklabels()[-1], visible=False)
def clearsetplot():
#plot lightpath curve
ax1.clear()
ax1.set_title('Light path', fontsize=12, pad=10)
ax1.set_xlabel('x [m]', fontsize=12)
ax1.set_ylabel('y [m]', fontsize=12)
ax1.tick_params(axis='both', which='major', labelsize=8)
ax2.clear()
ax2.set_title('Angular Deflection', fontsize=12, pad=10)
ax2.set_xlabel('x [m]', fontsize=12)
ax2.set_ylabel('deflection angle [rad]', fontsize=12)
ax2.tick_params(axis='both', which='major', labelsize=8)
ax3.clear()
ax3.set_title('Angular Distribution', fontsize=12, pad=10)
ax3.set_xlabel('x [m]', fontsize=12)
ax3.set_ylabel('angular change [rad/m]', fontsize=12)
ax3.tick_params(axis='both', which='major', labelsize=8)
#Display plot n startup.
updateplot(1)
#The following code checkes if slider changes. This line is looped automatic by pyplot
s_sm.on_changed(updateplot);
s_sr.on_changed(updateplot);
s_zoom.on_changed(updateplot);
check.on_clicked(updateplot);
Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …
from IPython.display import YouTubeVideo
YouTubeVideo('0xngLoNqM1g')