Python en Jacobi

[Professor Puntje]

Er zijn wel grenzen aan de zinnigheid van het gelijk 1 stellen. Als je bijvoorbeeld stelt dat r_s = r₀ = 1 dan stel je dat onze zon een zwart gat is. En vervolgens bekijk je wat er met licht gebeurt dat dat zwarte gat tot op de waarnemingshorizon nadert. Dat is niet waar we hier in geïnteresseerd zijn...

[OOOVincentOOO]

Volgens mij begrjip je mij niet. Om een programma of formule op correctheid te controleren dient je niet alle details vooraf in te vullen. Immers uiteindelijk [Wiki]:

$$\theta=\frac{4GM}{rc^2}$$

Het belang is eerst te controleren of de factor 4 aanwezig is.

Jij hebt blijkbaar je eigen aanpak en je begint wederom met de Schwarschild radius. Dan maak je het zelf volgens mij alleen maar onnodig ingewikkeld. Ik spreek alleen uit ervaring hoe ik dingen mijzelf en mij geleerd is door anderen. En ik kan verkeerd zijn.

[Professor Puntje]

Als je de formules eerst wilt controleren wordt het inderdaad een ander verhaal. Zie je een eenvoudige manier om dat te doen?

[OOOVincentOOO]

Ik probeer eerst jouw mooie samenvatting van functies te begrijpen.

$$r = \frac{1}{ e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \varphi + \sigma ; h)} \,\,\,\,\,\,\, (6)$$

Volgens mij is \(\sigma\) ook een constante? Of ben ik verkeerd?

Misschien kunnen we deze functie plotten als \(y=r \cdot sin(\varphi)\) en \(x=r \cdot cos(\varphi)\) (indien dat mag) en kijken of er uberhaubt een ellipsachtige orbit en/of hyperbool achtige curve's uitkomen?

Ik schrijf curve's meervoud omdat ik vermoed dat de functie periodiek is. Maar dat weet ik niet zeker. Soms is speculeren leuk!

[Professor Puntje]

Over σ zegt het artikel dit:

We may calculate the constant σ in the above formula by imposing some initial
conditions.

En sn is voor complexe argumenten zelfs dubbelperiodiek.

Voor de asymptoten van de lichtbaan moet r naar oneindig gaan, en dat gebeurt als de noemer naar nul gaat.

Ook moeten we de baan voor het toppen-onderzoek uiteindelijk in een xy-stelsel beschrijven.

[Professor Puntje]

Ik ga onderwijl toch maar even verder met de specifieke waarden van de diverse constanten voor een lichtstraal die rakelings langs de zon scheert. Ik ben namelijk heel benieuwd wat daar uit komt.

Code: Selecteer alles

 calculate (r_0 - r_1 + sqrt((r_0 - r_1)*(r_0 + 3*r_1)))/(2*r_1*r_0) for r_0 = 7*10^8  and r_1= 2.95*10^3 

(Met r₁ in plaats van r_s omdat WolframAlpha voor r_s moeilijk doet.

[Professor Puntje]

Code: Selecteer alles

 1/(7*10^8) 

[Professor Puntje]

Code: Selecteer alles

 calculate (r_0 - r_1 - sqrt((r_0 - r_1)*(r_0 + 3*r_1)))/(2*r_1*r_0) for r_0 = 7*10^8  and r_1= 2.95*10^3 

[Professor Puntje]

Code: Selecteer alles

 sqrt(( r_1*(e_1 - e_3)  )/(4)) for r_1 = 2.95*10^3 & e_1 = 0.000338983 & e_3 = -1.42857*10^(-9) 

[Professor Puntje]

Code: Selecteer alles

 sqrt(( e_2 - e_3)/(e_1 - e_3)) for e_1 = 0.000338983 & e_2 = 1/(7*10^8)    & e_3 = -1.42857*10^(-9) 

[Professor Puntje]

We hebben voor een lichtstraal die rakelings langs de zon scheert dus nu gevonden:

\(\)

\( \left. \begin{array} {lcrr} e_1 = 0,000338983 \, \mathrm{m}^{-1} \\ e_2 = (7 . 10^8)^{-1} \, \mathrm{m}^{-1} \\ e_3 = -1,42857 . 10^{-9} \, \mathrm{m}^{-1} \\ \tau = 0,500001 \\ h = 0,00290319 \end{array} \right \} \,\,\,\,\,\, (8) \)

[OOOVincentOOO]

Kijk eens of jij ermee verder komt. Heb het geprobeerd overzichtelijk te houden zonder veel toeters en bellen. Ik krijg een onzin grafiek. Maar ben nog niet vertrouwd met de wiskunde.

Code: Selecteer alles

import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sps
fig, ax1= plt.subplots(figsize=(15, 15))

M=1.989*10**30
G=6.67408*10**(-11)
Ro=7*10**8
Sigma=0
c=300000000

phi=np.linspace(0,20*np.pi, 100000)

esqrt=np.sqrt((Ro-2*M*G/c**2)*(Ro+6*M*G/c**2))

e1=(Ro-2*M*G/c**2+esqrt)/(4*M*G*Ro/c**2)
print('e1: ' + str(e1))

e2=1/Ro
print('e2: ' + str(e2))

e3=(Ro-2*M*G/c**2-esqrt)/(4*M*G*Ro/c**2)
print('e2: ' + str(e3))

Tau=np.sqrt(M*G/(c**2*2))
print('Tau: ' + str(Tau))

h=np.sqrt((e2-e3)/(e1-e3))
print('h: ' + str(h))

val=Tau*phi+Sigma
sn,_,_,_=sps.ellipj(val,h)
r=1/(e3+(e2-e3)*sn**2)

x=r*np.cos(phi)
y=r*np.sin(phi)

ax1.plot(x,y)

[Professor Puntje]

Als je niet met r_s werkt moet je net als de auteur van het artikel c=1 nemen. Die op 1 gestelde constanten zijn lastig om mee te werken, tenzij je dat gewend bent.

[OOOVincentOOO]

Ik heb het ingesteld zodat je: e1, e2, e3, Tau en h jouw waarden zijn. Dat kun je zien als je de code test! Ik print speciaal voor jouw de waarden erin gezet. Dus ik begrijp jouw opmerking niet.

Code: Selecteer alles

e1: 0.0003389895594409787
e2: 1.4285714285714286e-09
e2: -1.4285654083149677e-09
Tau: 27.156690520017346
h: 0.0029031631580164635

Wacht even ik zie dat Tau nog een foutje heeft.

[Professor Puntje]

Ben je van formule (6) uitgegaan?

Ik lees ook:

Code: Selecteer alles

M=1.989*10**30
G=6.67408*10**(-11)
Ro=7*10**8
Sigma=0
c=300000000

Wat doe je daarmee? Alleen sigma komt daarvan in de formule (6) terug.