pentachorons.jl

using Luxor, ColorSchemes

struct Point4D <: AbstractArray{Float64, 1}
    x::Float64
    y::Float64
    z::Float64
    w::Float64
end

Point4D(a::Array{Float64, 1}) = Point4D(a...)

Base.size(pt::Point4D) = (4, )

Base.getindex(pt::Point4D, i) = [pt.x, pt.y, pt.z, pt.w][i]

struct Point3D <: AbstractArray{Float64, 1}
    x::Float64
    y::Float64
    z::Float64
end

Base.size(pt::Point3D) = (3, )

function convert(Point, pt3::Point3D)
    k = 1/(K - pt3.z)
    return Point(pt3.x * k, pt3.y * k)
end

const K = 4.0

function convert(Point3D, pt4::Point4D)
    k = 1/(K - pt4.w)
    return Point3D(pt4.x * k, pt4.y * k, pt4.z * k)
end

function rotate4(A, matrixfunction)
    return map(A) do pt4
        Point4D(matrixfunction * pt4)
    end
end

function XY(θ)
    [cos(θ)   -sin(θ)  0        0;
    sin(θ)    cos(θ)   0        0;
    0         0        1        0;
    0         0        0        1]
    end

function XZ(θ)
    [cos(θ)   0        -sin(θ)  0;
    0         1        0        0;
    sin(θ)    0        cos(θ)   0;
    0         0        0        1]
    end

function XW(θ)
    [cos(θ)   0        0        -sin(θ);
    0         1        0        0;
    0         0        1        0;
    sin(θ)    0        0        cos(θ)]
    end

function YZ(θ)
    [1        0        0        0;
    0         cos(θ)   -sin(θ)  0;
    0         sin(θ)   cos(θ)   0;
    0         0        0        1]
    end

function YW(θ)
    [1        0        0        0;
    0         cos(θ)   0        -sin(θ);
    0         0        1        0;
    0         sin(θ)   0        cos(θ)]
    end

function ZW(θ)
    [1        0        0        0;
    0         1        0        0;
    0         0        cos(θ)   -sin(θ);
    0         0        sin(θ)   cos(θ)];
    end

function flatten(shape4)
    return map(pt3 -> convert(Point, pt3), map(pt4 -> convert(Point3D, pt4), shape4))
end

const n = -1/√5

const pentachoron = [Point4D(vertex...) for vertex in [
    [ 1.0,  1.0,  1.0, n],
    [ 1.0, -1.0, -1.0, n],
    [-1.0,  1.0, -1.0, n],
    [-1.0, -1.0,  1.0, n],
    [ 0.0,  0.0,  0.0, n + √5]]]

const pentachoronfaces = [
    [1, 2, 3],
    [1, 2, 4],
    [1, 2, 5],
    [1, 3, 4],
    [1, 3, 5],
    [1, 4, 5],
    [2, 3, 4],
    [2, 3, 5],
    [2, 4, 5],
    [3, 4, 5]]

XYW(a) = XY(a) * XW(a)
XZW(a) = XZ(a) * ZW(a)
YZW(a) = YZ(a) * ZW(a)
XZYW(a) = XZ(a) * YW(a)

XYW′(a) =  XW(a) * XY(a)
XZW′(a) =  ZW(a) * XZ(a)
YZW′(a) =  ZW(a) * YZ(a)
XZYW′(a) = YW(a) * XZ(a)

function snapshot(f, scalefactor)
    setlinejoin("bevel")
    setopacity(0.3)
    setline(1.0)
    pentachoron2D = flatten(
        rotate4(
            pentachoron,
                YZ(f * 2π) *
                #YW(f * 2π) *
                XZ(f * 2π) *  YW(f * 2π)
                #XY(f * 2π)
                ))
    for (n, face) in enumerate(pentachoronfaces)
        sethue(get(ColorSchemes.diverging_rainbow_bgymr_45_85_c67_n256,
            n/length(pentachoronfaces)))
        poly(scalefactor * pentachoron2D[face], :fillpreserve, close=true)
        sethue("black")
        strokepath()
    end
end

function drawtable(w, h, fname;
        rows = 20,
        cols = 20,
        sf = 300 # scale/fudge factor
        )
    Drawing(w, h, :png)
    origin()
    background("mistyrose4")
    t = Tiler(w, h, rows, cols, margin=30)
    tt = t.ncols * t.nrows
    for (pos, n) in t
        @layer begin
            translate(pos)
            snapshot(rescale(n, 1, tt, 0, 100), sf)
        end
    end
    finish()
    preview()
end

drawtable(600, 600, "/tmp/pentachorons.svg",
    rows=10,
    cols=10,
    sf=300)