Hallo,
Mooi, straks jouw methode voor de alpha limiet bekijken. Kijken of ik het ga snappen. Ben meer aan het programmeren en tweaken tussendoor.
Mijn code heb nog aangepast. Hopelijk werkt het nu wel. Misschien heb ik onbewust routines uit Jupyter notebook overgenomen. Ik gebruik Spyder vanuit Anaconda installatie.
De zon heeft nu altijd de vaste diameter. Ik heb een kleine circel eromheen getekend met afstand lichtstraal. Dan is het hopelijk beter te begrijpen.
Deflectie hoek berekening (enkele hoek): Via numerieke diff. en de maximale waarde genomen.
Het lijkt erop dat de Mikado toestanden ook de goede hoek geven. Maar dien hier meer te testen.
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 scipy.special as sps
from matplotlib.widgets import Slider, CheckButtons
fig, ax= plt.subplots(2,1,figsize=(9, 5),gridspec_kw={'height_ratios': [5, 1]})
#Do not display second dummy graph (position widgets)
ax[1].axis('off')
# Draw slider alpha
axcolor = 'lightgoldenrodyellow'
axalpha = plt.axes([0.25, 0.05,0.4, 0.04], facecolor=axcolor)
s_alpha = Slider(axalpha, 'zoom: ', 0, 1.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: ', 1, 10, 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
#The following function updateplot runs whenever slider changes value
def updateplot(val):
ax[0].clear()
ax[0].set_title('Light Deflection (Jacobi Elliptic)', fontsize=20)
ax[0].set_xlabel('x', fontsize=14)
ax[0].set_ylabel('y', fontsize=14)
ax[0].autoscale(enable=False)
#Solar masses/radius and magnification alpha
alpha = s_alpha.val
sm = s_sm.val
sr = s_sr.val
phi=np.linspace(np.pi/2-alpha*np.pi/2,np.pi/2+alpha*np.pi/2,1000)
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)
#print('e1: ' + str(e1))
e2=1/(sr*Ro)
#print('e2: ' + str(e2))
e3=(sr*Ro-sm*2*M*G-esqrt)/(sm*4*M*G*sr*Ro)
#print('e3: ' + str(e3))
Tau=np.sqrt(sm*M*G*(e1-e3)/2)
#print('Tau: ' + str(Tau))
h=(e2-e3)/(e1-e3)
#print('h: ' + str(h))
#Solve angle start conditions
angle=sps.ellipk(h)
#print('angle: ' + str(angle))
#Calculate Sigma Start Condition
Sigma=-Tau*np.pi/2-angle
#print('sigma: ' + str(Sigma))
v=Tau*phi+Sigma
sn,_,_,_=sps.ellipj(v,h)
r=1/(e3+(e2-e3)*(sn**2))
x=r*np.cos(phi)
y=r*np.sin(phi)
#plot curve
ax[0].plot(x,y)
ax[0].grid()
#Mean Deflection Angle
dx=np.diff(x)
dy=np.diff(y)
da=np.arctan(dy/dx)
da=np.abs(da)
dam=np.degrees(np.max(da))
text=('Deflection:\n' +
str('{:08.6f}'.format(dam)) + '$^{\circ}$\n' +
str('{:08.2f}'.format(dam*3600)) + "$''$")
ax[0].text(0.01, 0.75,text, transform=ax[0].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:
ax[0].axis('auto')
range=np.max(y)-np.min(y)
ax[0].set_ylim([np.min(y)-range*2/100,np.max(y)+range*2/100])
range=np.max(x)-np.min(x)
ax[0].set_xlim([np.min(x)-range*2/100,np.max(x)]+range*2/100)
else:
ax[0].axis('equal')
range=np.max(x)-np.min(x)
ax[0].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),sr*Ro,fill=False, ls='--',ec='black')
ax[0].add_patch(circle)
circle=plt.Circle((0,np.max(y)-sr*Ro),Ro, color='orange')
ax[0].add_patch(circle)
i=i+1
#Display plot n startup.
updateplot(1)
#The following code checkes if slider changes. This line is looped automatic by pyplot
s_alpha.on_changed(updateplot);
s_sm.on_changed(updateplot);
s_sr.on_changed(updateplot);
check.on_clicked(updateplot);
jij zult ook veel Mikado oplossingen gezien hebben! 