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)
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);
In [11]:
from IPython.display import YouTubeVideo
YouTubeVideo('0xngLoNqM1g')
Out[11]:
In [ ]: