# tetrahedrons in dim-dimensional space
using LinearAlgebra,Random,Distributions
# Juhani Kaukoranta 18.12.2023
function TetrahedronVolume(d12,d13,d14,d23,d24,d34)
    m = [0 1 1 1 1; 1 0 d12^2 d13^2 d14^2;1 d12^2 0 d23^2 d24^2; 1 d13^2 d23^2 0 d34^2; 1 d14^2 d24^2 d34^2 0]
    return sqrt(1/288*det(m))
end
function tahko(a,b,c)
    # laskee tahkon pinta-alan
    p = (a+b+c)/2 # tahkon sivut
    return sqrt(p*(p-a)*(p-b)*(p-c))
end
function TetrahedronArea(d12,d13,d14,d23,d24,d34)
    return sum([tahko(d12,d13,d23) tahko(d12,d14,d24) tahko(d23,d24,d34) tahko(d13,d14,d34)])
end
function mDballtetrahedrons(dim,n)
    # dim = ulottuvuus, n = tetraedrien lukumäärä  
    Volume = 0 # tetraedrien tilavuuksien summa
    Ala = 0 # tetraedrien pinta-alojen summa
    for i = 1 : n
        r = rand(4) .^(1/dim) # uniformjakauman dim-juuresta säteet
        u = randn(4,dim) # 4 kpl m-ulotteista normaalijakatunutta kärkipistettä 
        norm1 = norm(u[1,:]) # kärki P1 normi
        norm2 = norm(u[2,:]) # kärki P2 normi
        norm3 = norm(u[3,:]) # kärki P3 normi
        norm4 = norm(u[4,:]) # kärki P4 normi
        normi = [norm1;norm2;norm3;norm4]
        x = @. r * u / normi # 4 normitettua m-ulotteista muuttujaa
        P1 = x[1,:] # tetraedrin kärkipiste (m kpl koordinaatteja)
        P2 = x[2,:] # tetraedrin kärkipiste
        P3 = x[3,:] # tetraedrin kärkipiste 
        P4 = x[4,:] # teraedrin kärkipiste
        d12 = norm(P1-P2) # tetraedrit sivujen pituudet
        d13 = norm(P1-P3)
        d14 = norm(P1-P4)
        d23 = norm(P2-P3)
        d24 = norm(P2-P4)
        d34 = norm(P3-P4)
        Volume +=  TetrahedronVolume(d12,d13,d14,d23,d24,d34)
        Ala += TetrahedronArea(d12,d13,d14,d23,d24,d34)
    end   
    println(n," tetraedria satunnaisesti arvottuna ",dim,"-ulotteisen yksikköpallon sisälle")
    println("Tedraedrien tilavuuksien keskiarvo = ", Volume/n)
    println("tetraedrien kokonaispinta-alan keskiarvo =",Ala/n)
end