Python en Jacobi
- Berichten: 7.463
Re: Python en Jacobi
Ik wacht even af tot je de tau gecorrigeerd hebt en dan probeer ik het nog een keer. Ik kreeg wel een grafiek maar die zag er inderdaad raar uit.
- Berichten: 1.605
Re: Python en Jacobi
Ik zie nu waar ik verkeerd was met Tau.
Enkele uren geleden heb ik geprobeerd met genormaliseerde getallen. Dat lukt niet omdat dan de wortel negatief word. Bij enkele willekeurige instellingen krijg ik cirkels of concentrische cirkels.
Een hyperbool als afbuiging heb ik (nog) niet gezien.
Moet men eigenlijk voor \(r_{o}\) overal \(r\) invullen en dat oplossen voor \(r\)? Net zoals die in andere topics?
Voor formule (6) in jouw samenvatting schrijft in document "The equation of the orbit is thus written as:" en heeft het over orbit en niet over lightpath.
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('e3: ' + str(e3))
Tau=np.sqrt(M*G*(e1-e3)/(2*c**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)
Een hyperbool als afbuiging heb ik (nog) niet gezien.
Moet men eigenlijk voor \(r_{o}\) overal \(r\) invullen en dat oplossen voor \(r\)? Net zoals die in andere topics?
Voor formule (6) in jouw samenvatting schrijft in document "The equation of the orbit is thus written as:" en heeft het over orbit en niet over lightpath.
- Berichten: 7.463
Re: Python en Jacobi
Wat het euvel zou kunnen zijn is dat we σ nog niet netjes bepaald hebben. Eigenlijk moet je die constante oplossen door een randvoorwaarde te kiezen. Dat bepaalt dan de specifieke baan die je krijgt. Bijvoorbeeld zouden we kunnen kiezen dat r minimaal (d.w.z. r0) is voor φ = π/2. Dat wordt immers meestal zo gekozen.
- Berichten: 7.463
Re: Python en Jacobi
Ja - zo moet dat kunnen. Je kunt r als functie van σ schrijven door φ gelijk aan π/2 te stellen, en dan is de specifieke σ die de door ons gezochte lichtbaan oplevert zo'n σ waarvoor r gelijk aan r0 wordt. Als er meerdere van zulke σ's zijn, moeten we nog bekijken voor welke σ we bij φ = π/2 ook inderdaad een minimum van r hebben.
-
- Technicus
- Berichten: 1.162
Re: Python en Jacobi
Vincent
Tip:
Is hetzelfde als
Maar dat laatste schrijft wel veel makkelijker en overzichtelijker.
Tip:
Code: Selecteer alles
M=1.989*10**30
G=6.67408*10**(-11)
Ro=7*10**8
Code: Selecteer alles
M=1.989e30
G=6.67408e-11
Ro=7e8
- Berichten: 1.605
Re: Python en Jacobi
Dankjewel! Ik werk nauwelijks/niet met zulke getallen in te voeren!
- Berichten: 7.463
Re: Python en Jacobi
Even terug naar:
\(\)
\( r = \frac{1}{ e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \varphi + \sigma ; h)} \,\,\,\,\,\,\, (6) \)
\(\)
Nu bekijken we de lichtbaan waarvoor r minimaal (d.w.z. r0) is voor φ = π/2. Dat geeft:
\(\)
\( r_0 = \frac{1}{ e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \frac{\pi}{2} + \sigma ; h)} \)
\(\)
En wegens (2) dat:
\(\)
\( \frac{1}{e_2} = \frac{1}{ e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \frac{\pi}{2} + \sigma ; h)} \)
\(\)
\( e_2 = e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \frac{\pi}{2} + \sigma ; h) \)
\(\)
\( e_2 - e_3 = (e_2 - e_3) \mathrm{sn}^2( \tau \frac{\pi}{2} + \sigma ; h) \)
\(\)
\( \mathrm{sn}^2( \tau \frac{\pi}{2} + \sigma ; h) = 1 \,\,\,\,\,\,\, (9) \)
\(\)
- Berichten: 7.463
Re: Python en Jacobi
\(\)
\( \tau \frac{\pi}{2} + \sigma = \pm 1,57194 \)
\(\)
\( \sigma = -\tau \frac{\pi}{2} \pm 1,57194 \)
\(\)
En wegens (8):
\(\)
\( \sigma = -0,500001 \cdot \frac{\pi}{2} \pm 1,57194 \)
\(\)
\( \sigma = -0,78540 \pm 1,57194 \)
\(\)
Nu nog zien welke daarvan voldoet...- Berichten: 1.605
Re: Python en Jacobi
Als ik de functie verken krijg ik plot uit zoals toegevoegde. Het convergeert naar een cirkel.
Ik vermoed dat de formules van \((9.2)\) "Trajectory of Light Signal using Jacobi’s Elliptic Function" voor een omloopbaan zijn van licht. Bij de laatste formule staat ook: "The equation of the orbit is thus written as":
$$r = \frac{1}{ e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \varphi + \sigma ; h)} \,\,\,\,\,\,\, (6)$$
Voor deflectie moeten we volgens mij hoofdstuk \((6)\) gebruiken: "Deflection of Light in the Sun’s Gravitational Field.":
$$ \Delta \varphi= 4 \left[ \frac{ r_{0} }{ r_{0}-2MG } \sqrt{ \frac{r_{0}-2MG}{r_{0}+6MG }} \right ]^{1/4} F(a;k) - \pi$$
Waarbij \(F(a;k)\) de elliptische integraal is first kind.
De Jacobi’s Elliptic Function \(sn\) is periodic. Dus volgens mij gaat het om omloopbanen volgens mij. Vandaar dat ik denk dat \(9.2\) niet de correcte is voor deflektie. Tenzij je bewust de doelstelling veranderd hebt van deflektie naar omloopbaan?
- Berichten: 7.463
Re: Python en Jacobi
Volgens mij gaat hoofdstuk 9 van het artikel niet (specifiek) over gesloten lichtbanen, maar hoofdstuk 10 wel. En bij 10 wordt dat ook vermeld:
Hier een Python probeersel uitgaande van formule (6) met de gevonden constanten:
Mogelijk zitten hier nog fouten in, maar dit geeft mij wel voldoende hoop op de goede weg te zijn.
10. Period of Revolution
We shall now calculate the period of revolution of a test body describing a closed orbit around a heavy object of mass M producing the gravitational field.
Hier een Python probeersel uitgaande van formule (6) met de gevonden constanten:
Code: Selecteer alles
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import ellipj
t = np.linspace(np.pi/2 - 1,np.pi/2 + 1,1000)
sn,_,_,_=ellipj(0.500001*t+(-0.78540 + 1.57194),0.00290319)
x = (1/(-1.42857e-9 + ( (7e8)**-1 - -1.42857e-9)*sn*sn))*np.cos(t)
y = (1/(-1.42857e-9 + ( (7e8)**-1 - -1.42857e-9)*sn*sn))*np.sin(t)
plt.plot(x,y,c=(1,0.2,0.5),lw=1)
plt.title('lichtbaan')
plt.axis('on')
plt.show()
- Berichten: 1.605
Re: Python en Jacobi
Dat ziet er mooi uit! Mooi dat je ook al de grafieken begint te pimpen met kleur en lijnbreedte.
Ik zit achter telefoon. Volgen mij kan je kleuren ook aangeven met tekst: color="red". Maar kan excacte code nu niet controleren.
Erg interessant die Jacobi elliptics. Niet dat ik het allemaal begrijp maar dan leert men de taal een beetje spreken en herkennen.
Ik zit achter telefoon. Volgen mij kan je kleuren ook aangeven met tekst: color="red". Maar kan excacte code nu niet controleren.
Erg interessant die Jacobi elliptics. Niet dat ik het allemaal begrijp maar dan leert men de taal een beetje spreken en herkennen.
- Berichten: 7.463
Re: Python en Jacobi
Ja - ben nu nog even aan het natrekken of die twee sigma's exact dezelfde lichtbaan opleveren.
Maar als dit resultaat klopt kunnen we numeriek gaan differentiëren om te zien hoeveel toppen er nu echt - bij zo exact mogelijke berekening - verschijnen. Dat is immers de reden dat ik hier überhaupt aan begonnen ben.
Mijn programmaatje is een soort van monster van Frankenstein. Ik heb uit verschillende voorbeeld programmaatjes geknipt en geplakt en dat stapsgewijs aan ons vraagstuk aangepast.
Maar als dit resultaat klopt kunnen we numeriek gaan differentiëren om te zien hoeveel toppen er nu echt - bij zo exact mogelijke berekening - verschijnen. Dat is immers de reden dat ik hier überhaupt aan begonnen ben.
Mijn programmaatje is een soort van monster van Frankenstein. Ik heb uit verschillende voorbeeld programmaatjes geknipt en geplakt en dat stapsgewijs aan ons vraagstuk aangepast.
- Berichten: 7.463
Re: Python en Jacobi
Bron: https://mathworld.wolfram.com/JacobiEll ... tions.html
Hieruit leren we dat:
\(\)
\( \mathrm{sn}(u + 2 m K + 2 n i K' , k ) = (-1)^m \mathrm{sn}(u,k) \)
\(\)
\( (\mathrm{sn}(u + 2 m K + 2 n i K' , k ))^2 = ((-1)^m \mathrm{sn}(u,k))^2 \)
\(\)
\( \mathrm{sn}^2(u + 2 m K + 2 n i K' , k ) = \mathrm{sn}^2(u,k) \)
\(\)
Dus voor reële argumenten u heeft \( \mathrm{sn}^2(u,k) \) een periode 2K. Dit zelfde geldt dan wegens (6) ook voor r(φ).- Berichten: 1.605
Re: Python en Jacobi
Ik ben jouw herleiding nagegaan voor \(\sigma\) en snap jouw gedachten gang. Tevens kan ik jouw grafiek ook reproduceren.
Echter als ik de deflektie uitreken kom ik vele orden groter dan verwacht. Uit jouw grafiek komt ik uit op: 0.0014 rad (halve hoek). Men zou moeten verwachten 0.9 arc seconde. Wat bereken jij voor een deflectie?
Echter als ik de deflektie uitreken kom ik vele orden groter dan verwacht. Uit jouw grafiek komt ik uit op: 0.0014 rad (halve hoek). Men zou moeten verwachten 0.9 arc seconde. Wat bereken jij voor een deflectie?
- Berichten: 4.540
Re: Python en Jacobi
Mapleplot op basis van dit fragment uit een eerder gepost paper komt aardig overeen..