Merge branch 'main' of https://gitlab.rrz.uni-hamburg.de/bax8453/geodynamik-uebungen

AdamsWilliamsonCMB.py

@@ -26,7 +26,7 @@ def read_data(file) -> tuple:
 def mass(r: list, rho_r: list) -> float:
     '''Calculates difference of actual mass of earth and mass calculated by a given density model.
     If density model is not complete, the function returns the mass not yet incorporated in the model.'''
-    M_e=5.974*10**24
+    M_e=5.951e24
     result = M_e
     for i in range(len(r)-1):
         result -= (4/3)*np.pi*rho_r[i]*(r[i]**3-r[i+1]**3)
@@ -82,7 +82,7 @@ def calc_density(rho_0: int, rho_CMB: int, radius: list, vp: list, vs: list) ->
 
         # more discontinuities
         # elif depth[i] == 410 and depth[i-1] == 410:   # 410km discontinuity
-        #    rho_ex[i] = 3921.4285714285716
+        #    rho_ex[i] = 3931.7
         # elif depth[i] == 660 and depth[i-1] == 660:   # 660km discontinuity
         #    rho_ex[i] = 4238.7
         # elif depth[i] == 2740 and depth[i-1] == 2740:   # D" discontinuity
@@ -92,11 +92,11 @@ def calc_density(rho_0: int, rho_CMB: int, radius: list, vp: list, vs: list) ->
             # rho_im = density_expansion_implicit(rho_im, vp, vs, r, i)
     return rho_ex
 
-def best_inertia_model(file, n=31, rho_0_low=3500, rho_0_high=3700, rho_CMB_low=9600, rho_CMB_high=9900):  
+def best_inertia_model(file, n=31, rho_0_low=3600, rho_0_high=3700, rho_CMB_low=9450, rho_CMB_high=9950):  
     depth, density, vp, vs = read_data(file)
 
-    M_0 = 5.972e24  # Masse in 18km Tiefe
-    I_0 = 8.025e37  # Moment of Inertia
+    M_0 = 5.951e24  # Masse in 18km Tiefe
+    I_0 = 7.964e37  # Moment of Inertia, without the uppermost 18km
     vp *= 1000
     vs *= 1000
 
@@ -130,16 +130,16 @@ def best_inertia_model(file, n=31, rho_0_low=3500, rho_0_high=3700, rho_CMB_low=
     plt.figure()
     plt.plot(density_choice,radius_model, label='calculated density profile')
     plt.plot(density,radius_model, label= 'PREM: density profile')
-    plt.xlabel('density [kg/m³]')
-    plt.ylabel('radius [km]')
-    plt.legend()
+    plt.xlabel('density [kg/m³]', fontsize=20)
+    plt.ylabel('radius [m]', fontsize=20)
+    plt.legend(fontsize=20)
     plt.show()
 
-def best_mass_model(file, n=31, rho_0_low=3700, rho_0_high=3900, rho_CMB_low=9800, rho_CMB_high=10050):
+def best_mass_model(file, n=31, rho_0_low=3700, rho_0_high=3900, rho_CMB_low=9950, rho_CMB_high=10050):
     depth, density, vp, vs = read_data(file)
 
-    M_0 = 5.972e24  # Masse in 18km Tiefe
-    I_0 = 8.025e37  # Moment of Inertia
+    M_0 = 5.951e24  # Masse in 18km Tiefe
+    I_0 = 7.964e37  # Moment of Inertia, without the uppermost 18km
     vp *= 1000
     vs *= 1000
 
@@ -173,9 +173,9 @@ def best_mass_model(file, n=31, rho_0_low=3700, rho_0_high=3900, rho_CMB_low=980
     plt.figure()
     plt.plot(density_choice,radius_model, label='Calculated density profile')
     plt.plot(density,radius_model, label='PREM: density profile')
-    plt.xlabel('density [kg/m³]')
-    plt.ylabel('radius [km]')
-    plt.legend()
+    plt.xlabel('density [kg/m³]', fontsize=20)
+    plt.ylabel('radius [m]', fontsize=20)
+    plt.legend(fontsize=20)
     plt.show()
 
 
@@ -185,7 +185,7 @@ file = os.getcwd() + r"\VpVs-valuesak135.csv"
 
 print('Calculating optimal density model in regard to total mass. Depending on the grid size of the search domain this may take a while. \n')
 
-best_mass_model(file,n=20,rho_0_low=3600,rho_0_high=3800)
+best_mass_model(file,n=20,rho_0_low=3660,rho_0_high=3730)
 
 print('\nCalculating optimal density model in regard to inertia. Depending on the grid size of the search domain this may take a while. \n')