Triangles_mDunitball.jl

# arvotaan m-ulottesen yksikköpallon sisään kolmioita
# tutkitaan kuinka suuri osa kolmioista on tylppäkulmaisia
# m-ulotteisen pallon sisällä
# Juhani Kaukoranta 18.1.2023
using LinearAlgebra,AngleBetweenVectors,Random,Distributions

function mDballquadro(m,n)
    # m = ulottuvuus,n kolmioiden määrä
    #d = Normal(0,1)
    jananpituus = 0
    area1 = 0 # nelikulmion puoliskon ala = kolmion ala
    area2 = 0 # nelikulmien puolisko
    piiri = 0 # nelikulmion piiri
    kovera = 0 # laskee koverien (concave) nelikulmioiden määrän,yksi kulma > 180
    kolmionpiiri = 0
    tylppa = 0 # joku kolmion kulmista > 90 astetta
    acute = 0 # kaikki kulmat < 90
    suorakulmainen = 0

    for i = 1 : n
        d = Normal(0,1)
        r = rand(4) .^(1/m) # pallon säde kärkipisteille uniformajakaumasta m-juurikorjattuna
        u = rand(d,4,m) # m kpl muuttujia, kaikilla 4 arvoa
        norm1 = norm(u[1,:]) # muuttujien normit
        norm2 = norm(u[2,:])
        norm3 = norm(u[3,:])
        norm4 = norm(u[4,:])
        normi = [norm1;norm2;norm3;norm4]
        x = @. r * u / normi # 4 normitettua m-ulotteista muuttujaa
        P1 = x[1,:] # nelikulmion kärkipiste (m kpl koordinaatteja)
        P2 = x[2,:] # nelikulmion kärkipiste
        P3 = x[3,:] # nelikulmion kärkipiste
        P4 = x[4,:] # nelikulmion kärkipiste
        u1 = P1 .- P2 # kärjestä P1
        v1 = P1 .- P4 # kärjestä P1
        area1 += 0.5*sqrt(dot(u1,u1)*dot(v1,v1) - dot(u1,v1)^2) #P1 kolmion ala
        u2 = P2 .- P3 # kärjestä P3
        v2 = P4 .- P3 # kärjestä P3
        area2 += 0.5*sqrt(dot(u2,u2)*dot(v2,v2) - dot(u2,v2)^2) #P3 kolmion ala
        angleP1 = angle(P2 .- P1,P4 .- P1) # nelikulmion kärkikulmia
        angleP2 = angle(P3 .- P2,P1 .- P2)
        angleP3 = angle(P4 .- P3,P2 .- P3)
        angleP4 = angle(P1 .- P4,P3 .- P4)
        alfa = angleP1 # kolmion kulma
        beta = angle(P4-P2,P1-P2) # kolmion kulma
        gamma = angle(P2-P4,P1-P4) # kolmion kulma
        # jos yksikin kolmion kulmista > 90°,kolmio on tylppäkulmainen
        if alfa > pi/2 || beta > pi/2 || gamma > pi/2
            tylppa += 1
        end
        if abs(alfa-pi/2) < 1.0e-6 && abs(beta-pi/2)<1.0e-6  && abs(gamma-pi/2)<1.0e-6
            acute += 1
        end
        if abs(alfa-pi/2)<1.0e-6 || abs(beta-pi/2)<1.0e-6 || abs(gamma-pi/2)<1.0e-6
            suorakulmainen += 1
        end
        if abs(angleP1 + angleP2 +angleP3 + angleP4 - 2*pi) <1.0e-6
            kovera += 1 # nelikulmio on kovera (concave)
        end
        a, b, c, d = norm(u1), norm(u2), norm(v2),norm(v1)
        lavistaja = norm(P2-P4)
        piiri += a+b+c+d # nelikulmion piiri
        kolmionpiiri += a+b+lavistaja
        jananpituus += norm(u1) #
    end
    println("koveria (concave) nelikulmioita = ",kovera)
    println("koverien osuus =", kovera/n)
    println("nelikulmioiden pinta-alan keskiarvo = ",(area1+area2)/n)
    println("nelikulmioiden piirien keskiarvo = ",piiri/n)
    println("kolmioiden piirien keskiarvo = ",kolmionpiiri/n)
    println("kolmioiden pinta-alojen keskiarvo = ",area1/n)
    println("tylppäkulmaisten kolmioiden osuus = ",tylppa/n)
    println("teräväkulmaisten kolmioiden osuus = ",1-tylppa/n)
    println("suorakulmaisia kolmioita(radiaanitarkkuus 1.e-6) = ",suorakulmainen)
    println("janojen keskipituus = ",jananpituus/n)
end