CMB

AdamsWilliamsonCMB.py

+import numpy as np
+import matplotlib.pyplot as plt
+import csv
+import scipy.integrate as integrate
+import os
+
+
+# Einlesen der Daten
+def read_data(file):
+    with open(file, 'r') as csvfile:
+        csvreader = csv.reader(csvfile, delimiter=";")
+        next(csvreader)
+        depth = []
+        rho = []
+        vp = []
+        vs = []
+        for row in csvreader:
+            if row:
+                depth.append(float(row[0]))  # km
+                rho.append(float(row[1]) * 1000)  # kg/m^3
+                vp.append(float(row[2]))  # km/s
+                vs.append(float(row[3]))  # km/s
+    return np.column_stack([depth]), np.column_stack([rho]), np.column_stack([vp]), np.column_stack([vs])
+
+
+def mass(r, rho_r):
+    '''Calculates mass of earth inside sphere of radius r.'''
+    mass = (4/3)*np.pi*rho_r*r**3
+    return mass
+
+
+# Berechnung der Dichte
+def density_expansion(rho, vp, vs, r1, r2):
+    '''Returns Euler approximation for rho at r2.
+    Method: 
+    rho(r+h) = rho(r) + h*d(rho)/dr(r)'''
+    G = 6.67 * 10**(-11)  # Gravitationskonstante
+    phi = vp**2 - (4/3) * vs**2
+    rho_r = rho + (r2-r1) * (-G * mass(r1, rho) * rho / (r1**2*phi))
+    return rho_r
+
+
+file = os.getcwd() + r"\VpVs-valuesak135.csv"   
+
+# Einlesen der Daten
+depth, density, vp, vs = read_data(file)  # Tiefe, Dichte, vp, vs
+
+# r geht vom Mittelpunkt nach außen
+r = 6371 - depth
+
+# Startbedingungen
+rho_0 = 3000  # kg/m^3
+M_0 = 5.973 * 10**24  # Masse in 18km Tiefe
+
+rho = np.zeros(len(density))
+masse = np.zeros(len(density))
+masse[8] = M_0
+rho[8] = rho_0
+
+for i in range(9, len(rho)):
+    if depth[i] == 2891.5:
+        rho[i] = 9400
+    else:
+        rho[i] = density_expansion(rho[i-1], vp[i-1], vs[i-1],r[i-1], r[i])
+        masse[i] = mass(r[i-1], rho[i-1])
+
+print(rho)
+print(masse)
+
+# Berechnung von rho und M für jede Schicht
+# for i in range(8, len(rho) - 1):  # Die ersten Werte werden weg gelassen
+#     if r[i, 0] != r[i-1, 0]:
+#         M[i+1] = mass(r[i, 0], r[i-1, 0], rho[i])
+#         rho[i+1] = density(M[i+1], rho[i], vp[i, 0], vs[i, 0], r[i, 0], r[i-1, 0])
+#     else:
+#         M[i+1] = M[i]
+#         rho[i+1] = rho[i]
\ No newline at end of file