#http://kobayashi.hub.hit-u.ac.jp/topics/simplex.pdf
#Cayley-Menger matrix
#H.S.M. Coxeter, The Circumradius of the General Simplex, The Mathematical Gazette 15, 210 (1930) 229-231

import numpy as np
from math import factorial
from numpy import sqrt


def determinant(m):
    return np.linalg.det(m)

def distance_matrix(points):
    n = len(points)
    distances = np.zeros((n, n))
    for i in range(n):
        for j in range(i + 1, n):
            d = np.array(points[i]) - np.array(points[j])
            distances[i, j] = np.dot(d, d)
            distances[j, i] = distances[i, j]

    return distances

def cayley_menger_matrix(d_matrix):
    n = len(d_matrix) - 1
    cm_matrix = np.ones((n + 2, n + 2))
    cm_matrix[1:, 1:] = np.array(d_matrix)
    cm_matrix[0, 0] = 0
    
    return cm_matrix

def simplex_volume(n, det_cm):
    v2 = (-1)**(n + 1) / 2**n / factorial(n)**2 * det_cm

    return sqrt(v2)

def circumradius(det_d, det_cm):
	return sqrt(-det_d / 2.0 / det_cm)

def basetriangle_area(p1, p2, p3):
    d12 = np.array(p2) - np.array(p1)
    d13 = np.array(p3) - np.array(p1)
    return 0.5 * abs(d12[0] * d13[1] - d12[1] * d13[0])

def triangle_area(p1, p2, p3):
    ab = np.array(p2) - np.array(p1)
    ac = np.array(p3) - np.array(p1)
    cross_product = np.cross(ab, ac)
    S = 0.5 * np.linalg.norm(cross_product)

    return S

def plane_equation(p1, p2, p3):
    v1 = np.array(p2) - np.array(p1)
    v2 = np.array(p3) - np.array(p1)
    normal_vector = np.cross(v1, v2)
    a, b, c = normal_vector
    d = -np.dot(normal_vector, p1)
    return a, b, c, d

def point_plane_distance(a, b, c, d, p4):
    x4, y4, z4 = p4
    distance = abs(a * x4 + b * y4 + c * z4 + d) / np.sqrt(a**2 + b**2 + c**2)
    return distance

def tirgonalpyramid_height(p1, p2, p3, p4):
    a, b, c, d = plane_equation(p1, p2, p3)
    h = point_plane_distance(a, b, c, d, p4)

    return h

def trigonalpyramid_volume(p1, p2, p3, p4):
    S = triangle_area(p1, p2, p3)
    h = tirgonalpyramid_height(p1, p2, p3, p4)
    return S * h / 3.0


def main():
# 例: 3次元空間で4点が作る単体の体積と外接円の半径を計算
    vertices = [        # 正四面体
        [1, 1, 1],
        [-1, -1, 1],
        [-1, 1, -1],
        [1, -1, -1]
    ]
    '''
    vertices = [
        [0, 0, 0],
        [1, 1, 0],
        [-1, 1, 0],
        [0, 0, 1]
    ]
    '''

    ndim = len(vertices) - 1
    print()
    print("ndim=", ndim)
    print("vertices:\n", vertices)

    Sbase = triangle_area(*vertices[:3])
    htrip = tirgonalpyramid_height(*vertices)
    Vtrip = trigonalpyramid_volume(*vertices)
    print()
    print("底辺三角形の面積:", Sbase)
    print("4面体の高さ:", htrip)
    print("4面体の体積:", Vtrip)

    print()
    d_matrix  = distance_matrix(vertices)
    cm_matrix = cayley_menger_matrix(d_matrix)
    print("d_matrix:\n", d_matrix)
    print("cm_matrix:\n", cm_matrix)

    det_d  = determinant(d_matrix)
    det_cm = determinant(cm_matrix)
    print("det_d=", det_d)
    print("det_cm=", det_cm)

    volume = simplex_volume(ndim, det_cm)
    radius = circumradius(det_d, det_cm)

    print(f"単体の体積: {volume}")
    print(f"外接円の半径: {radius}")
    print()
    

if __name__ == "__main__":
    main()
        
