#! /usr/bin/env python3

import numpy as np
import numpy.linalg as linalg

K = np.array( [ [ 2, 1/np.sqrt(3), 4 ],
                [ 0, 2,            4 ],
                [ 0, 0,            1 ] ] )
S = linalg.inv(K.T) @ linalg.inv( K )

def innerproduct(a,b): return a.T @ S @ b

def angle(a,b): 
    num = innerproduct(a, b) 
    den = np.sqrt( innerproduct(a,a) * innerproduct(b,b) )
    print ( num.item() )
    print ( den.item() )
    return np.arccos( num.item() / den.item() )

A = np.array( [ 0,0,1 ] ).reshape(3,1)
B = np.array( [ 2,0,1 ] ).reshape(3,1)
C = np.array( [ 2,1,1 ] ).reshape(3,1)

A0 = K @ A
B0 = K @ B
C0 = K @ C

print ( "A' = ", A0.T )
print ( "B' = ", B0.T )
print ( "C' = ", C0.T )

AB = B0 - A0
AC = C0 - A0

print( "AB = ", AB.T )
print( "AC = ", AC.T )

alpha0 = angle( AB, AC )
beta0 = angle( B0-C0, B0-A0 )

print ( "alpha0 = ", alpha0, " = ", alpha0*180/np.pi, " degrees" )
print ( "beta0 = ", beta0, " = ", beta0*180/np.pi, " degrees" )