Demo entry 6731122

Homework III

   

Submitted by anonymous on Apr 10, 2018 at 12:37
Language: Python 3. Code size: 3.0 kB.

from scipy import special as ss
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm

gamma_basalt = 2700*9.8 # NOTE: basalt density
gamma_crust = 3300*9.8 # NOTE: crust density
E = 70*pow(10,9) # NOTE: young's modules
miu = 0.25 # NOTE: poisson's ratio
cthk = 100 # NOTE: crust thickness in km
T = cthk*1000 # NOTE: crust Thickness in m
R = 6350*pow(10,3) # NOTE: Earth's radius

gamma_prime = gamma_crust+E*T/pow(R,2)
D = E*pow(T,3)/(12*(1-pow(miu,2)))
l = pow(D/gamma_prime,0.25)
print ("gamma_prime=%.2f kg m^3, D=%.3g N m , l=%.3f km" % (gamma_prime/9.8, D ,l/1000))
N1 = 50 # NOTE: Height steps, Height=N1*dh
dh = 2000/N1 # NOTE: Height dh
N2 = 100 # NOTE: time steps, time=N2*dt
dt = 2/N2 # NOTE: time dh
N3 = 120 # NOTE: distance steps, dist=N3*dk
dk = 10000 # NOTE: distance dk
coe = (gamma_basalt*dh)/gamma_prime

A = np.zeros((N1,N2),dtype=np.float64)
a = np.zeros((N1,N2),dtype=np.float64)

for i in range(N1):
    for j in range(N2):
        A[i][j] =  240000*(dt*(j+1))+(60000/N1)*(dt*(j+1))*i+(30000/N1)*(dt*j)
        a[i][j] = A[i][j]/l

wi = np.zeros((N1,N2,N3),dtype=np.float64)
wo = np.zeros((N1,N2,N3),dtype=np.float64)
wult = np.zeros((N3,N2),dtype=np.float64)
print ("calculating crustal deflection change:")
for j in tqdm(range(N2)):
    for k in range(N3):
        for i in range(N1):
            wi[i][j][k] = coe* ( a[i][j]*ss.kerp(a[i][j])*ss.ber(dk*(k+1)/l) - a[i][j]*ss.keip(a[i][j])*ss.bei(dk*(k+1)/l) + 1 )
            wo[i][j][k] = coe* ( a[i][j]*ss.berp(a[i][j])*ss.ker(dk*(k+1)/l) - a[i][j]*ss.beip(a[i][j])*ss.kei(dk*(k+1)/l) )

            if (dk*(k+1) <= A[i][j]):
                wult[k][j] = wult[k][j] - wi[i][j][k]
            else:
                wult[k][j] = wult[k][j] - wo[i][j][k]

print ("\ncalculating topography change:")
topo = np.zeros((N3,N2),dtype=np.float64)
for j in tqdm(range(N2)):
    for k in range(N3):
        if (dk*(k+1)<=dt*(j+1)*240000):
            topo[k][j] = 2000+wult[k][j]
        elif (dk*(k+1)>=dt*(j+1)*300000):
            topo[k][j] = wult[k][j]
        else:
            topo[k][j] = 1000*(6*(j+1)-10*(k+1))/(0.6*(j+1))+wult[k][j]

plt.figure()
j=[0,100]
plt.xlim(0,100)
plt.ylim(-2250,250)
plt.plot(wult[14], label="x=150km", lw=1.5)
plt.plot(wult[29], label="x=300km", lw=1.5)
plt.plot(wult[44], label="x=450km", lw=1.5)
plt.xticks(np.arange(0,100,20), np.arange(0,2,step=0.4))
plt.xlabel("time(Ma)")
plt.ylabel("deflection(m)")
plt.legend()
plt.title("deflection-time curves -- crust thickness=%dkm" % cthk)
plt.savefig('deflection-time--ctk=%dkm.jpg' % cthk)

plt.figure()
j=[0,100]
plt.xlim(0,100)
plt.ylim(-750,1500)
plt.plot(topo[14], label="x=150km", lw=1.5)
plt.plot(topo[29], label="x=300km", lw=1.5)
plt.plot(topo[44], label="x=450km", lw=1.5)
plt.xticks(np.arange(0,100,20), np.arange(0,2,step=0.4))
plt.xlabel("time(Ma)")
plt.ylabel("topography(m)")
plt.legend()
plt.title("topography-time curves -- crust thickness=%dkm" % cthk)
plt.savefig('topography-time--ctk=%dkm.jpg' % cthk)

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).