import numpy as np

def axes_to_sdr(p_trend, p_plunge, t_trend, t_plunge):
    # Convert degrees to radians
    p_tr = np.radians(p_trend)
    p_pl = np.radians(p_plunge)
    t_tr = np.radians(t_trend)
    t_pl = np.radians(t_plunge)
    
    # Calculate P and T unit vectors in NED coordinates
    # x: North, y: East, z: Down
    P = np.array([
        np.cos(p_pl) * np.cos(p_tr),
        np.cos(p_pl) * np.sin(p_tr),
        np.sin(p_pl)
    ])
    
    T = np.array([
        np.cos(t_pl) * np.cos(t_tr),
        np.cos(t_pl) * np.sin(t_tr),
        np.sin(t_pl)
    ])
    
    # Normalize just to be sure
    P /= np.linalg.norm(P)
    T /= np.linalg.norm(T)
    
    # The normal vector n and slip vector u are 45 degrees between P and T
    # n = (T + P) / sqrt(2), u = (T - P) / sqrt(2) or vice versa
    # Let's compute both combinations to get the two nodal planes
    # Plane A
    nA = (T + P) / np.sqrt(2)
    uA = (T - P) / np.sqrt(2)
    
    # Plane B (swapping normal and slip)
    nB = uA.copy()
    uB = nA.copy()
    
    def get_sdr(n, u):
        # Make sure normal vector points upwards to respect Aki & Richards / standard convention where z is up?
        # Standard Aki & Richards convention: x=North, y=East, z=Up.
        # Let's convert n and u to North, East, Up (NEU)
        # In NED: n_neu = [n[0], n[1], -n[2]]
        n_neu = np.array([n[0], n[1], -n[2]])
        u_neu = np.array([u[0], u[1], -u[2]])
        
        # If normal points down (z < 0), flip both vectors to keep it pointing into the upper hemisphere
        if n_neu[2] < 0:
            n_neu = -n_neu
            u_neu = -u_neu
            
        # Dip
        dip = np.arccos(n_neu[2]) # since n is unit vector, n_z = cos(dip)
        
        # Strike
        if abs(n_neu[2]) > 0.99999: # Horizontal plane
            strike = 0.0
            dip = 0.0
            rake = np.arctan2(u_neu[1], u_neu[0])
        else:
            # Strike vector is horizontal and perpendicular to normal: s = (n_y, -n_x, 0)
            strike = np.arctan2(-n_neu[0], n_neu[1])
            if strike < 0:
                strike += 2 * np.pi
                
            # Rake
            # s vector components
            s_vec = np.array([np.sin(strike), np.cos(strike), 0.0]) # wait, standard right hand rule strike vector:
            # if strike is phi, strike vector is (sin phi, cos phi, 0) in NEU? 
            # Let's double check standard formula:
            # n = (-sin(strike)*sin(dip), cos(strike)*sin(dip), cos(dip))
            # So sin(strike)*sin(dip) = -n_x, cos(strike)*sin(dip) = n_y
            # Therefore strike = arctan2(-n_x, n_y) -> this matches.
            
            # Let's re-verify strike vector direction:
            # s = (cos(strike), sin(strike), 0) or (sin(strike), cos(strike), 0)
            # In standard seismology (Aki & Richards):
            # s = (cos(phi), sin(phi), 0)
            # Let's find phi: tan(phi) = n_y / -n_x => phi = arctan2(-n_x, n_y)? No, tan(phi) = sin(phi)/cos(phi).
            # If n_x = -sin(phi)*sin(delta) and n_y = cos(phi)*sin(delta)
            # then sin(phi) = -n_x/sin(delta), cos(phi) = n_y/sin(delta)
            # phi = arctan2(-n_x, n_y)
            phi = np.arctan2(-n_neu[0], n_neu[1])
            if phi < 0: phi += 2*np.pi
            
            s_vec = np.array([np.cos(phi), np.sin(phi), 0.0])
            # Slip vector u = cos(rake)*s + sin(rake)*upper_dip_vector
            # Or rake = arctan2(u . dip_dir_vector, u . s_vec)
            # dip_dir_vector = s_vec x n_neu?
            # Actually, standard formula for u in Aki & Richards:
            # u_x = cos(rake)*cos(phi) + sin(rake)*cos(delta)*sin(phi)
            # u_y = cos(rake)*sin(phi) - sin(rake)*cos(delta)*cos(phi)
            # u_z = sin(rake)*sin(delta)
            # So rake = arctan2(u_z / sin(delta), u_x*cos(phi) + u_y*sin(phi))
            u_z = u_neu[2]
            u_x_prime = u_neu[0]*np.cos(phi) + u_neu[1]*np.sin(phi)
            rake = np.arctan2(u_z / np.sin(dip), u_x_prime)
            
        return np.degrees(phi), np.degrees(dip), np.degrees(rake)

    s1, d1, r1 = get_sdr(nA, uA)
    s2, d2, r2 = get_sdr(nB, uB)
    return (s1, d1, r1), (s2, d2, r2)

res = axes_to_sdr(172.1, 1.4, 262.2, 6.3)

strike1= res[0][0]
dip1= res[0][1]
rake1= res[0][2]

print(res[0])
print(res[1])
print(f"{res[0][0]:8.1f} {res[0][1]:8.1f} {res[0][2]:8.1f}")
print(f"{res[1][0]:8.1f} {res[1][1]:8.1f} {res[1][2]:8.1f}")

print(f"{strike1:8.1f} {dip1:8.1f} {rake1:8.1f}")