# OMEGA local geometric object — copy-pasteable Float64 illustration
# of the existing Paper VII arc-and-feet construction. Not a theorem checker.
# Not a new source paper. LinearAlgebra and Printf only.
# λ = φ/Ω₀ and μ = π/Ω₀ are selected-measurement scalars, not mark coordinates.

using LinearAlgebra
using Printf

# ============================================================
# OMEGA local geometric object
# Copy-pasteable Julia script
# ============================================================

# ------------------------------------------------------------
# 1. Defining constants
# ------------------------------------------------------------

const a = 1.0                          # radius / global scale
const φ = (1 + sqrt(5.0)) / 2         # golden ratio
const Ω₀ = φ + π

const ε = π / 648000.0                # one arcsecond in radians
const δ = ε / 2.0

const s = sin(δ)
const c = cos(δ)

const λ = φ / Ω₀
const μ = π / Ω₀

# ------------------------------------------------------------
# 2. Canonical placement in R^3
# ------------------------------------------------------------

const α = [0.0, 0.0, 0.0]             # center
const u = [0.0, -1.0, 0.0]            # downward radial axis
const v = [1.0,  0.0, 0.0]            # horizontal axis
const n = cross(u, v)                 # plane normal

@assert isapprox(dot(u, v), 0.0; atol=1e-14)
@assert isapprox(norm(u), 1.0; atol=1e-14)
@assert isapprox(norm(v), 1.0; atol=1e-14)
@assert isapprox(norm(n), 1.0; atol=1e-14)

# ------------------------------------------------------------
# 3. Positive baseline levels
#
# x0 = s
# x1 = 1 + s
# x2 = φ^2 + s
# ------------------------------------------------------------

const x0 = s
const x1 = 1.0 + s
const x2 = φ^2 + s

@assert isapprox(1.0 + φ, φ^2; atol=1e-14)

# ------------------------------------------------------------
# 4. Marked points on the baseline y = -ac, z = 0
#
# p(level, σ) = α + a*(c*u + σ*level*v)
# where σ = -1 or +1
# ------------------------------------------------------------

mark(level, σ) = α + a * (c * u + σ * level * v)

const p0m = mark(x0, -1)
const p0p = mark(x0, +1)

const p1m = mark(x1, -1)
const p1p = mark(x1, +1)

const p2m = mark(x2, -1)
const p2p = mark(x2, +1)

# Ordered left-to-right list
const marked_points = [
    ("p2,-", p2m),
    ("p1,-", p1m),
    ("p0,-", p0m),
    ("p0,+", p0p),
    ("p1,+", p1p),
    ("p2,+", p2p),
]

# ------------------------------------------------------------
# 5. Major arc
#
# γ(t) = (a sin t, -a cos t, 0)
# for ε/2 ≤ t ≤ 2π - ε/2
#
# This is the retained large arc of the circle.
# The tiny bottom arc between the seed points is omitted.
# ------------------------------------------------------------

gamma(t) = [a * sin(t), -a * cos(t), 0.0]

const arc_start = ε / 2.0
const arc_end   = 2.0 * π - ε / 2.0

# ------------------------------------------------------------
# 6. Horizontal extensions
#
# E_+ = {(a r, -a c, 0) : s ≤ r ≤ s + φ^2}
# E_- = {(-a r, -a c, 0) : s ≤ r ≤ s + φ^2}
# ------------------------------------------------------------

E_plus(r)  = [ a * r, -a * c, 0.0]
E_minus(r) = [-a * r, -a * c, 0.0]

# ------------------------------------------------------------
# 7. Technical checks
# ------------------------------------------------------------

# Seed points lie on the circle of radius a
@assert isapprox(norm(p0m - α), a; atol=1e-14)
@assert isapprox(norm(p0p - α), a; atol=1e-14)

# Their angular separation is exactly ε
function angle_between(x, y)
    cx = dot(x, y) / (norm(x) * norm(y))
    cx = clamp(cx, -1.0, 1.0)
    return acos(cx)
end

const θ_seed = angle_between(p0m - α, p0p - α)
# acos of a cosine 1-O(ε²) is ill-conditioned in Float64; chord formula is tight.
@assert isapprox(θ_seed, ε; atol=1e-10)
@assert isapprox(2 * asin(norm(p0p - p0m) / (2 * a)), ε; atol=1e-14)

# All six marked points lie in the plane span{u,v}
for (_, p) in marked_points
    @assert abs(dot(p - α, n)) < 1e-12
end

# Reflection across the vertical axis x = 0
reflect_x(q) = [-q[1], q[2], q[3]]

@assert isapprox(reflect_x(p0p), p0m; atol=1e-12)
@assert isapprox(reflect_x(p1p), p1m; atol=1e-12)
@assert isapprox(reflect_x(p2p), p2m; atol=1e-12)

# Consecutive baseline separations, left to right
const sep_21 = norm(p1m - p2m)   # should be aφ
const sep_10 = norm(p0m - p1m)   # should be a
const sep_00 = norm(p0p - p0m)   # should be 2as
const sep_01 = norm(p1p - p0p)   # should be a
const sep_12 = norm(p2p - p1p)   # should be aφ

@assert isapprox(sep_21, a * φ; atol=1e-12)
@assert isapprox(sep_10, a; atol=1e-12)
@assert isapprox(sep_00, 2.0 * a * s; atol=1e-12)
@assert isapprox(sep_01, a; atol=1e-12)
@assert isapprox(sep_12, a * φ; atol=1e-12)

# Center-relative radii of the marked points
ρ0 = norm(p0p - α)   # same as p0m
ρ1 = norm(p1p - α)   # same as p1m
ρ2 = norm(p2p - α)   # same as p2m

# Closed-form checks
@assert isapprox(ρ0, a; atol=1e-12)
@assert isapprox(ρ1, a * sqrt(2.0 * (1.0 + s)); atol=1e-12)
@assert isapprox(ρ2, a * φ * sqrt(3.0 + 2.0 * s); atol=1e-12)

# Arc endpoints are the seed points; samples lie on the circle.
@assert isapprox(gamma(arc_start), p0p; atol=1e-12)
@assert isapprox(gamma(arc_end), p0m; atol=1e-12)
@assert isapprox(E_plus(s), p0p; atol=1e-12)
@assert isapprox(E_minus(s), p0m; atol=1e-12)
@assert isapprox(E_plus(x2), p2p; atol=1e-12)
@assert isapprox(E_minus(x2), p2m; atol=1e-12)

# ------------------------------------------------------------
# 8. Optional samplers for numerical rendering
# ------------------------------------------------------------

function sample_major_arc(N::Int=2000)
    ts = range(arc_start, arc_end; length=N)
    return [gamma(t) for t in ts]
end

function sample_extension_plus(N::Int=400)
    rs = range(s, x2; length=N)
    return [E_plus(r) for r in rs]
end

function sample_extension_minus(N::Int=400)
    rs = range(s, x2; length=N)
    return [E_minus(r) for r in rs]
end

const arc_samples = sample_major_arc(722)
for q in arc_samples
    @assert isapprox(norm(q - α), a; atol=1e-12)
end
@assert isapprox(arc_samples[1], p0p; atol=1e-12)
@assert isapprox(arc_samples[end], p0m; atol=1e-12)

# ------------------------------------------------------------
# 9. Summary printout
# ------------------------------------------------------------

println("============================================================")
println("OMEGA local geometric object")
println("============================================================")

@printf("a            = %.15f\n", a)
@printf("φ            = %.15f\n", φ)
@printf("Ω₀ = φ + π   = %.15f\n", Ω₀)
println()

@printf("ε            = %.15e rad\n", ε)
@printf("δ = ε/2      = %.15e rad\n", δ)
@printf("s = sin(δ)   = %.15e\n", s)
@printf("c = cos(δ)   = %.15f\n", c)
println()

println("Positive baseline levels:")
@printf("x0 = s       = %.15f\n", x0)
@printf("x1 = 1 + s   = %.15f\n", x1)
@printf("x2 = φ^2 + s = %.15f\n", x2)
println()

println("Six marked points (left to right):")
for (name, p) in marked_points
    @printf("%-4s = (% .15f, % .15f, % .15f)\n", name, p[1], p[2], p[3])
end
println()

println("Consecutive separations along the baseline:")
@printf("|p2,- - p1,-| = %.15f   (expected aφ)\n", sep_21)
@printf("|p1,- - p0,-| = %.15f   (expected a)\n",  sep_10)
@printf("|p0,+ - p0,-| = %.15e   (expected 2as)\n", sep_00)
@printf("|p1,+ - p0,+| = %.15f   (expected a)\n",  sep_01)
@printf("|p2,+ - p1,+| = %.15f   (expected aφ)\n", sep_12)
println()

println("Center-relative radii:")
@printf("ρ0 = %.15f\n", ρ0)
@printf("ρ1 = %.15f\n", ρ1)
@printf("ρ2 = %.15f\n", ρ2)
println()

println("Seed angle check:")
@printf("θ_seed       = %.15e rad\n", θ_seed)
@printf("difference   = %.15e\n", abs(θ_seed - ε))
println()

println("Normalized scalar shares:")
@printf("λ = φ/(φ+π)  = %.15f\n", λ)
@printf("μ = π/(φ+π)  = %.15f\n", μ)
@printf("λ + μ        = %.15f\n", λ + μ)
println()

println("Geometric definition:")
println("  Major arc: γ(t) = (a sin t, -a cos t, 0),   ε/2 ≤ t ≤ 2π - ε/2")
println("  Extensions: E±(r) = (±a r, -a c, 0),        s ≤ r ≤ s + φ^2")
println("  Object = major arc ∪ left extension ∪ right extension")
println("============================================================")
println("LOCAL OBJECT PASS: seed pair at radius a; major arc ∪ E₋ ∪ E₊;")
println("consecutive separations aφ, a, 2as, a, aφ; not a diameter figure.")