Zou je kunnen uitleggen hoe je komt tot:
dφdx=12∂∂Rln((1−rsr)2rsrx2r2+1−rsr)
- Hoe en waarom kom je tot de ln?
- Als je mathpages en andere litt. zie je dat c(r) de component is in deze diff.
- Volgens mijn geheugen heb ik immer een wortel gezien.
Als ik numeriek reken:
Dus volgens jouw afleiding is dan:
c(r)=12ln((1−rsr)2rsrx2r2+1−rsr)
Observatie: verkeerde deflectie hoek (3.5") (gelijk aan jouw eindformule zie enkele berichten terug). Factor twee verschil. Op y=2 (dwz. 2R is het correct).
Volgens mijn inzicht moet het dit zijn:
c(r)=12√((1−rsr)2rsrx2r2+1−rsr)
Observatie: correcte deflectie hoek (1.75")
Volens jouw eindformule en de formule voor
c(r) krijg ik de verkeerde afbuighoek. Methode in deze reply en hier (met jouw eindformule):
OOOVincentOOO schreef: ↑vr 17 sep 2021, 10:32
Mijn analyse jouw aanpak:
Met de eindformule en de formule voor
c(r). Factor 2 te veel. Zie ook edit hieronder
Discussie:
Waarom heb geen recht iets te willen begrijpen? Op twee verschillende manieren toon ik een foute afbuighoek aan volgens jouw formules.
Ik kan ook een fout maken in mijn berekening, wellicht kan jij zelf de deflektie hoek uitrekenen. Dan ben je niet afhankelijk van mijn resultaten.
Volgens mij voeren wij degelijke wetenschap uit. We zijn momenteel jouw afleiding aan het toetsen op correctheid.
Fouten maken bestaat niet in wetenschap. Als ik al mijn klad papier zou posten zou je schrikken! Maar dat verbaasd jouw misschien helemaal niet.
EDIT:
Indien ik dit doe:
c(r)=14ln((1−rsr)2rsrx2r2+1−rsr)
Krijg ik ook het goede resultaat.
Code: Selecteer alles
#Open pyplot in separate interactive window
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'qt5')
#https://www.mathpages.com/rr/s8-09/8-09.htm
#https://www.mathpages.com/rr/s6-03/6-03.htm
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['axes.titlepad'] = 20
widths = [5,5,5]
heights = [5,5]
fig= plt.figure(figsize=(30,10))
gs=fig.add_gridspec(2,3,width_ratios=widths, height_ratios=heights)
ax1=fig.add_subplot(gs[0,0])
ax2=fig.add_subplot(gs[0,1])
ax3=fig.add_subplot(gs[0,2])
ax4=fig.add_subplot(gs[1,0])
ax5=fig.add_subplot(gs[1,1])
ax6=fig.add_subplot(gs[1,2])
#Physical constants
M=1.989e30
G=6.67408e-11
c=3e8
R=696340000
#schwarzschild Radius
Rs=4*M*G/c**2
#Function c(r)
def cxy(x,y):
#puntje met log
cr=0.5*np.log((1-Rs/np.sqrt(x**2+y**2) )**2/(Rs/np.sqrt(x**2+y**2) * x**2/(x**2+y**2) +1 - Rs/np.sqrt(x**2+y**2) ))
#volgens mathpages en litteratuur
cr=0.5*np.sqrt((1-Rs/np.sqrt(x**2+y**2) )**2/(Rs/np.sqrt(x**2+y**2) * x**2/(x**2+y**2) +1 - Rs/np.sqrt(x**2+y**2) ))
return cr
size=2500
#Create 2D mesh and calculate c(r)
x=np.linspace(-30,30,size)
x=R*x
y=np.linspace(1,6,size)
y=R*y
X,Y =np.meshgrid(x,y)
z=cxy(X,Y)
#Function c(r) as function of x and y
cf1=ax1.contourf(x/R,y/R,z,levels=15, cmap='seismic',alpha=1)
fig.colorbar(cf1,ax=ax1, location="bottom")
ax1.set_title(r'$c(\sqrt{x^2+y^2})$',fontsize=20)
ax1.set_xlabel(r'$x$',fontsize=15)
ax1.set_ylabel(r'$y$',fontsize=15)
#Gradient for d/dy this is dTheta/dx
dcdy = np.gradient(z, axis=0)/(y[5]-y[4])
cf2=ax2.contourf(x/R,y/R,dcdy,levels=15, cmap='seismic',alpha=1)
fig.colorbar(cf2,ax=ax2, location="bottom")
ax2.set_title(r'$\partial c / \partial y = d \Theta /dx$',fontsize=20)
ax2.set_xlabel(r'$x$',fontsize=15)
ax2.set_ylabel(r'$y$',fontsize=15)
#Integrate Theta over z
intz=np.cumsum(dcdy,axis=1)*(x[5]-x[4])
cf3=ax3.contourf(x/R,y/R,np.degrees(intz)*3600,levels=15, cmap='seismic',alpha=1)
fig.colorbar(cf3,ax=ax3, location="bottom")
ax3.set_title(r'$\int \Theta \ dx$',fontsize=20)
ax3.set_xlabel(r'$x$',fontsize=15)
ax3.set_ylabel(r'$y$',fontsize=15)
#Function c(r) as function of x and y
for i in range(np.size(y)):
if 5*(i/size)%1==0:
cr=z[i]
ax4.plot(x/R,cr,color="black", linewidth=0.5, label="x=" + str(x[i]))
ax4.annotate("y=" + str(round(y[i]/R,2)), xy=( x[np.argmin(cr)]/R , np.min(cr)),ha='center', xycoords='data', xytext=(-100 , -5), textcoords='offset points',rotation = 0,arrowprops=dict(arrowstyle='-',lw=0.2))
ax4.set_xlabel(r'$x$',fontsize=15)
ax4.set_ylabel(r'$c(\sqrt{x^2+y^2})$',fontsize=15)
ax4.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
#Gradient for d/dy this is dTheta/dx
p=0
for i in range(np.size(y)):
if 5*(i/size)%1==0:
dphidy=dcdy[i]
ax5.plot(x/R,dphidy,color="black", linewidth=0.5, label="x=" + str(x[i]))
ax5.annotate("y=" + str(round(y[i]/R,2)), xy=( x[np.argmax(dphidy)]/R , np.max(dphidy)),ha='center', xycoords='data', xytext=(((-1)**p)*100 ,-5), textcoords='offset points',rotation = 0,arrowprops=dict(arrowstyle='-',lw=0.2))
ax5.set_xlabel(r'$x$',fontsize=15)
ax5.set_ylabel(r'$d \Theta /dx$',fontsize=15)
p+=1
print(p)
#Integrate Theta over z
for i in range(np.size(y)):
if 5*(i/size)%1==0:
phidx=np.degrees(intz[i])*3600
ax6.plot(x/R,phidx,color="black", linewidth=0.5)
ax6.annotate("y=" + str(round(y[i]/R,2)), xy=( x[np.argmax(phidx)]/R , np.max(phidx)),ha='center', xycoords='data', xytext=(-50 ,-10), textcoords='offset points',rotation = 0,arrowprops=dict(arrowstyle='-',lw=0.2))
ax6.set_xlabel(r'$x$',fontsize=15)
ax6.set_ylabel(r'$\int \Theta \ dx$ [arcsec.]',fontsize=15)
theta=4*G*M/(R*c**2)
theta=np.degrees(theta)*3600
ax6.plot([x[0]/R,x[-1]/R],[theta,theta],color="red", linewidth=0.5, label=r"$\dfrac{4GM}{Rc^2}$")
ax6.legend(loc="lower right")