Python en Jacobi

[OOOVincentOOO]

Interessant met r0 < rs = zwart gat. De versie met berekende limieten (ipv. Newton-Rapson) geeft inderdaad de oplossing tot de benaderings radius.

Zonder jouw inzicht op de as limieten (de roots) hadden we via Newton-Rapson al veel minder Mikado oplossingen. Op dit moment ben ik nog geen Mikado oplossing gezien met berekende limieten dus niet door itteratie.

Zware massa lichtstraal op \(r_{zon}\):

Zware massa kleine radius (zon als referentie):

Limiet massa kleine radius (zon als referentie):

Via mybinder:
https://mybinder.org/v2/gh/oooVincentoo ... ptic.ipynb

Python Console (hoop dat deze werkt bij jouw):

Code: Selecteer alles

#Version 17.0 Schwarzschild, Jacobi Elliptic
#https://sciencetalk.nl/forum/viewtopic.php?f=85&t=212427

#Open pyplot in separate interactive window
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'qt5')

import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sps
from scipy import optimize
from matplotlib.widgets import Slider, CheckButtons                    

fig= plt.figure(figsize=(15, 6))  

ax1 = plt.subplot2grid((3, 3), (0, 0),rowspan=2)
ax2 = plt.subplot2grid((3, 3), (0, 1),rowspan=2)
ax3 = plt.subplot2grid((3, 3), (0, 2),rowspan=2)
ax4 = plt.subplot2grid((3, 3), (2, 0),colspan=3)

ax4.axis('off')

plt.suptitle('Schwarzschild, Jacobi Elliptic', fontsize=22)

# Draw slider alpha
axcolor = 'lightgoldenrodyellow'
axalpha = plt.axes([0.25, 0.05,0.4, 0.04], facecolor=axcolor)
s_zoom = Slider(axalpha, 'zoom: ', 0, 1, valinit=1)

axsm = plt.axes([0.25, 0.1,0.4, 0.04], facecolor=axcolor)
s_sm = Slider(axsm, 'solar masses: ', 1,100000, valinit=1)

radius = plt.axes([0.25, 0.15,0.4, 0.04], facecolor=axcolor)
s_sr = Slider(radius, 'radius: ', 0.1, 5, valinit=1)

rax = plt.axes([0.75, 0.05, 0.15, 0.15])
check = CheckButtons(rax, ('xy equal', 'sun'), (False, False))

#Set constants note G is divided by c^2    
global M
global G
global Ro

M=1.989e30
G=6.67408e-11/((3e8)**2)
Ro=696340000

def e1(sm,sr):
    esqrt=np.sqrt((sr*Ro-sm*2*M*G)*(sr*Ro+sm*6*M*G))
    e1=(sr*Ro-sm*2*M*G+esqrt)/(sm*4*M*G*sr*Ro)
    return e1

def e2(sr):
    e2=1/(sr*Ro)
    return e2

def e3(sm,sr):
    esqrt=np.sqrt((sr*Ro-sm*2*M*G)*(sr*Ro+sm*6*M*G))
    e3=(sr*Ro-sm*2*M*G-esqrt)/(sm*4*M*G*sr*Ro)
    return e3   
    
def r(phi):
    sm = s_sm.val
    sr = s_sr.val
     
    e1v=e1(sm,sr)
    e2v=e2(sr)
    e3v=e3(sm,sr)
    
    tau=np.sqrt(sm*M*G*(e1v-e3v)/2)
    h = (e2v-e3v)/(e1v-e3v)
    sigma = -tau*np.pi/2 - sps.ellipk(h)  
    sn,_,_,_=sps.ellipj(tau*phi+sigma, h)
    r=e3v+(e2v-e3v)*sn**2
    return r

def alpha():
    #Determine the root of solution
    sm = s_sm.val
    sr = s_sr.val
     
    e1v=e1(sm,sr)
    e2v=e2(sr)
    e3v=e3(sm,sr)
    
    h = (e2v-e3v)/(e1v-e3v)   
    phi=np.arcsin(np.sqrt(-e3v/(e2v-e3v)))    
    angle= sps.ellipkinc(phi,  h)  
    
    tau=np.sqrt(sm*M*G*(e1v-e3v)/2)
    sigma = -tau*np.pi/2 - sps.ellipk(h)  
    
    alphav = (angle-sigma) /tau-2*np.pi
 
    return alphav

#The following function updateplot runs whenever slider changes value
def updateplot(val):
    
    clearsetplot()

    #Solar masses/radius and magnification alpha
    zoom = s_zoom.val
    sm = s_sm.val
    sr = s_sr.val
   
    #Find root angle for non-solutions   
    alphamax=alpha()
    print('Direct:' + str(alphamax))
    
    #Set angle range and calculate radius
    phi= np.linspace(alphamax/2,np.pi-alphamax/2,500)
    phi=phi[int(0.5*(1-zoom)*500):int(0.5*(1+zoom)*500)]
   
    rv=1/r(phi)
    
    x=rv*np.cos(phi)
    y=rv*np.sin(phi)
    
    #plot lightpath curve
    ax1.plot(x,y)
    ax1.grid()
            
    #Angular Deflection
    dx=np.diff(x)
    dy=np.diff(y)
    dydx=dy/dx
    da=np.arctan(dydx)
    dam=np.max(da)
    ax2.plot(x[:-1],-da)
    ax2.grid()
    
    #Angular Distribution
    ddy=np.diff(da)
    ddydx=ddy/dx[:-1]
    ax3.plot(x[:-2],-ddydx)        
    ax3.grid()       
     
    text=(  str('{:10.6f}'.format(dam))      +   ' $[rad]$\n' +
            str('{:10.6f}'.format(np.degrees(dam)))      +   '$^{\circ}$\n' +
            str('{:10.6f}'.format(np.degrees(dam)*3600)) +  "$''$")         
    ax1.text(0.01, 0.75,text, transform=ax1.transAxes, fontsize=10)
    
    #scaling method and draw sun
    plotoptions(x,y,sr)
    fig.canvas.draw_idle()
  
#This function set scaling axis and drawes sun
def plotoptions(x,y,sr):
    i=0
    for r in check.get_status():
       
        if i==0:
            #Set scaling auto or equal
            if r==False:
                ax1.axis('auto')
                range=np.max(y)-np.min(y)
                ax1.set_ylim([np.min(y)-range*2/100,np.max(y)+range*2/100])
                range=np.max(x)-np.min(x)
                ax1.set_xlim([np.min(x)-range*2/100,np.max(x)]+range*2/100)
                ax2.set_xlim([np.min(x)-range*2/100,np.max(x)]+range*2/100)
                ax3.set_xlim([np.min(x)-range*2/100,np.max(x)]+range*2/100)
                
                
            else:
                ax1.axis('equal')
                range=np.max(x)-np.min(x)
                ax1.set_xlim([np.min(x)-range*2/100,np.max(x)+range*2/100])
                ax2.set_xlim([np.min(x)-range*2/100,np.max(x)+range*2/100])
                ax3.set_xlim([np.min(x)-range*2/100,np.max(x)+range*2/100])
                                
        if i==1:
            #Draw Sun yes/no
              if r==True:
                circle=plt.Circle((0,np.max(y)-sr*Ro),Ro, color='orange',label='Label')
                ax1.add_patch(circle)
                
                circle=plt.Circle((0,np.max(y)-sr*Ro),sr*Ro,fill=False, ls='--',ec='black')
                ax1.add_patch(circle)
                
        i=i+1

def clearsetplot():
    ax1.clear()
    ax1.set_title('Light path', fontsize=14)
    ax1.set_xlabel('x [m]', fontsize=14)
    ax1.set_ylabel('y [m]', fontsize=14)
    ax1.ticklabel_format(axis='both', style='sci', scilimits=(0,0))
    
    ax2.clear()
    ax2.set_title('Angular Deflection', fontsize=14)
    ax2.set_xlabel('x [m]', fontsize=14)
    ax2.set_ylabel('deflection angle  [rad]', fontsize=14)
    ax2.ticklabel_format(axis='both', style='sci', scilimits=(0,0))
    
    ax3.clear()
    ax3.set_title('Angular Disribution', fontsize=14)
    ax3.set_xlabel('x [m]', fontsize=14)
    ax3.set_ylabel('angular Change [rad/m]', fontsize=14)
    ax3.ticklabel_format(axis='both', style='sci', scilimits=(0,0))
       
#Display plot n startup.
updateplot(1)

#The following code checkes if slider changes. This line is looped automatic by pyplot
s_sm.on_changed(updateplot);
s_sr.on_changed(updateplot);
s_zoom.on_changed(updateplot);
check.on_clicked(updateplot);

[Professor Puntje]

Werkt nu allebei.

[OOOVincentOOO]

Even nog een update. Het is niet de Schwarzschild radius het is de Photon sphere [Wiki]!

Schwardschild Radius (zwarte cirkel):
$$r_s=\frac{2GM}{c^2}$$
Photon Sphere (grijze cirkel):
$$r_p=\frac{3GM}{c^2}$$
Heb nog wat kleine aanpassingen gedaan:

Zon 100000 maal zo zwaar:

Op iets kleinere afstand (zon als schaal ingetekend):

Limiet bijna aan photonsphere (zon als schaal ingetekend):

De limiet is dus net iets groter dan Photonspere kan dat? Het kan ook komen door naukeurigheid intekenen circel (het zijn grote getallen).

https://mybinder.org/v2/gh/oooVincentoo ... ptic.ipynb

Code: Selecteer alles

#Version 18.0 Schwarzschild, Jacobi Elliptic
#https://sciencetalk.nl/forum/viewtopic.php?f=85&t=212427

#Open pyplot in separate interactive window
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'qt5')

import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sps
from scipy import optimize
from matplotlib.widgets import Slider, CheckButtons                    

fig= plt.figure(figsize=(15, 6))  

ax1 = plt.subplot2grid((3, 3), (0, 0),rowspan=2)
ax2 = plt.subplot2grid((3, 3), (0, 1),rowspan=2)
ax3 = plt.subplot2grid((3, 3), (0, 2),rowspan=2)
ax4 = plt.subplot2grid((3, 3), (2, 0),colspan=3)

ax4.axis('off')

plt.suptitle('Schwarzschild, Jacobi Elliptic', fontsize=22)

# Draw slider alpha
axcolor = 'lightgoldenrodyellow'
axalpha = plt.axes([0.25, 0.05,0.4, 0.04], facecolor=axcolor)
s_zoom = Slider(axalpha, 'zoom: ', 0, 1, valinit=1)

axsm = plt.axes([0.25, 0.1,0.4, 0.04], facecolor=axcolor)
s_sm = Slider(axsm, 'solar masses: ', 1,100000, valinit=1)

radius = plt.axes([0.25, 0.15,0.4, 0.04], facecolor=axcolor)
s_sr = Slider(radius, 'radius: ', 0.1, 5, valinit=1)

rax = plt.axes([0.75, 0.05, 0.15, 0.15])
check = CheckButtons(rax, ('xy equal', 'sun'), (False, False))

#Set constants note G is divided by c^2    
global M
global G
global Ro

M=1.989e30
G=6.67408e-11/((3e8)**2)
Ro=696340000

def e1(sm,sr):
    esqrt=np.sqrt((sr*Ro-sm*2*M*G)*(sr*Ro+sm*6*M*G))
    e1=(sr*Ro-sm*2*M*G+esqrt)/(sm*4*M*G*sr*Ro)
    return e1

def e2(sr):
    e2=1/(sr*Ro)
    return e2

def e3(sm,sr):
    esqrt=np.sqrt((sr*Ro-sm*2*M*G)*(sr*Ro+sm*6*M*G))
    e3=(sr*Ro-sm*2*M*G-esqrt)/(sm*4*M*G*sr*Ro)
    return e3   
    
def r(phi):
    sm = s_sm.val
    sr = s_sr.val
     
    e1v=e1(sm,sr)
    e2v=e2(sr)
    e3v=e3(sm,sr)
    
    tau=np.sqrt(sm*M*G*(e1v-e3v)/2)
    h = (e2v-e3v)/(e1v-e3v)
    sigma = -tau*np.pi/2 - sps.ellipk(h)  
    sn,_,_,_=sps.ellipj(tau*phi+sigma, h)
    r=e3v+(e2v-e3v)*sn**2
    return r

def alpha():
    #Determine the root of solution
    sm = s_sm.val
    sr = s_sr.val
     
    e1v=e1(sm,sr)
    e2v=e2(sr)
    e3v=e3(sm,sr)
    
    h = (e2v-e3v)/(e1v-e3v)   
    phi=np.arcsin(np.sqrt(-e3v/(e2v-e3v)))    
    angle= sps.ellipkinc(phi,  h)  
    
    tau=np.sqrt(sm*M*G*(e1v-e3v)/2)
    sigma = -tau*np.pi/2 - sps.ellipk(h)  
    
    alphav = (angle-sigma) /tau-2*np.pi
 
    return alphav

#The following function updateplot runs whenever slider changes value
def updateplot(val):
    
    clearsetplot()

    #Solar masses/radius and magnification alpha
    zoom = s_zoom.val
    sm = s_sm.val
    sr = s_sr.val
   
    #Find root angle for non-solutions   
    alphamax=alpha()
    
    #Set angle range and calculate radius
    phi= np.linspace(alphamax/2,np.pi-alphamax/2,500)
    phi=phi[int(0.5*(1-zoom)*500):int(0.5*(1+zoom)*500)]
   
    rv=1/r(phi)
    
    x=rv*np.cos(phi)
    y=rv*np.sin(phi)
    
    #plot lightpath curve
    ax1.plot(x,y)
    ax1.grid()
            
    #Angular Deflection
    dx=np.diff(x)
    dy=np.diff(y)
    dydx=dy/dx
    da=np.arctan(dydx)
    dam=np.max(da)
    ax2.plot(x[:-1],-da)
    ax2.grid()
    
    #Angular Distribution
    ddy=np.diff(da)
    ddydx=ddy/dx[:-1]
    ax3.plot(x[:-2],-ddydx)        
    ax3.grid()       
     
    text=(  str('{:10.6f}'.format(dam))      +   ' $[rad]$\n' +
            str('{:10.6f}'.format(np.degrees(dam)))      +   '$^{\circ}$\n' +
            str('{:10.6f}'.format(np.degrees(dam)*3600)) +  "$''$")         
    ax1.text(0.01, 0.75,text, transform=ax1.transAxes, fontsize=10)
    
    #scaling method and draw sun
    plotoptions(x,y,sr,sm)
    fig.canvas.draw_idle()
  
#This function set scaling axis and drawes sun
def plotoptions(x,y,sr,sm):
    i=0
    for r in check.get_status():
       
        if i==0:
            #Set scaling auto or equal
            if r==False:
                ax1.axis('auto')
                range=np.max(y)-np.min(y)
                ax1.set_ylim([np.min(y)-range*2/100,np.max(y)+range*2/100])
                range=np.max(x)-np.min(x)
                ax1.set_xlim([np.min(x)-range*2/100,np.max(x)]+range*2/100)
                ax2.set_xlim([np.min(x)-range*2/100,np.max(x)]+range*2/100)
                ax3.set_xlim([np.min(x)-range*2/100,np.max(x)]+range*2/100)
                
                
            else:
                ax1.axis('equal')
                range=np.max(x)-np.min(x)
                ax1.set_xlim([np.min(x)-range*2/100,np.max(x)+range*2/100])
                ax2.set_xlim([np.min(x)-range*2/100,np.max(x)+range*2/100])
                ax3.set_xlim([np.min(x)-range*2/100,np.max(x)+range*2/100])
                                
        if i==1:
            #Draw Sun yes/no
              if r==True:
                circle=plt.Circle((0,0),Ro, color='orange',label='Label')
                ax1.add_patch(circle)
                
                circle=plt.Circle((0,0),sr*Ro,fill=False, ls='--',ec='black')
                ax1.add_patch(circle)

                circle=plt.Circle((0,0),sm*3*G*M, color='gray',label='Label')
                ax1.add_patch(circle)

                
                circle=plt.Circle((0,0),sm*2*G*M, color='black',label='Label')
                ax1.add_patch(circle)                
                
        i=i+1

def clearsetplot():
    ax1.clear()
    ax1.set_title('Light path', fontsize=14)
    ax1.set_xlabel('x [m]', fontsize=14)
    ax1.set_ylabel('y [m]', fontsize=14)
    ax1.ticklabel_format(axis='both', style='sci', scilimits=(0,0))
    
    ax2.clear()
    ax2.set_title('Angular Deflection', fontsize=14)
    ax2.set_xlabel('x [m]', fontsize=14)
    ax2.set_ylabel('deflection angle  [rad]', fontsize=14)
    ax2.ticklabel_format(axis='both', style='sci', scilimits=(0,0))
    
    ax3.clear()
    ax3.set_title('Angular Disribution', fontsize=14)
    ax3.set_xlabel('x [m]', fontsize=14)
    ax3.set_ylabel('angular Change [rad/m]', fontsize=14)
    ax3.ticklabel_format(axis='both', style='sci', scilimits=(0,0))
       
#Display plot n startup.
updateplot(1)

#The following code checkes if slider changes. This line is looped automatic by pyplot
s_sm.on_changed(updateplot);
s_sr.on_changed(updateplot);
s_zoom.on_changed(updateplot);
check.on_clicked(updateplot);

[OOOVincentOOO]

De limiet voor de photon sphere uit de Jacobi Elliptic heb ik niet bepaald. Echter ik heb de betreffende radius terug gerekend en precies ingevuld in het model.

$$M=1,989e+30 \ [kg]$$
$$G/c^2=7,416e-28 \ [m/kg]$$
$$R_{sun}=6,963e+08 \ [m]$$

Massa factor: 100.000:
$$r_p=\frac{3GM}{c^2}=3(10e+4)(1,989e+30)(7,416e-28)=4,425e+08 \ [m]$$
$$\frac{r_p}{R_{sun}}=\frac{4,425e+08}{6,963e+08}=0.6356$$

Deze radius fractie kan ingevuld worden in het model. Dit kan door middel van de slider of in de code op de plaats (let de fractie is een klein beetje groter gekozen anders kom er "nul" uit de berekening):

Code: Selecteer alles

radius = plt.axes([0.25, 0.15,0.4, 0.04], facecolor=axcolor)
s_sr = Slider(radius, 'radius: ', 0.1, 5, valinit=0.636)

The photo sphere wordt precies gedefinieerd door de Jacobi Elliptic aanpak. Het wiskundige bewijs heb ik niet. Met de \(sn\) functie en de parameters \(e1\), \(e2\), \(e3\) is het moeilijk te overzien. In eerste instantie dacht ik ook dat de Schwarzschild radius de limiet zou zijn.

Ik ben meer een praktisch gericht met mijn wiskunde en denk makkelijk functies te herkennen. Merk op dat in de onderste plot een soort van \(arcsin\) verdeling te zien is (onderste figuur plaatje rechts: angular distribution merk de U=vorm op).

Speculatie: Dit type verdeling bij de photon sphere: de angular distribution \((arcsin)\) komt vaak ter sprake bij random walks / brownian motion. Gezien deze verdeling doet mij vermoeden dat deze omloop (photon sphere) instabiel (cumulatief divergeerd voor quanta welke stapjes \(x\) nemen) en statistisch van nature is. Maar dit is puur speculeren van iemand met geen ervaring in fysica van zwarte gaten noch GR.

[Professor Puntje]

Het is alweer een tijdje geleden dat ik een boek over zwarte gaten heb doorgewerkt, dus ik heb dat niet allemaal meer paraat. Op het moment ben ik ook nog even druk met de afronding van het piekenverhaal (waarom vindt MathPages die twee pieken en wij er maar één).

[OOOVincentOOO]

Dit zal een van mijn laatste berichten zijn. Hieronder een filmpje met object diameter zon en toenemende massa. Op een gegeven moment verandererd de deflectie hoek verdeling van "Gauss" (omgekeerd klok vormig) achtig naar "arcsine" distributie achtig (U vorm). Gedurende deze transitie ontstaan er kortstondig twee pieken.

Van mijn ervaring tot nu toe zou ik zeggen dat de twee pieken ontstaan afhankelijk welke benadering men neemt. In het integraal pad zou ik zeggen dat er oplossingen zijn met twee pieken met sommige benaderingen. Dus afhankelijk welke benadering men neemt betreffende de radius en massa.

Nu laat ik het misschien beter over aan anderen. Voorlopig ben ik voldaan met hetgeen ik geleerd heb. Excuses voor de vele plaatjes, maar dit is ten behoeve van discussie.

Het is beter de discussie "Twee Pieken" verder te voeren vanaf betreffende topic:
viewtopic.php?f=66&t=212342

[OOOVincentOOO]

Het probleem is dat er altijd details zijn waarmee men niet tevreden is. Ik heb een minder schokkerige/onrustige video gemaakt met hogere frame rate.

Men/ik blijf altijd details zien waarmee ik niet vertreden ben. En er komen altijd nieuwe details boven water.

Observaties:
Er zijn meerdere transities te herkennen:

• De "Gauss" verdeling veranderd. De buigpunten waar normaal "Standaard deviatie" zit worden positief in buiging [0:22]

• Verdeling begint te veranderen van concaaf naar convex [0:26]. Waarbij de verdeling van bol naar hol gaat.

• Verdeling splits op en veranderd langzaam in "x^2" parabool distributie [0:28]. Merk op lineaire verdeling Angular Deflection.

• Transitie storing: enkele frames is ruis te zien op de verdeling. Vanaf dit punt veranderd verdeling in "arcsin" distributie "U" vorm. [0:33]. Oorzaak ruis niet begrepen.

• Halve cirkel 180° op photon sphere (3 seconden frame bevroren). Enkelzijdige "arcsin". [0:35].

• Hele cirkel 360° afbuiging (3 seconden frame bevroren). Dubbelzijdige "arcsin". Dit is laatste frame, buiten deze zijn er geen oplossingen op basis van de methode.

[OOOVincentOOO]

Het laatste filmpje was beter maar nog niet tot mijn genoegen. Men moet op veel details letten een stabiel filmpje te genereren. Hierbij een verbeterde versie Jacobi met Python:

Het principe is dat de massa van de zon steeds groter wordt. Met behulp van het Schwarzschild geometry met Jacobi elliptics word het licht pad bepaald. Wanneer de massa vergroot bereikt de lichtstraal uiteindelijk de photo sphere.

Dit filmpje heb ik ingezoomd op 90%. Zo ziet men de zon met standaard massa bij aanvang, hier is Schwarzschild radius verwaarloosbaar. De massa van de zon zal natuurlijk nooit echt toenemen. Dit is slechts een berekening wat een lichtstraal doet in nabijheid van een (grote) massa.



De code te vinden in comments youtube of:
https://mybinder.org/v2/gh/oooVincentoo ... ptic.ipynb

Graag wil nog weten wat er gebeurt bij transitie [0:xx] (edit: andere tijd index) hier ontstaat "ruis". Dit heeft waarschijnlijk te doen met de Jacobi elliptic functies.

Deze bevindingen zijn nog steeds geldig:

Observaties:
Er zijn meerdere transities te herkennen:

• De "Gauss" verdeling veranderd. De buigpunten waar normaal "Standaard deviatie" zit worden positief in buiging [0:22]

• Verdeling begint te veranderen van concaaf naar convex [0:26]. Waarbij de verdeling van bol naar hol gaat.

• Verdeling splits op en veranderd langzaam in "x^2" parabool distributie [0:28]. Merk op lineaire verdeling Angular Deflection.

• Transitie storing: enkele frames is ruis te zien op de verdeling. Vanaf dit punt veranderd verdeling in "arcsin" distributie "U" vorm. [0:33]. Oorzaak ruis niet begrepen.

• Halve cirkel 180° op photon sphere (3 seconden frame bevroren). Enkelzijdige "arcsin". [0:35].

• Hele cirkel 360° afbuiging (3 seconden frame bevroren). Dubbelzijdige "arcsin". Dit is laatste frame, buiten deze zijn er geen oplossingen op basis van de methode.

[Professor Puntje]

In verband met je experimenten: kijk eens wat er met e₁ , e₃ en h gebeurt voor r₀ = r_s.

[OOOVincentOOO]

Die afstand valt niet in de oplossingen. Alles bevind zich buiten photon sphere zover ik weet. Eigenlijk logisch het is een model voor licht.

Echter er gebeurt een transitie van verdeling Gauss naar arcsin U achtig vlak voor photon spere. Via twee bulten / pieken. En een klein gebied waar ruis optreed.

De getallen e1, e3 en h weet ik niet. Die zou ik dan dienen te plotten / weer tegeven.

[Professor Puntje]

Weet je voor welke waarden van r₀/r_s de ruis en transitie optreden? Dan kan ik kijken of er daar iets vreemds met de formules gebeurt.

[OOOVincentOOO]

Daar ben ik nog niet aan toegekomen.

[OOOVincentOOO]

\(r/r_s=2\): Dit is het moment waarbij de verdeling van "Gauss" klokvorm naar hole U vorm gaat. Kort hierna treden twee pieken op [0:37] s.

\(r/r_p=1\): Er treedt kortstondig ruis op wanneer we bij de photon sphere zijn. Het bereik x is erg klein en het traject lijkt even horizontaal. Maar kort hierna zijn we bij de photon sphere. Hier is de eccentricity \(h^2=1\) een cirkel dus.
$$\boxed{ r = \frac{1}{ e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \varphi + \sigma ; h^2)}} \tag{6}$$
We bekijken alleen de oplossing tussen de nulpunten van de noemer van formule \((6)\). Ik denk dat we die criteria iets te strikt hanteren. Ik heb zojuist getest voor bereik: \([0,\pi]\) en krijg dezelfde oplossingen de lichtstraal loopt dan door tot x-as.

[Professor Puntje]

In de buurt van h^2 = 1 kun je problemen met de nauwkeurigheid krijgen. Zie:

https://docs.scipy.org/doc/scipy/refere ... llipj.html

En de elliptische integralen uit hetzelfde pakket.