glen/dyna3
1
0
Fork 0
Code Issues 11 Pull Requests 1 Actions Projects Releases 2 Wiki Activity

Clean up example of three mutually tangent spheres

Browse Source
This commit is contained in:
Aaron Fenyes 2024-02-15 17:16:37 -08:00
parent f2000e5731
commit 3170a933e4
2 changed files with 4 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
engine-proto/Engine.Algebraic.jl
View File

@ -198,7 +198,6 @@ function realize(ctx::Construction{T}) where T
end
(Generic.Ideal(coordring, eqns), eqns)
## [test] (nothing, eqns)
end
end

67
engine-proto/Engine.jl
View File

@ -23,63 +23,19 @@ using GLMakie
CoeffType = Rational{Int64}
##a = Engine.Point{CoeffType}()
##s = Engine.Sphere{CoeffType}()
##a_on_s = Engine.LiesOn{CoeffType}(a, s)
##ctx = Engine.Construction{CoeffType}(elements = Set([a]), relations= Set([a_on_s]))
##ideal_a_s = Engine.realize(ctx)
##println("A point on a sphere: $(Engine.dimension(ideal_a_s)) degrees of freedom")
##b = Engine.Point{CoeffType}()
##b_on_s = Engine.LiesOn{CoeffType}(b, s)
##Engine.push!(ctx, b)
##Engine.push!(ctx, s)
##Engine.push!(ctx, b_on_s)
##ideal_ab_s, eqns_ab_s = Engine.realize(ctx)
##freedom = Engine.dimension(ideal_ab_s)
##println("Two points on a sphere: $freedom degrees of freedom")
spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
tangencies = [
Engine.AlignsWithBy{CoeffType}(
spheres[n],
spheres[mod1(n+1, length(spheres))],
CoeffType([1, 1, 1][n])
CoeffType(1)
)
for n in 1:3
]
##tangencies = [
##Engine.LiesOn{CoeffType}(points[1], spheres[2]),
##Engine.LiesOn{CoeffType}(points[1], spheres[3]),
##Engine.LiesOn{CoeffType}(points[2], spheres[3]),
##Engine.LiesOn{CoeffType}(points[2], spheres[1]),
##Engine.LiesOn{CoeffType}(points[3], spheres[1]),
##Engine.LiesOn{CoeffType}(points[3], spheres[2])
##]
ctx_tan_sph = Engine.Construction{CoeffType}(elements = spheres, relations = tangencies)
ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph)
##small_eqns_tan_sph = eqns_tan_sph
##small_eqns_tan_sph = [
## eqns_tan_sph;
## spheres[2].coords - [1, 0, 0, 0, 1];
## spheres[3].coords - [1, 0, 0, 0, -1];
##]
##small_ideal_tan_sph = Generic.Ideal(base_ring(ideal_tan_sph), small_eqns_tan_sph)
freedom = Engine.dimension(ideal_tan_sph)
println("Three mutually tangent spheres, with two fixed: $freedom degrees of freedom")
##points = [Engine.Point{CoeffType}() for _ in 1:3]
##spheres = [Engine.Sphere{CoeffType}() for _ in 1:2]
##ctx_joined = Engine.Construction{CoeffType}(
## elements = Set([points; spheres]),
## relations= Set([
## Engine.LiesOn{CoeffType}(pt, sph)
## for pt in points for sph in spheres
## ])
##)
##ideal_joined, eqns_joined = Engine.realize(ctx_joined)
##freedom = Engine.dimension(ideal_joined)
##println("$(length(points)) points on $(length(spheres)) spheres: $freedom degrees of freedom")
println("Three mutually tangent spheres: $freedom degrees of freedom")
# --- test rational cut ---
@ -90,7 +46,7 @@ vbls = Variable.(symbols(coordring))
system = CompiledSystem(System(eqns_tan_sph, variables = vbls))
norm2 = vec -> real(dot(conj.(vec), vec))
rng = MersenneTwister(6071)
n_planes = 36
n_planes = 6
samples = []
for _ in 1:n_planes
real_solns = solution.(Engine.Numerical.real_samples(system, freedom, rng = rng))
@ -100,22 +56,7 @@ for _ in 1:n_planes
end
end
end
println("$(length(samples)) sample solutions:")
##for soln in samples
## ## display([vbls round.(soln, digits = 6)]) ## [verbose]
## k_sq = abs2(soln[1])
## if abs2(soln[end-2]) > 1e-12
## if k_sq < 1e-12
## println(" center at infinity: z coordinates $(round(soln[end], digits = 6)) and $(round(soln[end-1], digits = 6))")
## else
## sum_sq = soln[4]^2 + soln[7]^2 + soln[end-2]^2 / k_sq
## println(" center on z axis: r² = $(round(1/k_sq, digits = 6)), x² + y² + h² = $(round(sum_sq, digits = 6))")
## end
## else
## sum_sq = sum(soln[[4, 7, 10]] .^ 2)
## println(" center at origin: r² = $(round(1/k_sq, digits = 6)); x² + y² + z² = $(round(sum_sq, digits = 6))")
## end
##end
println("Found $(length(samples)) sample solutions")
# show a sample solution
function show_solution(ctx, vals)