# OMEGA marked reference geometry — unit-scale binary64 materialization
# Paper VII v1.6, DOI 10.5281/zenodo.22233573.
#
# The symbolic formulas in the preprint define the exact object.
# This program is a Float64 evaluation, numerical check, and Plots reconstruction
# of that object in the canonical frame. It is not a theorem checker.
# LinearAlgebra.rank and isapprox verify the binary64 image; they do not prove
# Theorem 3.4, Proposition 8.1, Proposition 9.1, or Definition 11.1.
#
# Canonical frame: α = 0, u = (0,−1,0), v = (1,0,0), n = u × v.
# Foot parameter is ℓ (not λ: in this archive λ = φ/Ω₀ already).
# Requires: LinearAlgebra, Plots, Printf. plotlyjs() is interactive; gr() is a
# static fallback.
# CC BY 4.0

using LinearAlgebra
using Plots
using Printf

# 1. Scale and binary64 constants (bindings are constant; values are Float64)
const a = 1.0
const φ = (1.0 + sqrt(5.0)) / 2.0
const φ² = φ + 1.0
const ϵ = π / 648000.0  # 1 arcsecond in radians
const δ = ϵ / 2.0
const s = sin(δ)
const c = cos(δ)

const x₀ = s
const x₁ = 1.0 + s
const x₂ = φ² + s

# 2. Canonical placement in oriented 3-space
const α = [0.0, 0.0, 0.0]
const u = [0.0, -1.0, 0.0]
const v = [1.0, 0.0, 0.0]
const n = cross(u, v)      # (0, 0, 1)
const w = -u               # in-plane, opposite the bisector

# 3. Marks p_{j,σ} = α + a (c u + σ x_j v)
p(x_val, σ) = α + a * (c * u + σ * x_val * v)
const d₂ = p(x₂, -1.0)
const d₁ = p(x₁, -1.0)
const d  = p(x₀, -1.0)
const b  = p(x₀,  1.0)
const b₁ = p(x₁,  1.0)
const b₂ = p(x₂,  1.0)
const marks = [d₂, d₁, d, b, b₁, b₂]
const mark_names = ("d₂", "d₁", "d", "b", "b₁", "b₂")

# 4. Fundamental bearings and the complete 15-pair fan
const θ₀ = δ
const θ₁ = π / 4 + ϵ / 4
const θ₂ = atan((φ² + s) / c)
const γ  = θ₂ - θ₁
const β  = [-θ₂, -θ₁, -θ₀, θ₀, θ₁, θ₂]
const A_Ω = [abs(β[i] - β[j]) for i in 1:6, j in 1:6]
const angular_pairs = [(i, j, A_Ω[i, j]) for i in 1:5 for j in i+1:6]
@assert length(angular_pairs) == 15

# 5. Center-relative Gram: Ḡ = c² 1 1ᵀ + z zᵀ = U Uᵀ, U = [c 1 | z]
#    For general a, G_Ω = a² Ḡ. Here a = 1, so G_Ω = Ḡ.
#    Symbolic proof of rank 2: ∑zᵢ = 0 ⇒ UᵀU = diag(6c², ∑zᵢ²), both > 0.
const z = [-x₂, -x₁, -x₀, x₀, x₁, x₂]
const ones6 = ones(6)
const Ḡ = (c^2) * (ones6 * transpose(ones6)) + (z * transpose(z))
const G_Ω = (a^2) * Ḡ
const Ufac = [c * ones6 z]
const UtU = transpose(Ufac) * Ufac
const rank_G = rank(Ḡ)

# 6. Outer metric scaffold (Definition 11.1; adopted transcription)
const m_chord = α + a * c * u
const S_space = m_chord + a * φ² * w
const E_echo  = d₂ - a * φ² * w
const L_loco  = b₂ - a * φ² * w
const T_minus = d₂ + a * φ² * w
const T_plus  = b₂ + a * φ² * w
const scaffold_loop = [T_minus, T_plus, L_loco, E_echo, T_minus]

# Coordinate-based widths (not re-typed closed forms)
const W = norm(T_plus - T_minus)
const H = norm(T_minus - E_echo)
const L_chord = norm(b - d)
const scaffold_atol = 64 * eps(Float64) * max(W, 1.0)

# Closed forms, for comparison only. The identity check uses W, H, L_chord above.
const W_closed = 2 * a * (φ² + s)
const H_closed = 2 * a * φ²

# 7. Parametric kernel: major arc and golden feet (parameter ℓ)
const arc_t = range(δ, stop=2π - δ, length=500)
const arc_points = [α + a * (cos(t) * u + sin(t) * v) for t in arc_t]
const left_line  = [d - a * ℓ * v for ℓ in range(0.0, stop=φ², length=100)]
const right_line = [b + a * ℓ * v for ℓ in range(0.0, stop=φ², length=100)]
const rays = [[α, q] for q in marks]

# 8. Binary64 recoveries. These check the Float64 image, not the proofs.
println("--- OMEGA marked reference geometry (binary64, a = 1) ---")
@printf("rank(Ḡ) = %d  (symbolic proof: UᵀU = diag(6c², ∑ zᵢ²) with both entries > 0)\n", rank_G)
@printf("UᵀU[1,1] = %.16e   6c² = %.16e\n", UtU[1, 1], 6 * c^2)
@printf("UᵀU[2,2] = %.16e   ∑zᵢ² = %.16e\n", UtU[2, 2], sum(abs2, z))
@printf("UᵀU[1,2] = %.16e   (should be ~ 0 by ∑ zᵢ = 0)\n", UtU[1, 2])
@printf("G_Ω = a² Ḡ  with a = %.1f, so G_Ω = Ḡ\n", a)
@printf("θ₀ = %.12f arcsec\n", rad2deg(θ₀) * 3600)
@printf("θ₁ = %.12f arcsec  = 45°00′00.25″\n", rad2deg(θ₁) * 3600)
@printf("θ₂ = %.12f arcsec  = %s\n", rad2deg(θ₂) * 3600, "69°05′41.496848751…″")
@printf("γ  = %.12f arcsec\n", rad2deg(γ) * 3600)
@printf("15 unordered center angles: %d\n", length(angular_pairs))
@printf("constructed W = ‖T₊ − T₋‖ = %.16f\n", W)
@printf("constructed H = ‖T₋ − E‖  = %.16f\n", H)
@printf("constructed L_chord = ‖b − d‖ = %.16f\n", L_chord)
@printf("closed-form W = 2a(φ²+s) = %.16f\n", W_closed)
@printf("closed-form H = 2aφ²     = %.16f\n", H_closed)
@printf("W − H − L_chord = %.3e   (atol = %.3e)\n", W - H - L_chord, scaffold_atol)
@printf("constructed W ≈ H + L_chord: %s\n",
        string(isapprox(W, H + L_chord; rtol=0.0, atol=scaffold_atol)))
println()
println("Fifteen pairwise center angles |βᵢ − βⱼ| (degrees):")
for (i, j, ang) in angular_pairs
    @printf("  ∠ %s α %s = %.9f°\n", mark_names[i], mark_names[j], rad2deg(ang))
end

extract_xyz(pts) = ([p[1] for p in pts], [p[2] for p in pts], [p[3] for p in pts])

plotlyjs()
fig = plot(
    aspect_ratio = :equal,
    xlabel = "X (v)",
    ylabel = "Y (−u)",
    zlabel = "Z (n)",
    title = "OMEGA kernel, angular fan, and outer scaffold (binary64)",
)

for ray in rays
    plot!(fig, extract_xyz(ray)...; label="", color=:gray, lw=1, alpha=0.55)
end
plot!(fig, extract_xyz(arc_points)...; label="Major arc C(t)", color=:blue, lw=2)
plot!(fig, extract_xyz(left_line)...; label="Left foot", color=:black, lw=2)
plot!(fig, extract_xyz(right_line)...; label="Right foot", color=:black, lw=2)
plot!(fig, extract_xyz(scaffold_loop)...; label="Outer scaffold", color=:gray, linestyle=:dash, lw=1.4)
scatter!(fig, extract_xyz(marks)...; label="Baseline marks", color=:red, ms=4)
scatter!(fig, [α[1]], [α[2]], [α[3]]; label="Center α", color=:black, ms=5)
scatter!(fig, extract_xyz([S_space, E_echo, L_loco])...; label="Space, Echo, Loco", color=:orange, ms=5)

display(fig)