Trägheitsmoment und Umrechnung in kg eingefügt

AdamsWilliamsonCMB.py

@@ -54,6 +54,9 @@ def mass_correction(rho, r):
         res += (4/3) * pi * rho[i] * -(r[i]**3-r[i-1]**3)
     return res
 
+def moment_of_inertia(r1,r2,M):
+    return 2/5 * M *(r1**5 - r2**5)/(r1**3 - r2**3)
+
 
 
 file = os.getcwd() + r"\VpVs-valuesak135.csv"   
@@ -62,7 +65,7 @@ file = os.getcwd() + r"\VpVs-valuesak135.csv"
 depth, density, vp, vs = read_data(file)  # Tiefe, Dichte, vp, vs
 
 # r geht vom Mittelpunkt nach außen
-r = 6371 - depth
+r = (6371 - depth) * 1000
 
 # Startbedingungen
 rho_0 = 3000  # kg/m^3
@@ -80,14 +83,23 @@ for i in range(9, len(rho_ex)):
         rho_ex[i] = 9900
         rho_im[i] = 9900
     else:
-        rho_ex[i] = density_expansion_explicit(rho_ex[i-1], vp[i-1], vs[i-1], r[i-1], r[i])
-        rho_im[i] = density_expansion_implicit(rho_im[i-1], vp[i-1], vs[i-1], r[i-1], r[i])
+        rho_ex[i] = density_expansion_explicit(rho_ex[i-1], vp[i-1]*1000, vs[i-1]*1000, r[i-1], r[i])
+        rho_im[i] = density_expansion_implicit(rho_im[i-1], vp[i-1]*1000, vs[i-1]*1000, r[i-1], r[i])
         masse[i] = mass(r[i-1], rho_ex[i-1])
 
+# Berechnung des Trägkeitsmomnts
+J=np.zeros(len(masse))
+for i in range(len(masse)-1):
+    if r[i] !=r[i+1]:
+        J[i]=moment_of_inertia(r[i],r[i+1],masse[i])
+        
+J[-1]=moment_of_inertia(r[i],0,masse[i])
+J_ges=sum(J)
 
 print("Earth mass explicit: ", mass_correction(rho_ex, r))
 print("Earth mass implicit: ", mass_correction(rho_im, r))
 print("Earth mass exact: ", mass_correction(density, r))
+print(J_ges)
 
 
 plt.plot(depth,rho_ex)