Python en Jacobi

[Professor Puntje]

Waarschijnlijk gaat het wel hiermee:

Want we hebben:

Code: Selecteer alles

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)
h = np.sqrt((e2 - e3)/(e1 - e3))
m = h*h
print((-e3)/(e2-e3))
print(np.sqrt((-e3)/(e2-e3)))
print(m*(-e3)/(e2-e3))

0.49999894643
0.7071060362
4.21425019424e-06

De precieze formule voor de totale afbuiging is minder eenvoudig dan ik gehoopt had. Dus die zal ik met de nodige zorgvuldigheid stapje voor stapje afleiden.

[OOOVincentOOO]

Ik kijk jouw resultaten morgen verder. Die integraal is toch de Elliptic integraal? Volgens mij staat die in hoofdstuk 9 waar hij refereerd naar Davis [4] dacht ik.

Ik heb verder geprutst een interactieve versie te maken. Voor sommige extremen krijg oplossingen buiten boundary's. Zullen we dat het Mikado effect noemen? Insider joke jij zult ook veel Mikado oplossingen gezien hebben!

Python is niet echt de programmeertaal waarmee ik interfaces kan maken. Maar met een beetje improviseren. Ik hoop dat he werkt. Misschien kan jij nog limieten vinden voor het zoomen (de alpha).

Code: Selecteer alles

#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 matplotlib.widgets import Slider, CheckButtons
fig, ax= plt.subplots(2,1,figsize=(9, 5),gridspec_kw={'height_ratios': [5, 1]})                      
                                          
                      
# Draw slider alpha
axcolor = 'lightgoldenrodyellow'
axalpha = plt.axes([0.25, 0.05,0.4, 0.04], facecolor=axcolor)
s_alpha = Slider(axalpha, 'alpha: ', 0, 2, 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, 'sun radius: ', 1, 10, valinit=1)

rax = plt.axes([0.75, 0.05, 0.1, 0.15])
check = CheckButtons(rax, ('x=y scale', 'sun'), (False, False))

#Set constants note G is divided by c^2    
M=1.989e30
G=6.67408e-11/((3e8)**2)
Ro=7e8   

ax[1].axis('off')
 
#The following function updateplot runs whenever slider changes value
def updateplot(val):
    
    ax[0].clear()
    
    #Solar masses/radius and magnification alpha
    alpha = s_alpha.val
    sm = s_sm.val
    sr = s_sr.val
    
    phi=np.linspace(np.pi/2-alpha*np.pi/2,np.pi/2+alpha*np.pi/2,1000)

    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)
    #print('e1: ' + str(e1))
    
    e2=1/(sr*Ro)
    #print('e2: ' + str(e2))
    
    e3=(sr*Ro-sm*2*M*G-esqrt)/(sm*4*M*G*sr*Ro)
    #print('e3: ' + str(e3))
    
    Tau=np.sqrt(sm*M*G*(e1-e3)/2)
    #print('Tau: ' + str(Tau))
    
    h=(e2-e3)/(e1-e3)
    #print('h: ' + str(h))
    
    #Solve angle start conditions
    angle=sps.ellipk(h)
    #print('angle: ' + str(angle))
    
    #Calculate Sigma Start Condition
    Sigma=-Tau*np.pi/2-angle
    #print('sigma: ' + str(Sigma))
     
    v=Tau*phi+Sigma
    sn,_,_,_=sps.ellipj(v,h)
    r=1/(e3+(e2-e3)*(sn**2))
    
    x=r*np.cos(phi)
    y=r*np.sin(phi)
    deflection, = ax[0].plot(x,y)

    i=0
    for r in check.get_status():
       
        if i==0:
            #Set scaling auto or equal
            if r==True:
                plt.autoscale() 
            else:
                ax[0].axis('equal')
        else:
            #Draw Sun yes/no
              if r==True:
                circle=plt.Circle((0,np.max(y)-sr*Ro),sr*Ro, color='orange')
                ax[0].add_patch(circle)
        
        i=i+1
    
    ax[0].grid()   
    fig.canvas.draw_idle()
    
updateplot(1)

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

[OOOVincentOOO]

Update code. De scaling van x=y equal loopt nog niet probleemloos. Hier een verbeterde versie:

Code: Selecteer alles

#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 matplotlib.widgets import Slider, CheckButtons
fig, ax= plt.subplots(2,1,figsize=(9, 5),gridspec_kw={'height_ratios': [5, 1]})                        
                                          
                      
# Draw slider alpha
axcolor = 'lightgoldenrodyellow'
axalpha = plt.axes([0.25, 0.05,0.4, 0.04], facecolor=axcolor)
s_alpha = Slider(axalpha, 'alpha: ', 0, 2, 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, 'sun radius: ', 1, 10, valinit=1)

rax = plt.axes([0.75, 0.05, 0.1, 0.15])
check = CheckButtons(rax, ('x=y scale', '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=7e8   

ax[1].axis('off')
 
#The following function updateplot runs whenever slider changes value
def updateplot(val):
    
    ax[0].clear()
    
    #Solar masses/radius and magnification alpha
    alpha = s_alpha.val
    sm = s_sm.val
    sr = s_sr.val
    
    phi=np.linspace(np.pi/2-alpha*np.pi/2,np.pi/2+alpha*np.pi/2,1000)

    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)
    #print('e1: ' + str(e1))
    
    e2=1/(sr*Ro)
    #print('e2: ' + str(e2))
    
    e3=(sr*Ro-sm*2*M*G-esqrt)/(sm*4*M*G*sr*Ro)
    #print('e3: ' + str(e3))
    
    Tau=np.sqrt(sm*M*G*(e1-e3)/2)
    #print('Tau: ' + str(Tau))
    
    h=(e2-e3)/(e1-e3)
    #print('h: ' + str(h))
    
    #Solve angle start conditions
    angle=sps.ellipk(h)
    #print('angle: ' + str(angle))
    
    #Calculate Sigma Start Condition
    Sigma=-Tau*np.pi/2-angle
    #print('sigma: ' + str(Sigma))
     
    v=Tau*phi+Sigma
    sn,_,_,_=sps.ellipj(v,h)
    r=1/(e3+(e2-e3)*(sn**2))
    
    x=r*np.cos(phi)
    y=r*np.sin(phi)
    deflection, = ax[0].plot(x,y)

    i=0
    for r in check.get_status():
       
        if i==0:
            #Set scaling auto or equal
            if r==False:
                plt.autoscale()                 
            else:
                ax[0].axis('equal')
        else:
            #Draw Sun yes/no
              if r==True:
                circle=plt.Circle((0,np.max(y)-sr*Ro),sr*Ro, color='orange')
                ax[0].add_patch(circle)
        
        i=i+1
    
    ax[0].grid()   
    fig.canvas.draw_idle()
    
updateplot(1)

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

[Professor Puntje]

Ik heb verder geprutst een interactieve versie te maken. Voor sommige extremen krijg oplossingen buiten boundary's. Zullen we dat het Mikado effect noemen? Insider joke jij zult ook veel Mikado oplossingen gezien hebben!

Ja dat is goed - heel toepasselijk!

Die inzoom functie is ook leuk.

[Professor Puntje]

Hoe los je φ uit onderstaande vergelijking voor gegeven waarden e₂, e₃, τ, σ en h in Python numeriek op?

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

[CoenCo]

Met Newton-Rhapson o.i.d.
Kijk eens of je hier wat mee kan: https://docs.scipy.org/doc/scipy/refere ... ewton.html

[OOOVincentOOO]

Ik kijk vanavond een naar jouw bericht. Zo op eerste gezicht lijkt dat toch op jouw methode om de startpositie uit te rekenen?

Update van versie. De knoppen xy equal en sun on/off lijken nu te werken. Het is wel een zoekwerk tot de simpelste oplossing te komen. Alle mislukkingen ingewikkelde (soms mooie) code is nu weg. Nu een beetje aan details werken: labels toevoegen, misschien deflectiehoek weergeven.

Code: Selecteer alles

#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 matplotlib.widgets import Slider, CheckButtons                    

fig = plt.figure(constrained_layout=False)
gs = fig.add_gridspec(6, 1)
ax = fig.add_subplot(gs[0 :4])

                                 
                      
# Draw slider alpha
axcolor = 'lightgoldenrodyellow'
axalpha = plt.axes([0.25, 0.05,0.4, 0.04], facecolor=axcolor)
s_alpha = Slider(axalpha, 'zoom: ', 0, 1.5, 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, 'sun radius: ', 1, 10, 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=7e8   

#The following function updateplot runs whenever slider changes value
def updateplot(val):
    
    ax.clear()
    ax.set_title('Light Deflection (Jacobi Elliptic)')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.autoscale(enable=False)
    
    #Solar masses/radius and magnification alpha
    alpha = s_alpha.val
    sm = s_sm.val
    sr = s_sr.val
    
    phi=np.linspace(np.pi/2-alpha*np.pi/2,np.pi/2+alpha*np.pi/2,1000)

    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)
    #print('e1: ' + str(e1))
    
    e2=1/(sr*Ro)
    #print('e2: ' + str(e2))
    
    e3=(sr*Ro-sm*2*M*G-esqrt)/(sm*4*M*G*sr*Ro)
    #print('e3: ' + str(e3))
    
    Tau=np.sqrt(sm*M*G*(e1-e3)/2)
    #print('Tau: ' + str(Tau))
    
    h=(e2-e3)/(e1-e3)
    #print('h: ' + str(h))
    
    #Solve angle start conditions
    angle=sps.ellipk(h)
    #print('angle: ' + str(angle))
    
    #Calculate Sigma Start Condition
    Sigma=-Tau*np.pi/2-angle
    #print('sigma: ' + str(Sigma))
     
    v=Tau*phi+Sigma
    sn,_,_,_=sps.ellipj(v,h)
    r=1/(e3+(e2-e3)*(sn**2))
    
    x=r*np.cos(phi)
    y=r*np.sin(phi)

    #plot curve
    ax.plot(x,y)
    ax.grid()   

    #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:
                ax.axis('auto')
                range=np.max(y)-np.min(y)
                ax.set_ylim([np.min(y)-range*2/100,np.max(y)+range*2/100])
                range=np.max(x)-np.min(x)
                ax.set_xlim([np.min(x)-range*2/100,np.max(x)]+range*2/100)
            else:
                ax.axis('equal')
                range=np.max(x)-np.min(x)
                ax.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),sr*Ro, color='orange')
                ax.add_patch(circle)
        i=i+1

#Display plot n startup.
updateplot(1)

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

[Professor Puntje]

@ OOOVincentOOO

Je laatste versie van je programma geeft bij mij een foutmeiding:

Code: Selecteer alles

TypeError: __init__() got an unexpected keyword argument 'constrained_layout'

[Professor Puntje]

Een vermoedelijk exacte formule voor de totale lichtafbuiging is formule (6.16) in Solomon Antoniou's artikel.

Het oerwoud aan verschillende notaties voor elliptische integralen en functies is mij een doorn in het oog, en heeft in dit topic al tot de nodige fouten geleid. Hopelijk zijn we nu op de goede weg...

[Professor Puntje]

Code: Selecteer alles

from scipy import optimize
import numpy as np
from scipy.special import ellipj
from scipy.special import ellipk


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)
h = np.sqrt((e2 - e3)/(e1 - e3))
m = h*h
sigma = -tau*(np.pi/2) + ellipk(m)


t = np.linspace(-0,5*np.pi,0,100000)
sn,_,_,_=ellipj(tau*t + sigma , m)

def f(t):

    return e3 + (e2 - e3)*sn**2

root = optimize.newton(f,-8e-6)

print(root)

Geeft:

Code: Selecteer alles

RuntimeError: Failed to converge after 50 iterations, value is []

Wat doe ik fout?

[CoenCo]

De berekening van sn ( en alle bijbehorende parameters) moet in f(t) gebeuren. Nu is f(t) namelijk helemaal geen functie van t.

probeer dit eens:

Code: Selecteer alles

from scipy import optimize
import numpy as np
from scipy.special import ellipj
from scipy.special import ellipk





def f(t):
    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)
    h = np.sqrt((e2 - e3)/(e1 - e3))
    m = h*h
    sigma = -tau*(np.pi/2) + ellipk(m)
    sn,_,_,_=ellipj(tau*t + sigma , m)
    return e3 + (e2 - e3)*sn**2

root = optimize.newton(f,-8e-6)

print(root)

[Professor Puntje]

Mooi! Dank - dat werkt.

[Professor Puntje]

Goed - dit is het dan voor zover de lichtbaan aangaat:

Code: Selecteer alles

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




# Berekening totale lichtafbuigng alpha

def f(t):
    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)
    h = np.sqrt((e2 - e3)/(e1 - e3))
    m = h*h
    sigma = -tau*(np.pi/2) + ellipk(m)
    sn,_,_,_=ellipj(tau*t + sigma , m)
    return e3 + (e2 - e3)*sn**2

root = optimize.newton(f,-1e-6)
alpha = -2*root



# Berekening lichtbaan

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)
h = np.sqrt((e2 - e3)/(e1 - e3))
m = h*h
sigma = -tau*(np.pi/2) + ellipk(m)

alpha99 = 0.99*alpha # beveiliging tegen singulariteiten
t = np.linspace(-alpha99/2,np.pi+alpha99/2,1000)

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

x = (1/(e3 + (e2 - e3)*sn**2))*np.cos(t)
y = (1/(e3 + (e2 - e3)*sn**2))*np.sin(t)

plt.plot(x,y,c=(1,0.2,0.5),lw=1)
plt.title('lichtbaan')
plt.axis('on')
plt.show()


# De berekende constanten

print("root =", root)
print("alpha =", alpha)
print("e1 =", e1)
print("e2 =", e2)
print("e3 =", e3)
print("tau =", tau)
print("h =", h)
print("m =", m)
print("sigma =", sigma)
print("alpha =", alpha)

[OOOVincentOOO]

Hallo,

Mooi, straks jouw methode voor de alpha limiet bekijken. Kijken of ik het ga snappen. Ben meer aan het programmeren en tweaken tussendoor.

Mijn code heb nog aangepast. Hopelijk werkt het nu wel. Misschien heb ik onbewust routines uit Jupyter notebook overgenomen. Ik gebruik Spyder vanuit Anaconda installatie.

De zon heeft nu altijd de vaste diameter. Ik heb een kleine circel eromheen getekend met afstand lichtstraal. Dan is het hopelijk beter te begrijpen.

Deflectie hoek berekening (enkele hoek): Via numerieke diff. en de maximale waarde genomen.

Het lijkt erop dat de Mikado toestanden ook de goede hoek geven. Maar dien hier meer te testen.

Code: Selecteer alles

#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 matplotlib.widgets import Slider, CheckButtons                    

fig, ax= plt.subplots(2,1,figsize=(9, 5),gridspec_kw={'height_ratios': [5, 1]})  

#Do not display second dummy graph (position widgets)
ax[1].axis('off')                            
                      
# Draw slider alpha
axcolor = 'lightgoldenrodyellow'
axalpha = plt.axes([0.25, 0.05,0.4, 0.04], facecolor=axcolor)
s_alpha = Slider(axalpha, 'zoom: ', 0, 1.5, 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, 'sun radius: ', 1, 10, 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=7e8   

#The following function updateplot runs whenever slider changes value
def updateplot(val):
    
    ax[0].clear()
    ax[0].set_title('Light Deflection (Jacobi Elliptic)', fontsize=20)
    ax[0].set_xlabel('x', fontsize=14)
    ax[0].set_ylabel('y', fontsize=14)
    ax[0].autoscale(enable=False)

    #Solar masses/radius and magnification alpha
    alpha = s_alpha.val
    sm = s_sm.val
    sr = s_sr.val

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

    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)
    #print('e1: ' + str(e1))
    
    e2=1/(sr*Ro)
    #print('e2: ' + str(e2))
    
    e3=(sr*Ro-sm*2*M*G-esqrt)/(sm*4*M*G*sr*Ro)
    #print('e3: ' + str(e3))
    
    Tau=np.sqrt(sm*M*G*(e1-e3)/2)
    #print('Tau: ' + str(Tau))
    
    h=(e2-e3)/(e1-e3)
    #print('h: ' + str(h))
    
    #Solve angle start conditions
    angle=sps.ellipk(h)
    #print('angle: ' + str(angle))
    
    #Calculate Sigma Start Condition
    Sigma=-Tau*np.pi/2-angle
    #print('sigma: ' + str(Sigma))
     
    v=Tau*phi+Sigma
    sn,_,_,_=sps.ellipj(v,h)
    r=1/(e3+(e2-e3)*(sn**2))
    
    x=r*np.cos(phi)
    y=r*np.sin(phi)

    #plot curve
    ax[0].plot(x,y)
    ax[0].grid()   
    
    #Mean Deflection Angle
    dx=np.diff(x)
    dy=np.diff(y)
    da=np.arctan(dy/dx)
    da=np.abs(da)
    dam=np.degrees(np.max(da))
    
    text=('Deflection:\n' +
            str('{:08.6f}'.format(dam))      +   '$^{\circ}$\n' +
            str('{:08.2f}'.format(dam*3600)) +  "$''$")
    ax[0].text(0.01, 0.75,text, transform=ax[0].transAxes, fontsize=14)
    
    

    #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:
                ax[0].axis('auto')
                range=np.max(y)-np.min(y)
                ax[0].set_ylim([np.min(y)-range*2/100,np.max(y)+range*2/100])
                range=np.max(x)-np.min(x)
                ax[0].set_xlim([np.min(x)-range*2/100,np.max(x)]+range*2/100)
            else:
                ax[0].axis('equal')
                range=np.max(x)-np.min(x)
                ax[0].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),sr*Ro,fill=False, ls='--',ec='black')
                ax[0].add_patch(circle)
                
                circle=plt.Circle((0,np.max(y)-sr*Ro),Ro, color='orange')
                ax[0].add_patch(circle)
                
        i=i+1

#Display plot n startup.
updateplot(1)

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

[Professor Puntje]

Ja - je programma werkt nu bij mij ook.

Als je al begonnen bent met dy/dx van de grafiek te bepalen zou ik dat stukje code graag overnemen, want een plaatje van de afgeleide van de grafiek hebben we ook nodig voor het pieken onderzoek.