sciencetalk

Code: Selecteer alles

%1.75 arcsec = 0.000008484239419416949 rad

rSchwarz=2950; %Schwarzild straall zon
rAfstandZon=696340000; %straal zon

fun = @(r) 1./((r.^2).*(1./(rAfstandZon^2)-(1-rSchwarz./r)./r.^2 ).^0.5);

phi = 2*integral(fun, rAfstandZon,Inf)

Code: Selecteer alles

#Open pyplot in separate interactive window
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'qt5')

#https://www.mathpages.com/rr/s8-09/8-09.htm
#https://www.mathpages.com/rr/s6-03/6-03.htm

import numpy as np
from scipy.special import comb
import matplotlib.pyplot as plt

widths = [5,5,5,5]
heights = [5]

fig= plt.figure(figsize=(20,5))

gs=fig.add_gridspec(1,4,width_ratios=widths, height_ratios=heights)

ax1=fig.add_subplot(gs[0,0])
ax2=fig.add_subplot(gs[0,1])
ax3=fig.add_subplot(gs[0,2])
ax4=fig.add_subplot(gs[0,3])


#Physical constants
M=1.989e30
G=6.67408e-11
c=3e8
Ro=696340000

#Schwarzshield radius
Rs=2*M*G/c**2

#set radius q=Ro/R range
q=np.linspace(0.005,1,1000)
q=q[1:-1]

#calculate radius (reciprocal)
r=Ro/q

#leading term power series and element b**n
a=1/np.sqrt(1-q**2)
b=2*(1-q**3)/(1-q**2)*M*G/(Ro*c**2)
dp=np.zeros(np.size(q))

dq=q[10]-q[9]

#Loop over every element power series
for n in range(5):

    i=n+1
    
    #determine fractions
    #https://oeis.org/A046161
    #(2n choose n)/2^n
    com=comb(2*i, i)/4**i
    dp=dp+a*com*b**i

    #plot angular change as function q and r note: dq=-(q^2/Ro)dr
    ax1.plot(q,dp)
    ax1.set_ylabel('dphi/dq')
    ax1.set_xlabel('q')
    
    ax3.plot(r,dp*(q**2/Ro))
    ax3.set_ylabel('dphi/dr')
    ax3.set_xlabel('r')

    #Integrate angular deflection     
    phi=np.cumsum(dp)*dq
    phisec=np.degrees(phi)*3600
    phisec=phisec[-1]-phisec
       
    ax2.plot(q,phisec)
    ax2.set_ylabel('dphi/dq')
    ax2.set_xlabel('q')
    
    ax4.plot(r,phisec)
    ax4.set_ylabel('dphi/dr')
    ax4.set_xlabel('r')