Python en Jacobi

[OOOVincentOOO]

Dat heb ik gedaan in mijn simulatie.

Om de deflectie hoek te krijgen moet de x-as 320 maal vergroot worden.

Ik ben niet overtuigd dat dit een decimalen probleem is.

[Professor Puntje]

Geeft jouw simulatie nu ondanks de preciezere berekening precies dezelfde lichtbaan als de mijne?

[Professor Puntje]

We gebruiken allebei de functie "ellipj".
Zie: https://docs.scipy.org/doc/scipy/refere ... llipj.html

Deze functie heeft de parameter "m". Maar welke parameter is dat?
Volgens Wikipedia zijn er drie opties: de parameter m, de modulus k en de modulaire hoek α.
Zie: https://en.wikipedia.org/wiki/Jacobi_el ... s#Notation

En welke parameter gebruikt Solomon Antoniou? Hij schrijft h in plaats van k of m. Op zijn formules is dit topic gebaseerd. Ik ben er onwetend van de drie mogelijke notaties vanuit gegaan dat zowel WolframAlpha, als Solomon Antoniou, als ellipj de modulus k gebruiken, maar dat hoeft niet zo te zijn. Ook dit is dus nog een mogelijke bron van fouten.

[CoenCo]

De scipy module geeft aan dat het een wrapper is voor cephes. Verdere documentatie staat dus hier: http://www.netlib.org/cephes/doubldoc.html#ellpj

En zie ook hier:
https://docs.scipy.org/doc/scipy/refere ... ial.ellipk

[OOOVincentOOO]

Hallo pp,

Ik heb nog terug gebladerd naar hoofdstuk 6. Hier staat de elliptic integral (zie ook: [Wiki]). Hier staat \(k^2\).

Nu heb ik een beetje gestoeid. Ik heb in de formule \(sn(\varphi;h^2)\) heb ik \(h^{2}\) ingevuld. De constante krijg ik nog niet goed berekend. Het figuur is net niet symmetrisch. Maar orde grote krijg ik volgens mij wel de juiste deflektie hoek.

Code: Selecteer alles

Wolfram Alpha:
solve (sn(x,0.0029031631580164635^2)^2=1

\(x=-1.5708\)

Dit levert net geen symmetrisch figuur. Na wat tweaking wel symmetrisch:

\(x=-1.5708013\)

Zodat

\(\sigma=-\pi/4-1.5708013\)

Dan kom ik op de deflectie van een zijde uit op:

\(\Delta \varphi \approx \frac{7}{1.75} \cdot 10^{-6}=4 \cdot 10^{-6} \ [rad]\)
\(\Delta \varphi \approx 0.000229 \ [^{\circ}]\)
\(\Delta \varphi \approx 0.000229*3600=0.825059 \ [arcsec]\)

Maar dit moet volgens mij eigenlijk 0.45 arcsec zijn, ik weet de definitie van de deflektiehoek niet. Volgens mij is de totale deflektie 0.9 arcsec.

Gaarne mijn methode en berekening controleren. Jij bent beter in het speur en verifieer werk. Denk jij dat dit mogelijk is?

Code: Selecteer alles

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

M=1.989e30
G=6.67408e-11/((3e8)**2)
Ro=7e8
#Sigma=-np.pi/4+1.57194
Sigma=-np.pi/4-1.5708013

print('sigma: ' + str(Sigma))

phi=np.linspace(0,np.pi,1000)

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

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

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

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

Tau=np.sqrt(M*G*(e1-e3)/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**2)
r=1/(e3+(e2-e3)*(sn**2))


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

ax1.grid()
ax1.plot(x,y)

[OOOVincentOOO]

In Wiki staat meer over notaties k en k² [Wiki]:

Van Wiki:
"Notational variants: There are still other conventions for the notation of elliptic integrals employed in the literature. The notation with interchanged arguments, F(k,φ), is often encountered; and similarly E(k,φ) for the integral of the second kind. Abramowitz and Stegun substitute the integral of the first kind, F(φ,k), for the argument φ in their definition of the integrals of the second and third kinds, unless this argument is followed by a vertical bar: i.e. E(F(φ,k) | k2) for E(φ | k2). Moreover, their complete integrals employ the parameter k² as argument in place of the modulus k, i.e. K(k²) rather than K(k). And the integral of the third kind defined by Gradshteyn and Ryzhik, Π(φ,n,k), puts the amplitude φ first and not the "characteristic" n."

[Professor Puntje]

Ziet er interessant uit! Ga even bekijken of dat werkt.

[Professor Puntje]

Weer even terug naar:

\(\)

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

\(\)

Het is goed mogelijk dat we h² in plaats van h moeten nemen. Immers heeft men het bij ellipj over de parameter m in plaats van k. Daarom bepaal ik nu σ uitgaande van de aangepaste formule (6*):

\(\)

\( r = \frac{1}{ e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \varphi + \sigma ; h^2)} \,\,\,\,\,\,\, (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^2)} \)

\(\)

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^2)} \)

\(\)

\( e_2 = e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \frac{\pi}{2} + \sigma ; h^2) \)

\(\)

\( e_2 - e_3 = (e_2 - e_3) \mathrm{sn}^2( \tau \frac{\pi}{2} + \sigma ; h^2) \)

\(\)

\( \mathrm{sn}^2( \tau \frac{\pi}{2} + \sigma ; h^2) = 1 \)

\(\)

Aangezien sn ook voor reële argumenten periodiek is zijn daar oneindig veel oplossingen voor. We nemen een simpel geval en bekijken:

\(\)

\( \mathrm{sn}( \tau \frac{\pi}{2} + \sigma ; h^2) = 1 \)

\(\)

Proberenderwijs met Spyder en ellipj kom ik op:

\(\)

\( \sigma = -\tau \frac{\pi}{2} + 1,57079965 \)

\(\)

[Professor Puntje]

Code: Selecteer alles

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import ellipj

alpha = 8.1e-6

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)
h2 = (e2 - e3)/(e1 - e3)
sigma = -tau*(np.pi/2) + 1.57079965

t = np.linspace(-alpha/2,np.pi+alpha/2,1000)

sn,_,_,_=ellipj(tau*t + sigma , h2)

x = (1/(e3 + (e2 - e3)*sn*sn))*np.cos(t)
y = (1/(e3 + (e2 - e3)*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()

Je kunt hiermee de hoek alpha stapsgewijze vergroten totdat je aan een waarde komt die groter dan de totale lichtafbuiging is, en de grafiek onzinnig wordt. Ook kan zo waarschijnlijk sigma nog iets preciezer bepaald worden.

[OOOVincentOOO]

Inderdaad ik heb geplot van: 0 tot pi. Linker grafiek om hoek te bepalen. Jouw sigma is netter bepaald. Ik had Tau op 0.5 staan. Volgens mij kom jij op vergelijkbare resultaten dan ik heb net jouw code laten lopen en uitgerekend.

Het document refereert naar (Ref: jouw orginele post):

Op wiki staat integraal vorm [Wiki]:

Beide hebben het over \(h^2/k^{2}\). Kan jij inschatten of we op de correcte weg zitten, orde grote zitten we goed (denk ik)? Ik snap alleen niet waar de volgende kleine fout zit, de deflektie zou: 0.9 arcsec dienen te zijn. Maar genoeg voor mij morgen weer werken.

Als we alles correct werkend hebben kunnen we een slidebar invoeren voor jouw alpha. Dan kan je interactief inzoomen.

[Professor Puntje]

Code: Selecteer alles

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import ellipj

alpha = 8.42e-6

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)
h2 = (e2 - e3)/(e1 - e3)
sigma = -tau*(np.pi/2) + 1.570799636678

t = np.linspace(-alpha/2,np.pi+alpha/2,1000)

sn,_,_,_=ellipj(tau*t + sigma , h2)

x = (1/(e3 + (e2 - e3)*sn*sn))*np.cos(t)
y = (1/(e3 + (e2 - e3)*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()

Draai dit eens. Zie ook dit voor de afbuiging in radialen: https://physics.stackexchange.com/quest ... by-the-sun

[OOOVincentOOO]

Oke, dankjewel. Ik bekijk morgen de code. Je heb blijkbaar weer een of andere (mooie) black magic toegepast. Dat kost meestal tijd bij mij om te verwerken.

Morgen een keer [Stacks] standaard deflektie formule narekenen. Als ik google "deflection of light by sun" krijg ik steeds: 0.9 arcsec als resultaat. De integraal schijnt 1.75 arcsec. te zijn volgens jouw link, en komt blijkbaar meer overeen met onze resultaten.

[Professor Puntje]

Inderdaad - morgen verder. Maar ik kan mij niet voorstellen dat we nu nog fout zitten. De door ons gevonden waarde voor de totale afbuiging komt nu vrijwel exact overeen met de alom bekende waarde. Maak je in mijn progje alpha nóg groter dan krijg je een onzinnige grafiek, wat ook logisch is omdat je dan r(φ) voor hoeken φ berekent waar het licht helemaal niet komt.

[OOOVincentOOO]

@CoenCo,

Zojuist jouw links er nog op nageslagen [Scipy]. Hier lijkt de factor \(m\) in elleptic integral inderdaad zonder kwadraat.

Een vraag als je tijd en zin hebt. Via Jupyter notebook kan ik relatief eenvoudig een pyplot grafiek maken interactief met widget's. En kan een slidebar toevoegen waarna de grafiek automatisch geupdate wordt. Ook kan ik via Jupyter notebook makkelijk coordinaten uitlezen en inzoomen in grafieken.

Echter via Python spyder console krijg ik dat niet makkelijk voor elkaar (heb vorige week te vergeefs geprobeerd). Heb jij enige ideeen op een heel simpele manier een Pyplot grafiek interactief te maken in combinatie met een widget?

Dan kunnen we misschien de grafiek inzoomen door 'alpha' te veranderen met een slidebar. He zou tevens handig zijn als de Pyplot interactief is zodat we posities interactief kunnen uitlezen in grafiek.

[Professor Puntje]

Verder is het nog interessant om y numeriek naar x te differentiëren wat tan(β(x)) oplevert, waarin β(x) de hellingshoek van de grafiek is. Hoe moet dat?