Python en Jacobi
- Berichten: 7.463
Re: Python en Jacobi
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.
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.
- Berichten: 1.605
Re: Python en Jacobi
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.
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);
- Berichten: 7.463
Re: Python en Jacobi
Dat geeft bij mij een lege grafiek. Foutcode:
Code: Selecteer alles
TypeError: __init__() got an unexpected keyword argument 'figure'
- Berichten: 1.605
Re: Python en Jacobi
Probeer de onderstaande eens, hier heb ik tight fit gedisabled: "#spec.tight_layout(fig)".
Hier een extreme oplossing met half rondje om zwaar object:
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);
- Berichten: 7.463
Re: Python en Jacobi
Geeft weer een lege grafiek, en:
Code: Selecteer alles
TypeError: __init__() got an unexpected keyword argument 'figure'
- Berichten: 1.605
Re: Python en Jacobi
Misschien heb jij een oude versie mathplot lib. Omdat jij mint geinstalleerd hebt (niet veel gedownload en geactualiseerd denk ik).
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):
https://stackoverflow.com/a/53454820
Code: Selecteer alles
import matplotlib
print(matplotlib.__version__)
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
- Berichten: 7.463
Re: Python en Jacobi
Code: Selecteer alles
import matplotlib
print(matplotlib.__version__)
2.1.1
- Berichten: 1.605
Re: Python en Jacobi
Dan verschillen de versies nogal. Heb je geprobeerd:
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?
Code: Selecteer alles
pip install --upgrade matplotlib
@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?
Laatst gewijzigd door OOOVincentOOO op do 10 jun 2021, 12:33, 1 keer totaal gewijzigd.
- Berichten: 7.463
Re: Python en Jacobi
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.
Laatst gewijzigd door Professor Puntje op do 10 jun 2021, 12:37, 1 keer totaal gewijzigd.
- Berichten: 7.463
Re: Python en Jacobi
Nog een vraag over de code, wat betekenen de [:-1] en [:-2]?
- Berichten: 1.605
Re: Python en Jacobi
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
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
Laatst gewijzigd door OOOVincentOOO op do 10 jun 2021, 12:40, 1 keer totaal gewijzigd.
- Berichten: 7.463
Re: Python en Jacobi
Is dat om de arrays voor deling van gelijke lengte te maken?
- Berichten: 1.605
Re: Python en Jacobi
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.
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.
-
- Technicus
- Berichten: 1.163
Re: Python en Jacobi
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.
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.
- Berichten: 1.605
Re: Python en Jacobi
@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?
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?
OOOVincentOOO schreef: ↑do 10 jun 2021, 12:29 Dan verschillen de versies nogal. Heb je geprobeerd:
In spyder copieren en plakken in console venster rechts onder.Code: Selecteer alles
pip install --upgrade matplotlib
@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?