Re: Twee pieken of toch maar één?
Geplaatst: do 22 jul 2021, 16:15
Onder voorbehoud (hier zitten wellicht nog fouten in):
Graag commentaar van ervaren Python gebruikers.
Code: Selecteer alles
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import ellipj
from scipy.special import ellipk
fig, ax= plt.subplots(4,1,figsize=(4, 6))
# Berekening lichtbaan
c = 3e8
r0 = 7e8
rs = 2.95e3
e1 = (r0 - rs + np.sqrt((r0 - rs)*(r0 + 3*rs)))/(2*rs*r0)
e2 = 1/(r0)
e3 = (r0 - rs - np.sqrt((r0 - rs)*(r0 + 3*rs)))/(2*rs*r0)
tau = np.sqrt((rs*(e1 - e3))/4)
h = np.sqrt((e2 - e3)/(e1 - e3))
m = h*h
sigma = -tau*(np.pi/2) + ellipk(m)
t = np.linspace(np.pi/2-1.4,np.pi/2+1.4,1000)
sn,_,_,_=ellipj(tau*t + sigma , m)
x = (1/(e3 + (e2 - e3)*sn**2))*np.cos(t)
y = (1/(e3 + (e2 - e3)*sn**2))*np.sin(t)
# Berekening diverse afgeleiden
dx=np.diff(x)
dy=np.diff(y)
dydx=dy/dx
r=np.sqrt(x**2 + y**2)
S=(x[:-1] + y[:-1]*dydx)/r[:-1]
T=(S*(x[:-1]/r[:-1]) - 1)/y[:-1]
dx2dt2=(( (1 - (rs/r[:-1]))**2 )/( S**2 + (r[:-1]**2)*(1 - (rs/r[:-1]))*(T**2) ))*(c**2)
ddydx=np.diff(dydx)
ddydxdx=ddydx/dx[:-1]
# Plots
ax[0].plot(x,y,c=(1,0.2,0.5),lw=1)
ax[0].title.set_text('y = f(x)')
ax[1].plot(x[:-1],dydx,c=(1,0.2,0.5),lw=1)
ax[1].title.set_text('dy/dx')
ax[2].plot(x[:-1],dx2dt2,c=(1,0.2,0.5),lw=1)
ax[2].title.set_text('dx2/dt2')
ax[3].plot(x[:-2],ddydxdx,c=(1,0.2,0.5),lw=1)
ax[3].title.set_text('d2y/dx2')
plt.tight_layout()
plt.show()