#!/usr/bin/env python
# coding: utf-8

# # Problem 3 (30 points)

# In[1]:


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))


# In[2]:


S = np.diag([60, 50, 45])

S_G = compute_Sg(S, angles=(180,0,90)); S_G


# In[3]:


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}$**

# In[4]:


n, ns, nd = compute_unit_vectors(45, 65); 
print(n)
print(ns)
print(nd)


# Now we compute the normal and shear stresses on the plane

# In[5]:


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


# In[6]:


tau_s = np.dot(np.dot(S_G, n), ns); tau_s


# In[7]:


tau_d = np.dot(np.dot(S_G, n), nd); tau_d


# In[8]:


tau_mag = np.sqrt(tau_d ** 2 + tau_s ** 2); tau_mag


# In[9]:


Pc = (0.6 * S_n - tau_mag) / 0.6; Pc


# In[ ]: