• Over
  • Forummenu
    • Reageer
    Gebruikersavatar
    Professor Puntje
    Artikelen: 0
    Berichten: 7.592
    Lid geworden op: vr 23 okt 2015, 23:02

    Re: Python en Jacobi

    Installeren gaat op Linux Mint via de Software Manager. Dat ziet er zo uit:
    spyder
    Daar kan ik dan uit kiezen. Vaak is zijn dat oudere versies, maar het installeren van nieuwere versies is naar mijn ervaring een hachelijke zaak. Dat gaat wel maar leidt vaak tot problemen met andere programma's of dependencies. Ik ben er weken mee zoet geweest om de zo ontstane schade weer te herstellen. Dan maar liever niet het nieuwste van het nieuwste.
    sensor
    Artikelen: 0
    Berichten: 338
    Lid geworden op: vr 27 jan 2012, 11:42

    Re: Python en Jacobi

    Werkt het python project in een virtual environment bijvoorbeeld met venv of virtualenv. Op die manier kun je eigen veilige python omgeving maken. Op et moment dat er iets verkeerd geïnstalleerd staat dan kun je desbetreffende environment verwijderen zonder dat de rest van het systeem vastloopt. Overigens werkt spyder niet met virtualenv. Ik gebruik als IDE vscode.
    Gebruikersavatar
    OOOVincentOOO
    Artikelen: 0
    Berichten: 1.645
    Lid geworden op: ma 29 dec 2014, 14:34

    Re: Python en Jacobi

    mmm, probeer deze versie eens:

    Code: Selecteer alles

    #Version 6.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 matplotlib.gridspec as gridspec
import scipy.special as sps
from scipy import optimize
from matplotlib.widgets import Slider, CheckButtons                    

fig= plt.figure(figsize=(15, 6))  
widths = [5, 5, 5]
heights = [5, 1]
spec = gridspec.GridSpec(ncols=3, nrows=2,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, Jabobi 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.01, 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, '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=7e8   

def e1(sm,sr):
    if (sr*Ro>sm*2*M*G):
        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)
    else:
        e1=0.1
    return e1

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

def e3(sm,sr):
    if (sr*Ro>sm*2*M*G):
        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)
    else:
        e3=0
    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

   
    
#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
    root = optimize.newton(r,0,disp=False)
    alphamax =np.abs(root) # beveiliging tegen singulariteiten

    #Set angle range and calculate radius
    phi = np.linspace(np.pi/2-zoom*np.pi/2,np.pi/2+zoom*np.pi/2,500)
    #mask=(phi<np.pi-alphamax) | (phi>alphamax)
    #phi=phi[mask]
    #print(alphamax)
    if (alphamax<4*np.pi) and (alphamax<np.max(phi)):
    
        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('{:08.6f}'.format(dam))      +   ' $[rad]$\n' +
                str('{:08.6f}'.format(np.degrees(dam)))      +   '$^{\circ}$\n' +
                str('{:08.2f}'.format(np.degrees(dam)*3600)) +  "$''$")
        ax1.text(0.01, 0.8,text, transform=ax1.transAxes, fontsize=14)
        
        #scaling method and draw sun
        plotoptions(x,y,sr)
    
    else:
         #When array is empty disply 'no solution'
         text='No solution'
         ax1.text(0.5, 0.5,text, transform=ax1.transAxes, fontsize=20,color='red',ha='center')       
         x=[0,0]
         y=[0,0]
         
    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')
                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_zoom.on_changed(updateplot);
s_sm.on_changed(updateplot);
s_sr.on_changed(updateplot);
check.on_clicked(updateplot);
    Gebruikersavatar
    Professor Puntje
    Artikelen: 0
    Berichten: 7.592
    Lid geworden op: vr 23 okt 2015, 23:02

    Re: Python en Jacobi

    Ik heb nu ook al werkende Python code:

    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


fig, ax= plt.subplots(3,1,figsize=(4, 6))

# 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)
t = np.linspace(np.pi/2-1.4,np.pi/2+1.4,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)

ax[0].plot(x,y,c=(1,0.2,0.5),lw=1)
ax[0].title.set_text('lichtbaan')

dx=np.diff(x)
dy=np.diff(y)
dydx=dy/dx
afb=-np.arctan(dydx)

ax[1].plot(x[:-1],afb,c=(1,0.2,0.5),lw=1)
ax[1].title.set_text('lichtafbuiging')

dafb=np.diff(afb)
aanw=dafb/dx[:-1]

ax[2].plot(x[:-2],aanw,c=(1,0.2,0.5),lw=1)
ax[2].title.set_text('aanwas afbuiging')


plt.tight_layout()
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)

    Dat geeft de plots:

    plots

    Tenzij er nog fouten in onze aanpak gevonden worden, beschouw ik dit topic daarmee als succesvol ten einde gebracht. Een recentere versie van Spyder lijkt mij daarvoor niet nodig.
    Gebruikersavatar
    OOOVincentOOO
    Artikelen: 0
    Berichten: 1.645
    Lid geworden op: ma 29 dec 2014, 14:34

    Re: Python en Jacobi

    Heb je mijn laatste code getest? Lijkt mij wel zo leuk.

    Edit:
    Indien je geen python geinstalleerd hebt kun je veel code laten open op:
    https://trinket.io/embed/python3
    Echter mijn versie niet.
    Laatst gewijzigd door OOOVincentOOO op do 10 jun 2021, 13:19, 2 keer totaal gewijzigd.
    Gebruikersavatar
    Professor Puntje
    Artikelen: 0
    Berichten: 7.592
    Lid geworden op: vr 23 okt 2015, 23:02

    Re: Python en Jacobi

    Heb ik net getest, maar die geeft bij mij ook lege grafiekjes.
    Gebruikersavatar
    OOOVincentOOO
    Artikelen: 0
    Berichten: 1.645
    Lid geworden op: ma 29 dec 2014, 14:34

    Re: Python en Jacobi

    Even een andere methode geprobeerd. En deze versie lucky 7.0?
    Figure_1 v7

    Code: Selecteer alles

    #Version 7.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, Jabobi 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.01, 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, '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=7e8   

def e1(sm,sr):
    if (sr*Ro>sm*2*M*G):
        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)
    else:
        e1=0.1
    return e1

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

def e3(sm,sr):
    if (sr*Ro>sm*2*M*G):
        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)
    else:
        e3=0
    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

   
    
#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
    root = optimize.newton(r,0,disp=False)
    alphamax =np.abs(root) # beveiliging tegen singulariteiten

    #Set angle range and calculate radius
    phi = np.linspace(np.pi/2-zoom*np.pi/2,np.pi/2+zoom*np.pi/2,500)
    #mask=(phi<np.pi-alphamax) | (phi>alphamax)
    #phi=phi[mask]
    #print(alphamax)
    if (alphamax<4*np.pi) and (alphamax<np.max(phi)):
    
        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('{:08.6f}'.format(dam))      +   ' $[rad]$\n' +
                str('{:08.6f}'.format(np.degrees(dam)))      +   '$^{\circ}$\n' +
                str('{:08.2f}'.format(np.degrees(dam)*3600)) +  "$''$")
        ax1.text(0.01, 0.75,text, transform=ax1.transAxes, fontsize=14)
        
        #scaling method and draw sun
        plotoptions(x,y,sr)
    
    else:
         #When array is empty disply 'no solution'
         text='No solution'
         ax1.text(0.5, 0.5,text, transform=ax1.transAxes, fontsize=20,color='red',ha='center')       
         x=[0,0]
         y=[0,0]
         
    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')
                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_zoom.on_changed(updateplot);
s_sm.on_changed(updateplot);
s_sr.on_changed(updateplot);
check.on_clicked(updateplot);
    Gebruikersavatar
    Professor Puntje
    Artikelen: 0
    Berichten: 7.592
    Lid geworden op: vr 23 okt 2015, 23:02

    Re: Python en Jacobi

    Nog steeds lege grafiekjes, met de foutcode:

    Code: Selecteer alles

        root = optimize.newton(r,0,disp=False)

TypeError: newton() got an unexpected keyword argument 'disp'
    Als je "disp=False" weglaat krijg je wel curves.
    Gebruikersavatar
    OOOVincentOOO
    Artikelen: 0
    Berichten: 1.645
    Lid geworden op: ma 29 dec 2014, 14:34

    Re: Python en Jacobi

    Dat stukje disp is om lege oplossingen op te vangen. Wellicht ook niet beschikbaar in jouw versie. Volgende keer beter backward compatible programmeren, goede les voor mij.

    EDIT:
    Ook interatief? Met sliders?

    Heb je uberhaubt code geprobeerd van mij in de tussentijd? Want die disp staat al lange tijd in mijn code. Ik krijg het idee dat je mijn code's helemaal niet hebt laten lopen? Ik bestudeer ook jouw bijdragen.

    Hier een aangepaste versie:

    Code: Selecteer alles

    #Version 8.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, Jabobi 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.01, 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, '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=7e8   

def e1(sm,sr):
    if (sr*Ro>sm*2*M*G):
        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)
    else:
        e1=0.1
    return e1

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

def e3(sm,sr):
    if (sr*Ro>sm*2*M*G):
        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)
    else:
        e3=0
    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

   
    
#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
    root = optimize.newton(r,0)
    alphamax =np.abs(root) # beveiliging tegen singulariteiten

    #Set angle range and calculate radius
    phi = np.linspace(np.pi/2-zoom*np.pi/2,np.pi/2+zoom*np.pi/2,500)
    #mask=(phi<np.pi-alphamax) | (phi>alphamax)
    #phi=phi[mask]
    #print(alphamax)
    if (alphamax<4*np.pi) and (alphamax<np.max(phi)):
    
        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('{:08.6f}'.format(dam))      +   ' $[rad]$\n' +
                str('{:08.6f}'.format(np.degrees(dam)))      +   '$^{\circ}$\n' +
                str('{:08.2f}'.format(np.degrees(dam)*3600)) +  "$''$")
        ax1.text(0.01, 0.75,text, transform=ax1.transAxes, fontsize=14)
        
        #scaling method and draw sun
        plotoptions(x,y,sr)
    
    else:
         #When array is empty disply 'no solution'
         text='No solution'
         ax1.text(0.5, 0.5,text, transform=ax1.transAxes, fontsize=20,color='red',ha='center')       
         x=[0,0]
         y=[0,0]
         
    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')
                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_zoom.on_changed(updateplot);
s_sm.on_changed(updateplot);
s_sr.on_changed(updateplot);
check.on_clicked(updateplot);
    Gebruikersavatar
    Professor Puntje
    Artikelen: 0
    Berichten: 7.592
    Lid geworden op: vr 23 okt 2015, 23:02

    Re: Python en Jacobi

    Ja - ik heb steeds jouw codes uitgeprobeerd. Maar zoals ik al schreef ging het mij in dit topic om de pieken, en die vraag is nu opgelost. Wat ik nu nog uitprobeer doe ik enkel om jou een plezier te doen.
    Gebruikersavatar
    OOOVincentOOO
    Artikelen: 0
    Berichten: 1.645
    Lid geworden op: ma 29 dec 2014, 14:34

    Re: Python en Jacobi

    Mag ik U dan vragen of de sliders ook werken? Dat zou ik wel zo prettig vinden.

    EDIT:
    Tevens kan met het geprogrammeerde tooltje de extremen verkend worden (zoals de 360 lus). Zonder elkaars bijdragen was er niets tot stand gekomen!

    Dit topic: Python en Jacobi

    Gr,

    Vince
    Gebruikersavatar
    Professor Puntje
    Artikelen: 0
    Berichten: 7.592
    Lid geworden op: vr 23 okt 2015, 23:02

    Re: Python en Jacobi

    Ja - die werken maar je krijgt soms onzinnige mikado oplossingen.
    Gebruikersavatar
    Marko
    Artikelen: 0
    Berichten: 10.605
    Lid geworden op: vr 03 nov 2006, 23:08

    Re: Python en Jacobi

    Professor Puntje schreef: ↑do 10 jun 2021, 11:16 Mooi - moet nog even nalopen wat je precies doet, maar dat ziet er slecht uit voor de twee, of drie of meer pieken. Het is waarschijnlijk gewoon netjes maar één enkele piek in het midden.

    Als dit klopt dan zijn de twee pieken van MathPages geen gevolg van het gebruik van de xy-coördinaten want die hebben wij hier zelf ook gebruikt. De enige andere oorzaak die ik dan nog kan bedenken zijn de op MathPages toegepaste benaderingen. Die twee pieken zijn dan artefacten en geen fysisch reële verschijnselen, zij behoeven dus ook geen fysische verklaring.
    Of de ene piek die je hier vindt is een artefact van de benaderingen die hier worden toegepast.
    Cetero censeo Senseo non esse bibendum
    Gebruikersavatar
    Professor Puntje
    Artikelen: 0
    Berichten: 7.592
    Lid geworden op: vr 23 okt 2015, 23:02

    Re: Python en Jacobi

    Graag wat specifieker. Aan vage suggesties hebben we niets.
    Gebruikersavatar
    OOOVincentOOO
    Artikelen: 0
    Berichten: 1.645
    Lid geworden op: ma 29 dec 2014, 14:34

    Re: Python en Jacobi

    Hier nog een update versie 10.0.

    De limieten waartussen de oplossing bestaan is hier bepaald:
    Professor Puntje schreef: ↑wo 09 jun 2021, 12:57 Hoe los je φ uit onderstaande vergelijking voor gegeven waarden e2, e3, τ, σ en h in Python numeriek op?
    \( e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \varphi + \sigma ; h^2) = 0 \)
    Deze code was al geintegreerd door ons beide. Ik heb nu een optimalisatie gemaakt in de code. De Mikado oplossingen komen nu niet voor als het goed is. Voorheen kon een groter bereik ingesteld worden.

    Hier twee oplossingen voor extreme massa's met een kleine straal. In een geval sterk vervormde verdeling. Of dit fysisch mogelijk is weet ik niet.

    Werkt deze code bij jouw Puntje?
    Version 10
    10

    Code: Selecteer alles

    #Version 10.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, Jabobi 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=7e8   

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

#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
    root = optimize.newton(r,0)
    alphamax =np.abs(root) # beveiliging tegen singulariteiten

    #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('{:08.6f}'.format(dam))      +   ' $[rad]$\n' +
            str('{:08.6f}'.format(np.degrees(dam)))      +   '$^{\circ}$\n' +
            str('{:08.2f}'.format(np.degrees(dam)*3600)) +  "$''$")
    ax1.text(0.01, 0.75,text, transform=ax1.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:
                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);
    Reageer

    Terug naar “Informatica en programmeren”

    Sciencetalk

    Hét onafhankelijke platform voor het delen van wetenschappelijke kennis en informatie

    • info@sciencetalk.nl
    • 0648478721