Tetrahedrons_in_unitball.jl

# tetrahedrons in dim-dimensional unitball
# Juhani Kaukoranta 21.1.2023
using LinearAlgebra,Random,Distributions
function mDballtetrahedrons(dim,n)
    # dim = ulottuvuus, n = tetraedrien lukumäärä
    # lasketaan yksikköpallon tetraedrien keskitilavuus ja keskiala

    Volume = 0 # kerää tetrahedrien tilavuuksien summaa
    kokonaisala = 0 # kerää kokonaispinta-alan summaa
    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 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

    for i = 1 : n
        d = Normal(0,1)
        r = rand(4) .^(1/dim) # uniformjakauman m-juuresta säteet
        u = rand(d,4,dim) # dim kpl muuttujia, kaikilla 4 arvoa
        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)
        tahkot = [tahko(d12,d13,d23) tahko(d12,d14,d24) tahko(d23,d24,d34) tahko(d13,d14,d34)]
        kokonaisala += sum(tahkot)
        #korkeudet =  3*Volume/tahkot


    end

    println(n," tetraedria satunnaisesti arvottuna ",dim,"-ulotteisen yksikköpallon sisälle")
    println("Tedraedrien tilavuuksien keskiarvo = ", Volume/n)
    println("tetraedrien kokonaispinta-alan keskiarvo =",kokonaisala/n)

end