Tetraedriunitball.jl

# arvotaan 3-ulotteisen yksikköpallon sisään tetraedeja
# lasketaan tetraedrien keskitilavuus
# Juhani Kaukoranta 20.1.2023
using LinearAlgebra,Random,Distributions
function mDballtetrahedrons(n)
    # n = tetraedrien lukumäärä
    # lasketaan tetraedrien keskitilavuus yksikköpallossa
    m = 3 # ulottuvuus, dimensio
    tilavuus = 0 # tetraedrien tilavuuksien summa
    V = 0 # determinantin avulla laskettu
    kokonaisala = 0 # laskee pinta-alan
    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

    for i = 1 : n
        d = Normal(0,1)
        r = rand(4) .^(1/m) # uniformjakauman m-juuresta säteet
        u = rand(d,4,m) # m 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, huippu
        # P1P2P3 muodostavat "pohjakolmion", P4 huippu
        a = norm(P2 .- P1) # pohjakolmion sivu
        b = norm(P3 .- P2) # pohjakolmion sivu
        c = norm(P3 .- P1) # pohjakolmion sivu
        d = norm(P4 .- P1)
        e = norm(P4 .- P2)
        f = norm(P4 .- P3)
        pohjanala = tahko(a,b,c) # Heron, pohjakolmion pinta-ala
        sivutahkot = tahko(a,e,d)+tahko(b,f,e)+tahko(c,d,f)
        kokonaisala += pohjanala + sivutahkot
        #N = cross((P2-P1),(P3-P1)) # pohjakolmion normaali
        #korkeus = abs(dot((P4-P1),N))/norm(N) #tetraedrin korkeus
        #tilavuus += pohjanala * korkeus/3 # tetraedrin tilavuuksien summa
        M = [P1-P4 P2-P4 P3-P4] # determinanttia varten
        determ = abs(det(M))/6
        V += determ # tilavuus
    end

    #println("tetraedrien tilavuuksien keskiarvo =",tilavuus/n)
    println("tilavuuksien keskiarvo = ", V/n)
    println("tetraedrien kokonaispinta-alan keskiarvo =",kokonaisala/n)

end