Python en Jacobi

[Professor Puntje]

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.

[OOOVincentOOO]

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.

[Professor Puntje]

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. r₀) is voor φ = π/2. Dat wordt immers meestal zo gekozen.

[Professor Puntje]

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 r₀ wordt. Als er meerdere van zulke σ's zijn, moeten we nog bekijken voor welke σ we bij φ = π/2 ook inderdaad een minimum van r hebben.

[CoenCo]

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.

[OOOVincentOOO]

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

[Professor Puntje]

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. r₀) 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) \)

\(\)

[Professor Puntje]

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

[OOOVincentOOO]

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?

[Professor Puntje]

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 890 keer bekeken

Mogelijk zitten hier nog fouten in, maar dit geeft mij wel voldoende hoop op de goede weg te zijn.

[OOOVincentOOO]

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.

[Professor Puntje]

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.

[Professor Puntje]

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

[OOOVincentOOO]

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?

[ukster]

Mapleplot op basis van dit fragment uit een eerder gepost paper komt aardig overeen..