import numpy as np
def compute_Sg(S, angles=(0,0,0)):
alpha, beta, gamma = np.radians(angles)
Rg = np.array([[np.cos(alpha) * np.cos(beta),
np.sin(alpha) * np.cos(beta),
-np.sin(beta)],
[np.cos(alpha) * np.sin(beta) * np.sin(gamma) - np.sin(alpha) * np.cos(gamma),
np.sin(alpha) * np.sin(beta) * np.sin(gamma) + np.cos(alpha) * np.cos(gamma),
np.cos(beta) * np.sin(gamma)],
[np.cos(alpha) * np.sin(beta) * np.cos(gamma) + np.sin(alpha) * np.sin(gamma),
np.sin(alpha) * np.sin(beta) * np.cos(gamma) - np.cos(alpha) * np.sin(gamma),
np.cos(beta) * np.cos(gamma)]])
return np.dot(Rg.T, np.dot(S,Rg))
S = np.diag([60, 50, 45])
S_G = compute_Sg(S, angles=(180,0,90)); S_G
array([[ 6.00000000e+01, -1.83697020e-15, -3.74939946e-32], [ -1.83697020e-15, 4.50000000e+01, -3.06161700e-16], [ -3.74939946e-32, -3.06161700e-16, 5.00000000e+01]])
def compute_unit_vectors(strike, dip):
strike = np.radians(strike)
dip = np.radians(dip)
n = np.array([-np.sin(strike) * np.sin(dip), np.cos(strike) * np.sin(dip), -np.cos(dip) ])
ns = np.array([ np.cos(strike), np.sin(strike), 0 ])
nd = np.array([ -np.sin(strike) * np.cos(dip), np.cos(strike) * np.cos(dip), np.sin(dip) ])
return (n, ns, nd)
i) $strike = 45^{\circ}, dip = 65^{\circ}$
n, ns, nd = compute_unit_vectors(45, 65);
print(n)
print(ns)
print(nd)
[-0.64085638 0.64085638 -0.42261826] [ 0.70710678 0.70710678 0. ] [-0.29883624 0.29883624 0.90630779]
Now we compute the normal and shear stresses on the plane
S_G = np.array([[47.5, -12.5, 0],[-12.5, 47.5, 0],[0,0,40]])
S_n = np.dot(np.dot(S_G, n), n); S_n
56.427876096865383
tau_s = np.dot(np.dot(S_G, n), ns); tau_s
-7.1054273576010019e-15
tau_d = np.dot(np.dot(S_G, n), nd); tau_d
7.6604444311897772
tau_mag = np.sqrt(tau_d ** 2 + tau_s ** 2); tau_mag
7.6604444311897772
Pc = (0.6 * S_n - tau_mag) / 0.6; Pc
43.660468711549093