## Trajectory of Light Signal using Jacobi’s Elliptic Function.¶

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}$$
In [11]:
#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)

# 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)

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):

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')

circle=plt.Circle((0,0),sr*Ro,fill=False, ls='--')

circle=plt.Circle((0,0),sm*3*G*M, color='gray',label='$r_p$: photon sphere')

circle=plt.Circle((0,0),sm*2*G*M, color='black',label='$r_s$: Schwarzschild radius')

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_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_xlabel('x [m]', fontsize=12)
ax2.tick_params(axis='both', which='major', labelsize=8)

ax3.clear()
ax3.set_xlabel('x [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);

In [11]:
from IPython.display import YouTubeVideo