Python en Jacobi

[Professor Puntje]

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.

[OOOVincentOOO]

Hierbij enige nettere grafieken. Zou ji de code bj jouw kunnen testen. Ik heb wederom een "tight fit" moeten toepassen. Ik gebruik nu een andere methode dan voorheen. Hoepelijk werkt deye wel bij jouw.

Ik heb het aantal phi punten tot 500 gereduceerd om beter vloeiend dynamisch te zijn.

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 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, figure=fig,width_ratios=widths, height_ratios=heights)

spec.tight_layout(fig)

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

plt.suptitle('Schwarzshield, 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, 'sun 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);

[Professor Puntje]

Dat geeft bij mij een lege grafiek. Foutcode:

Code: Selecteer alles

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

[OOOVincentOOO]

Probeer de onderstaande eens, hier heb ik tight fit gedisabled: "#spec.tight_layout(fig)".

Hier een extreme oplossing met half rondje om zwaar object:

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 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, figure=fig,width_ratios=widths, height_ratios=heights)

#spec.tight_layout(fig)

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, 'sun radius: ', 0.1, 2, 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);

[Professor Puntje]

Geeft weer een lege grafiek, en:

Code: Selecteer alles

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

[OOOVincentOOO]

Misschien heb jij een oude versie mathplot lib. Omdat jij mint geinstalleerd hebt (niet veel gedownload en geactualiseerd denk ik).

Code: Selecteer alles

import matplotlib
print(matplotlib.__version__)

Mijn versie is: 3.1.3

Misschien een update doen volgens mij kan je de volgende regel plakken in console Spyder Python type (wel alle pyplot vensters sluiten):

Code: Selecteer alles

pip install --upgrade matplotlib

https://stackoverflow.com/a/53454820

[Professor Puntje]

Code: Selecteer alles

import matplotlib
print(matplotlib.__version__)

2.1.1

[OOOVincentOOO]

Dan verschillen de versies nogal. Heb je geprobeerd:

Code: Selecteer alles

pip install --upgrade matplotlib

In spyder copieren en plakken in console venster rechts onder.

@CoenCo,
Misschien kan jij beter adviseren. Is er een instructie om Spyder geheel te actualiseren tot laatste versie. Ik weet niet wat Puntje geinstalleerd heeft? Spyder standalone?

[Professor Puntje]

In ben in Linux Mint wat huiverig om zaken buiten de repositories om te updaten of upgraden. Dat geeft vaak meer problemen dan voordelen. Mijn huidige spyder komt uit de repositories.

[Professor Puntje]

Nog een vraag over de code, wat betekenen de [:-1] en [:-2]?

[OOOVincentOOO]

Misschien dat @CoenCo kan helpen. Ik probeer de versie intussen backward compatible te maken.

a is een array:

Output element uit array b:
b=a[2] derde element uit array
b=a[n] n'de element
b=a[-1] laatste element in array

Ouput array's c:
c=a[2:6] elementen 2 tot 6
c=a[2:] alle elementen van 2 en verder
c=a[:-1] alle elementen behalve laatste
c=a[:-6] alle elementen behalve laatste 6

[Professor Puntje]

Is dat om de arrays voor deling van gelijke lengte te maken?

[OOOVincentOOO]

Je kunt niet een plot maken met twee verschillende lengten voor x en y.

Door np.diff te nemen neem je automatisch het verschil tussen elementen a[n+1]-a[n] de output array is een element kleiner.

Dus: dy=np.diff(y) is een element kleiner kan x.

[CoenCo]

Er is op een of andere manier python op de pc van proffessor puntje geinstalleerd.
Ik vermoed via anaconda, maar het kan ook een andere package manager zijn.
Via dia package manager moet het mogelijk zijn om te updaten, die zal dan ook alle dependencies checken.

[OOOVincentOOO]

@CoenCo,

Volgens mij lukte het niet Anoconda op mint PC van puntje te krijgen.

Edit:
Volgens mij mag je gewoon onderstaande doen. Ik heb dat zojuist ook gedaan op windows PC. Volgens mij zoekt Python zelf de directory uit. Het commando word uitgevoerd in Spyder en niet in Linux. toch?

Dan verschillen de versies nogal. Heb je geprobeerd:

Code: Selecteer alles

pip install --upgrade matplotlib

In spyder copieren en plakken in console venster rechts onder.

@CoenCo,
Misschien kan jij beter adviseren. Is er een instructie om Spyder geheel te actualiseren tot laatste versie. Ik weet niet wat Puntje geinstalleerd heeft? Spyder standalone?