Dat integreren met de quad-functie werkt daar duidelijk niet en ik heb geen idee waarom niet.
Ik heb nu zelf een Simpson-integratie gedaan. Ik doe dat over twee verschillende deeltrajecten en tel er de restterm (integratie van \(\lambda^{-1,5}\) naar oneindig) nog bij op. Dat lijkt goed te werken.
Verder heb ik de aarde opgedeeld in een sferische kern met straal 3450 km met ρ=11000 en een mantel (ellipsoïde) met ρ=4500. Die straal komt uit een bekend plaatje met een grafiek van de dichtheid van de aarde link, de twee dichtheden heb ik met een natte vinger uit datzelfde plaatje geschat. Verder natuurlijk gekeken of ik goed in de buurt kwam bij de bekende potentialen aan het aardoppervlak. Niet perfect, maar een beter benadering dan met een uniforme dichtheid.
Voor de duidelijkheid: ik bereken de potentiaal voor de hele ellipsoide met de dichtheid van de mantel rhom en tel daar de potentiaal van een kern met met een dichtheid gelijk aan het verschil tussen de echte dichtheid van de kern en de mantel, rhok-rhom, bij op.
Voor de evenaar geeft dat 62701127 J/kg, voor de pool 62809006 J/kg.
Het verschil is 107878 J/kg, terwijl de potentaal door de rotatie \(\frac{1}{2}re^2\omega^2\) 106849 J/kg is.
Dat het zo mooi overeen komt is waarschijnlijk toeval, maar wel leuk.
Maar ook bij uniforme dichtheid (2 regels uncommenten) kom je op een verschil van 84082 J/kg, ook niet ver van die rotatiepotentiaal.
Code: Selecteer alles
import numpy as np
import scipy
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
rp = 6356752
re = 6378137
m = 5.97e24
G = 6.67e-11
w = 7.2722e-5
a0 = re
b0 = re
c0 = rp
rk=3450000
rhok=11000
rhom=4500
# uniforme dichtheid
# rhok=5520
# rhom=rhok
def fun(u):
return (x ** 2 / (a ** 2 + u)) + (y ** 2 / (b ** 2 + u)) + (z ** 2 / (c ** 2 + u)) - 1
def integrand(lam):
return (1 - (x ** 2 / (a ** 2 + lam)) - (y ** 2 / (b ** 2 + lam)) - (z ** 2 / (c ** 2 + lam))) / (
((a ** 2 + lam) * (b ** 2 + lam) * (c ** 2 + lam)) ** 0.5)
def simpson(fn,og,bg,punten):
punten=max(punten,3)
if punten%2==0:
punten+=1
n=int((punten-1)//2)
h=(bg-og)/(punten-1)
s4=0
s2=0
x=og+h
for k in range(n):
s4+=fn(x)
x+=2*h
x=og+2*h
for k in range(n-1):
s2+=fn(x)
x+=2*h
s=h*(fn(og)+fn(bg)+4*s4+2*s2)/3
return s
def integreer(fn,u):
bg1=1e16
bg2=1000000.0*a**2
punten=100000
s=simpson(fn,u,bg1,punten)+simpson(fn,bg1,bg2,punten)+2/np.sqrt(bg2)
return s
################# evenaar
print("evenaar")
x = 0
y = re
z = 0
t1 = np.linspace(0, 1e15,100000)
a=a0
b=b0
c=c0
plt.figure()
plt.plot(t1, integrand(t1))
plt.show()
rho=rhok-rhom
a=rk
b=rk
c=rk
# u_ev = scipy.optimize.fsolve(fun, 10)
# u_ev = u_ev[0]
# # print(u_ev)
u=y**2-b**2
pek=G*np.pi*rho*a*b*c*integreer(integrand, u)
rho=rhom
a=a0
b=b0
c=c0
pem=G*np.pi*rho*a*b*c*integreer(integrand, 0)
pe=pek+pem
print(pek,pem,pe)
################# noordpool
print("noordpool")
x = 0
y = 0
z = rp
rho=rhok-rhom
a=rk
b=rk
c=rk
# u_np = scipy.optimize.fsolve(fun, 10)
# u_np = u_np[0]
u=z**2-c**2
pnk=G*np.pi*rho*a*b*c*integreer(integrand, u)
rho=rhom
a=a0
b=b0
c=c0
pnm=G*np.pi*rho*a*b*c*integreer(integrand, 0)
pp=pnk+pnm
print(pnk,pnm,pp)
print("klassieke formules")
print((G * m / rp))
print((G * m / re))
print('Verschil pool-evenaar',pp-pe)
print('rotatiepotentiaal',0.5*rp**2*w**2)
