-
Notifications
You must be signed in to change notification settings - Fork 4
/
stress_disp_tor.py
35 lines (27 loc) · 990 Bytes
/
stress_disp_tor.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
import numpy as np
from earthfun import earthfun
def stress_disp_tor(r, WT_0, l, omega, imod, rho, mu, rspan):
"""
Finds the derivatives of W(r) and T(r)
Args:
r (int or float): Radius at which to evaluate the derivatives
WT_0 (numpy.ndarray): WT_0[0] is W(r) and WT_0[1] is T(r)
l (int):
omega (int or float): Angular frequency
imod (int):
rho:
mu:
rspan:
Returns:
numpy.ndarray: dWT. dWT[0] is the derivative of W(r), and dWT[1] is the derivative of T(r)
"""
dWT = np.empty(2)
# structural values at radius r: density and rigidity
# note: if imod == 0, then we use the rho and mu from spshell.ipynb
if imod != 0:
rho, mu = earthfun(r, rspan, imod)
# displacement (first row of equation 1)
dWT[0] = WT_0[0] / r + WT_0[1] / mu
# stress (second row of equation 2)
dWT[1] = ((l-1)*(l+2)*mu/(r*r) - rho*omega*omega)*WT_0[0] - 3*WT_0[1]/r
return dWT