Python en Jacobi

Moderators: jkien, Xilvo

Reageer
Gebruikersavatar
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.

Gebruikersavatar
Berichten: 1.605

Re: Python en Jacobi

Ik zie nu waar ik verkeerd was met Tau.

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)
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.

Gebruikersavatar
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.

Gebruikersavatar
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.151

Re: Python en Jacobi

Vincent
Tip:

Code: Selecteer alles

M=1.989*10**30
G=6.67408*10**(-11)
Ro=7*10**8
Is hetzelfde als

Code: Selecteer alles

M=1.989e30
G=6.67408e-11
Ro=7e8
Maar dat laatste schrijft wel veel makkelijker en overzichtelijker.

Gebruikersavatar
Berichten: 1.605

Re: Python en Jacobi

Dankjewel! Ik werk nauwelijks/niet met zulke getallen in te voeren!

Gebruikersavatar
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) \)
\(\)

Gebruikersavatar
Berichten: 7.463

Re: Python en Jacobi

sigma.png
Dus:
\(\)
\( \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...

Gebruikersavatar
Berichten: 1.605

Re: Python en Jacobi

Orbits.png
Ondertussen ben ik de eigenschappen aan het verkennen van de functie. Ikzelf vind het makkelijker eerst alles te normaliseren en niet allemaal decimale getallen in te vullen. Ik werk liever met uitdrukkingen in pi.

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?

Gebruikersavatar
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:
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()
probeersel.png
probeersel.png (14.29 KiB) 596 keer bekeken
Mogelijk zitten hier nog fouten in, maar dit geeft mij wel voldoende hoop op de goede weg te zijn.

Gebruikersavatar
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.

Gebruikersavatar
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. ;)

Gebruikersavatar
Berichten: 7.463

Re: Python en Jacobi

jacobi.png

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(φ).

Gebruikersavatar
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?

Gebruikersavatar
Berichten: 4.518

Re: Python en Jacobi

Mapleplot op basis van dit fragment uit een eerder gepost paper komt aardig overeen..
Lightdeflection by gravity trajectory1.png
Lightdeflection by gravity trajectory.png

Reageer