sciencetalk


https://sciencetalk.nl/forum/

Python en Jacobi

https://sciencetalk.nl/forum/viewtopic.php?t=212427

11 van 14

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 12:55
door Professor Puntje
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.

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 13:04
door sensor
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.

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 13:04
door OOOVincentOOO
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);

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 13:15
door Professor Puntje
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.

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 13:16
door OOOVincentOOO
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.

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 13:18
door Professor Puntje
Heb ik net getest, maar die geeft bij mij ook lege grafiekjes.

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 13:36
door OOOVincentOOO
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);

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 13:47
door Professor Puntje
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.

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 13:55
door OOOVincentOOO
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);

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 14:00
door Professor Puntje
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.

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 14:02
door OOOVincentOOO
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

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 14:08
door Professor Puntje
Ja - die werken maar je krijgt soms onzinnige mikado oplossingen.

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 14:16
door Marko
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.

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 14:32
door Professor Puntje
Graag wat specifieker. Aan vage suggesties hebben we niets.

Re: Python en Jacobi

Geplaatst: do 10 jun 2021, 17:01
door OOOVincentOOO
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);
Alle tijden zijn UTC+01:00
11 van 14

Powered by phpBB® Forum Software © phpBB Limited

Nederlandse vertaling door phpBB.nl.