Python en Jacobi

[Professor Puntje]

Ook hebben we geen bewijs dat de deflectie convergeerd met de Jacobi methode. Of weet jij een handigheid?

Dat kun je zien uit formule (6) of (6*). De voerstraal r gaat naar oneindig als de noemer naar 0 gaat. Dus dat levert een asymptoot. De Jacobi elliptische functie hoeft daarvoor enkel voor het toepasselijke argument continu te zijn.

[OOOVincentOOO]

Oke thanks, goed gezien, even een overdosis gehad de laatste dagen.

[Professor Puntje]

Maar goed, we hebben veel werk verzet.

[OOOVincentOOO]

Mijn ervaring: fouten maken (en je dom voelen) de diepe dalen en veel leren en dingen ontdekken, de hoge pieken! Leuke aanvullingen om samen tot resultaat te komen. Knap dat je hebt doorgezet en de essenties van Python goed te pakken hebt!

Mijn mooiste ervaringen:

• Werken gestandaardiseerd werkt niet: \(R=1\), \(M=1\) etc.

• Onzin grafieken indien phi bereik verkeerd (Mikado).

• Methode om phi bereik te bepalen eerst via Wolfram daarna in Python Jacobi integral.

• Deflektie klopt niet vele ordes verkeerd: variabele \(h\) dient \(h^{2} \)te zijn.

• Boundary voor oplossingen roots via Newton-Rapson.

• Grafiek NIET integreren maar differentiëren. Het differentiëren is eigenlijk vorm van integreren omdat eerst \(arctan\) genomen word.

• Deflectie hoek komt overeen op 4 decimalen.

• Voor het eerst coderen in Spyder normaal doe ik dit alleen in Jupyter notebook (meestal niet interactief zonder sliders).

EDIT:
Dankjwel @CoenCo voor support op juiste moment, en het zwijgen op juiste moment om ons zelf te laten ploeteren!

[Professor Puntje]

Dit was voor mij onverwachts een stoomcursus Python, terwijl ik eigenlijk met Octave bezig was. Maar ik heb er geen spijt van, want Python is ook een heel mooie programmeertaal. Fijn om achter de hand te hebben als de formules te ingewikkeld worden, en met de hand oplossen niet meer gaat. Er waren in dit topic inderdaad de nodige hobbels te nemen, maar het is ons gelukt. Dank ook aan CoenCo voor de tips.

[Marko]

@ Marko

Ja - dat klopt. De benaderingen die hier nog inzitten komen uit de toepassing van numerieke methoden (oftewel het gebruik van Python). Ik had de zaak ook liever helemaal analytisch aangepakt, maar dat is me niet gelukt. De formules werden daar te ingewikkeld voor. Hoe nauwkeurig is Python? Wellicht kunnen anderen met meer ervaring met dat programma daar beter op antwoorden. Uit de berichten hier op het Wetenschapsforum begreep ik dat het nauwkeurig genoeg is voor een simulatie als in dit topic.

Een eerder bericht ging over de numerieke precisie van de floating point bewerkingen in Python. Die is meer dan voldoende. Maar ik weet niet wat er in die ellipj functie gebeurt. Ergens las ik dat de functiewaarde via een, weliswaar tiendegraads, polynoom wordt berekend. Je kunt je afvragen of die polynoom inderdaad alle kenmerken van de functie reproduceert, en onder welke voorwaarden.

Wat betreft Python, dat is inderdaad een fijne programmeertaal. Relatief eenvoudig aan te leren maar enorm krachtig.

[Professor Puntje]

@ Marko

Het zou fijn zijn als iemand die goed thuis is in Python en elliptische functies daar antwoord op geeft. Voorlopig heb ik voldoende vertrouwen in de nauwkeurigheid van Python om het resultaat van dit topic serieus te nemen, maar het bewijs is nu inderdaad (nog) niet waterdicht.

[OOOVincentOOO]

Hallo,

Ik heb de code in Jupyternotebook gezet en op Github geplaats. Via onderstaande link zou iedereen theoretisch de interactieve grafieken kunnen laten lopen via de webbrowser.

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

Het laden zou eventjes kunnen duren (in sommige gevallen krijg je een build dan kan het laden enkele minuten duren). Het programma loopt dan in de cloud. De code kun je starten via het driehoekje "Run". Als het goed is moet je twee keer op run klikken: de eerste cell bevat Latex en de tweede cell de Python code. Met de rammel computer op het werk gaat het goed.

Het laden van de code duurt enkele seconden. Dus niet gaan button bashen, maar zover ik weet is mybinder en Jupyter best robuust.

Met de slide bars kun je de instellingen veranderen: inzoomen, solar masses en radius afstand lichtstraal. Met de knopjes: y=x (zet de as verdeling gelijk, 'de echte verhouding'. En met knopje: sun kun de de afmetingen zon zien en de ingestelde straal.

Excuses als de code slecht en onoverzichtelijk geprogrammeerd is. Ik heb mijn best gedaan om het zo eenvoudig mogelijk te doen en commentaar te plaatsen.

Ik heb het getest met: edge, chrome en firefox. Het zou leuk zijn als het bij jullie ook werkt.

Zo moet het eruit zien:

[Professor Puntje]

Even snel getest, je moet inderdaad een tijdje wachten maar daarna lijkt het dan ook bij mij goed te werken. Eén commentaar puntje: in plaats van (6) werken we nu met (6*) vanwege de notatie kwestie.

Kan er op GitHub ook commentaar geplaatst worden, dat lijkt mij dan interessant om te lezen.

[OOOVincentOOO]

Wij kunnen dit document updaten (ik weet niet hoe dat te doen voor meerdere gebruikers). Er ook een revisie beheer. De orginele link blijft dan steeds werken dat is wel makkelijk.

[Professor Puntje]

Met commentaar bedoel ik meer van andere gebruikers van GitHub, niet van mijzelf. Wat die ervan vinden, en of ze er nog fouten in zien.

[OOOVincentOOO]

Hallo,

Ik ben niet goed thuis in Github. Maar ik heb gisteren gezien dat er ook een forum is. Voorlopig nog geen interesse te socialiseren op een ander forum. Ik heb nog enige details uitgewerkt.

In mijn versie van de simulatie heb ik de Newton-Rapson eruit gehaald betreffende de hoek bereik \(\varphi\): na (flink) puzzelwerk kwam ik erachter dat deze bepaald kan worden met de incomplete Jacobi Elliptic integraal. Zover ik kan zien lijkt het te werken. Ik begrijp het kleine stapjes meer erg mooie wiskunde Jacopi Elliptics. Veel te verkennen in dit vakgebied.

Ik heb de interactieve versie geupdate (oude link werkt ook). De verie is nog niet gebuild voor iedere browser. Dus indien je pech hebt word eerst de build geladen in de cloud:
https://mybinder.org/v2/gh/oooVincentoo ... ptic.ipynb

De onine versie loopt een beetje traag. Wellicht dat ik nog ga bekijken of het ook makkelijk kan in Geogebra maar daar heb ik minder ervaring mee. Volgens mij zijn er minder standaard functies voorhanden en dient men de Jacobi integraal zelf te berekenen.

Onze resulaten heb ik deze manier samengevat:
==================================================================================
Trajectory of Light Signal using Jacobi’s Elliptic Function.

This workbook describes the: Trajectory of Light Signal using Jacobi’s Elliptic Function. The light deflection arround the sun is determined.

The theory behind the method is found in paper: 'Equations of Orbits, Deflection of Light, Mercury’s Perihelion Shift, Period of Revolution and Time Delay of Signals in Schwarzchild’s Geometry. Some closed-form solutions using Weinberg’s approach. New Version (19/08/2016) Solomon M. Antoniou'.[Researchgate]

This project was started based upon a topic on: [WF] by: Professor Puntje, CoenCo and oooVincentooo.

Here follows a summary of the result found on forum Wetenschapsforum. The mean set of formulas required are found in article. Where \(r_s\) is Schwarzschild radius and \(r_0\) radius lightbeam from center mass and \(\varphi\) is the angle following the light signal. The radius \(r\) with angle \(\varphi\) are computed with \((6)\). Finally the startcondition \(\sigma\) and \(\varphi\) range are computed:

$$e_1 = \frac{r_0 - r_s + \sqrt{(r_0 - r_s)(r_0 + 3r_s)}}{2 r_s r_0} \tag{1}$$
$$e_2 = \frac{1}{r_0} \tag{2}$$
$$e_3 = \frac{r_0 - r_s - \sqrt{(r_0 - r_s)(r_0 + 3r_s)}}{2 r_s r_0} \tag{3}$$
$$\tau = \sqrt{ \frac{ r_s (e_1 - e_3)}{ 4 }} \tag{4}$$
$$h = \sqrt{ \frac{ e_2 - e_3}{ e_1 - e_3 }} \tag{5}$$
$$\boxed{ r = \frac{1}{ e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \varphi + \sigma ; h^2)}} \tag{6}$$

Start conditions for \(r=r_{0}\) at \(\varphi=\pi/2\) determine \(\sigma\) substitude \((2)\):
$$r_{0} = \frac{1}{ e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \frac{\pi}{2} + \sigma ; h^2)}=\frac{1}{e_2} \tag{7}$$
After simplification:
$$\mathrm{sn}( \tau \frac{\pi}{2} + \sigma ; h^2) = 1 \tag{8}$$
Where \(\tau\) and \(h^2\), are given and the startcondition \(\sigma\) can be solved with standard form complete elliptic integral [Wiki]. In the simulation the following integral solver was used from scipy: [special.ellipk]. Definition of Jacobi elliptic integral:
$$K(h^2)=\int_{0}^{\frac{\pi}{2}} \frac{d t}{\sqrt{1-h^2 \sin^2(t) }} =\tau \frac{\pi}{2} + \sigma \tag{9} $$
Property elliptic integrals:
$$\mathrm{sn}(K,h^2)=\sin \left( \frac{\pi}{2} \right) \tag{10} $$
Not all solutions are valid, the root of \((6)\) van have zeros. The \(x\)-range is found by searching the roots of \((6)\).
$$e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \varphi + \sigma ; h^2) = 0 \tag{12}$$
$$ \mathrm{sn}( \tau \varphi + \sigma ; h^2) =\sqrt{ -\frac{e_3}{e_2 - e_3}} \tag{13}$$
$$\mathrm{sn}(K,h^2)=\sin \left( \phi \right) \tag{14} $$
$$\phi=\arcsin \left( \sqrt{ -\frac{e_3}{e_2 - e_3}} \right) \tag{15}$$
Where \(\tau\), \(h^2\), \(\phi\) and \(\sigma\) are given and can be solved with incomplete elliptic integral [scipy: special.ellipkinc]. The maximum angle \(\varphi\) can be calulated with:
$$K(\phi,h^2)=\int_{0}^{\phi} \frac{d t}{\sqrt{1-h^2 \sin^2(t) }} =\tau \varphi + \sigma \tag{16} $$
==================================================================================

Voor Python console (niet getest voor verouderde versies):

Code: Selecteer alles

#Version 15.0 Schwarzschild, Jacobi Elliptic
#https://www.wetenschapsforum.nl/viewtopic.php?f=85&t=212427

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

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


#%matplotlib widget

fig= plt.figure(figsize=(9, 5))  
widths = [5, 5, 5]
heights = [5, 1]
spec = gridspec.GridSpec(ncols=3, nrows=2, figure=fig,width_ratios=widths, height_ratios=heights)

ax1 = fig.add_subplot(spec[0, 0])
ax2 = fig.add_subplot(spec[0, 1])
ax3 = fig.add_subplot(spec[0, 2])  


plt.suptitle('Schwarzschild, Jacobi Elliptic', fontsize=12)
#plt.tight_layout()
#spec.tight_layout(fig)
fig.set_tight_layout(True)

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

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

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

rax = plt.axes([0.75, 0.01, 0.075, 0.14])
check = CheckButtons(rax, ('x=y', '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)
    
    #Check if solution is valid, sr*Ro radius must be bigger rs
    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 x range for all valid 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)) +  "$''$\n"    
         )
    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.max(x)-range*2/100,np.max(x)]+range*2/100)
                ax2.set_xlim([-np.max(x)-range*2/100,np.max(x)]+range*2/100)
                ax3.set_xlim([-np.max(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.max(x)-range*2/100,np.max(x)]+range*2/100)
                ax2.set_xlim([-np.max(x)-range*2/100,np.max(x)]+range*2/100)
                ax3.set_xlim([-np.max(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)
                
        i=i+1

def clearsetplot():
    ax1.clear()
    ax1.set_title('Light Path', fontsize=10, pad=15)
    ax1.set_xlabel('x [m]', fontsize=10)
    ax1.set_ylabel('y [m]', fontsize=10)
    ax1.ticklabel_format(axis='both', style='sci', scilimits=(0,0))
    
    ax2.clear()
    ax2.set_title('Angular Deflection', fontsize=10, pad=15)
    ax2.set_xlabel('x [m]', fontsize=10)
    ax2.set_ylabel('deflection angle  [rad]', fontsize=10)
    ax2.ticklabel_format(axis='both', style='sci', scilimits=(0,0))
    
    ax3.clear()
    ax3.set_title('Angular Distribution', fontsize=10, pad=15)
    ax3.set_xlabel('x [m]', fontsize=10)
    ax3.set_ylabel('angular change [rad/m]', fontsize=10)
    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]

mmm, ik zie net nog een fout. Volgens mij is (16) verkeerd!

EDIT:
Het bereik via (16) is een kleine beetje groter bij extreme solar masses voor radius verandering 100% match hoek bereik. Voor de rest komt Newt-Rapson en (16) overeen.

Maar nog veel vragen, waarom. Maar lijkt dis toch te werken.

[OOOVincentOOO]

Ik was even in paniek.
Formule (16) lijkt allemaal te werken, de radius lichtstraal komt bij grote massa's in de buurt van \(r_0\). Ik vermoed dat Newton-Rapson dan problemen krijgt met deze kleine getallen.

Merk op dat een type extrema nog niet afgevangen is: indien de radius erg klein gekozen word waarbij: \(r_0<r_s\). Dit levert blanco grafieken op de Mikado effecten heb ik (nog) niet gezien bij deze versie.

Nu even topic laten rusten, het is weekeinde!

Interactief Jupyter Notebook:
https://mybinder.org/v2/gh/oooVincentoo ... ptic.ipynb

Versie 16:
Voor python console gebaseerd op laatste werkende versie. De voorgaande 15.0 geeft error als de code tweemaal gestart word.

Code: Selecteer alles

#Version 16.0 Schwarzschild, Jacobi Elliptic
#https://www.wetenschapsforum.nl/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=12)
#plt.tight_layout()
#spec.tight_layout(fig)
 

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

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

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

rax = plt.axes([0.75, 0.01, 0.075, 0.14])
check = CheckButtons(rax, ('x=y', '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)
    
    #Check if solution is valid, sr*Ro radius must be bigger rs
    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 x range for all valid solutions  
    
    #Find root angle for non-solutions
    root = optimize.newton(r,0.3)
    alphamax2 =np.abs(root) # beveiliging tegen singulariteiten
    print('Newton:' + str(alphamax2))
    
    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)) +  "$''$\n"    
         )
    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.max(x)-range*2/100,np.max(x)]+range*2/100)
                ax2.set_xlim([-np.max(x)-range*2/100,np.max(x)]+range*2/100)
                ax3.set_xlim([-np.max(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.max(x)-range*2/100,np.max(x)]+range*2/100)
                ax2.set_xlim([-np.max(x)-range*2/100,np.max(x)]+range*2/100)
                ax3.set_xlim([-np.max(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)
                
        i=i+1

def clearsetplot():
    ax1.clear()
    ax1.set_title('Light Path', fontsize=10, pad=15)
    ax1.set_xlabel('x [m]', fontsize=10)
    ax1.set_ylabel('y [m]', fontsize=10)
    ax1.ticklabel_format(axis='both', style='sci', scilimits=(0,0))
    
    ax2.clear()
    ax2.set_title('Angular Deflection', fontsize=10, pad=15)
    ax2.set_xlabel('x [m]', fontsize=10)
    ax2.set_ylabel('deflection angle  [rad]', fontsize=10)
    ax2.ticklabel_format(axis='both', style='sci', scilimits=(0,0))
    
    ax3.clear()
    ax3.set_title('Angular Distribution', fontsize=10, pad=15)
    ax3.set_xlabel('x [m]', fontsize=10)
    ax3.set_ylabel('angular change [rad/m]', fontsize=10)
    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]

r₀ < r_s = zwart gat. Dan bestaat er ook geen lichtstraal die zo'n zon rakelings kan passeren.

Je versie 16 geeft in mijn oude Spyder lege grafieken, en:

Code: Selecteer alles

    raise AttributeError('Unknown property %s' % k)

AttributeError: Unknown property pad

In Jupyter werkt het bij mij wel.