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

Study mutually tangent spheres with two fixed

Browse Source
This commit is contained in:
Aaron Fenyes 2024-02-15 13:28:01 -08:00
parent e41bcc7e13
commit 291d5c8ff6
2 changed files with 62 additions and 55 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -186,19 +186,19 @@ function realize(ctx::Construction{T}) where T
# add relations to center, orient, and scale the construction # add relations to center, orient, and scale the construction
# [to do] the scaling constraint, as written, can be impossible to satisfy # [to do] the scaling constraint, as written, can be impossible to satisfy
# when all of the spheres have to go through the origin # when all of the spheres have to go through the origin
if !isempty(ctx.points) ##if !isempty(ctx.points)
append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3]) ## append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3])
end ##end
if !isempty(ctx.spheres) ##if !isempty(ctx.spheres)
append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4]) ## append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4])
end ##end
n_elts = length(ctx.points) + length(ctx.spheres) ##n_elts = length(ctx.points) + length(ctx.spheres)
if n_elts > 0 ##if n_elts > 0
push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts) ## push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts)
end ##end
## [test] (Generic.Ideal(coordring, eqns), eqns) (Generic.Ideal(coordring, eqns), eqns)
(nothing, eqns) ## [test] (nothing, eqns)
end end
end end

93
engine-proto/Engine.jl
View File

@ -39,15 +39,15 @@ CoeffType = Rational{Int64}
##freedom = Engine.dimension(ideal_ab_s) ##freedom = Engine.dimension(ideal_ab_s)
##println("Two points on a sphere: $freedom degrees of freedom") ##println("Two points on a sphere: $freedom degrees of freedom")
##spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
##tangencies = [ tangencies = [
## Engine.AlignsWithBy{CoeffType}( Engine.AlignsWithBy{CoeffType}(
## spheres[n], spheres[n],
## spheres[mod1(n+1, length(spheres))], spheres[mod1(n+1, length(spheres))],
## CoeffType(-1//1) CoeffType(-1)^n
## ) )
## for n in 1:3 for n in 1:3
##] ]
##tangencies = [ ##tangencies = [
##Engine.LiesOn{CoeffType}(points[1], spheres[2]), ##Engine.LiesOn{CoeffType}(points[1], spheres[2]),
##Engine.LiesOn{CoeffType}(points[1], spheres[3]), ##Engine.LiesOn{CoeffType}(points[1], spheres[3]),
@ -56,34 +56,41 @@ CoeffType = Rational{Int64}
##Engine.LiesOn{CoeffType}(points[3], spheres[1]), ##Engine.LiesOn{CoeffType}(points[3], spheres[1]),
##Engine.LiesOn{CoeffType}(points[3], spheres[2]) ##Engine.LiesOn{CoeffType}(points[3], spheres[2])
##] ##]
##ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies))
##ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph) ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph)
##freedom = Engine.dimension(ideal_tan_sph) ##small_eqns_tan_sph = eqns_tan_sph
##println("Three mutually tangent spheres: $freedom degrees of freedom") 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(small_ideal_tan_sph)
println("Three mutually tangent spheres, with two fixed: $freedom degrees of freedom")
points = [Engine.Point{CoeffType}() for _ in 1:3] ##points = [Engine.Point{CoeffType}() for _ in 1:3]
spheres = [Engine.Sphere{CoeffType}() for _ in 1:2] ##spheres = [Engine.Sphere{CoeffType}() for _ in 1:2]
ctx_joined = Engine.Construction{CoeffType}( ##ctx_joined = Engine.Construction{CoeffType}(
elements = Set([points; spheres]), ## elements = Set([points; spheres]),
relations= Set([ ## relations= Set([
Engine.LiesOn{CoeffType}(pt, sph) ## Engine.LiesOn{CoeffType}(pt, sph)
for pt in points for sph in spheres ## for pt in points for sph in spheres
]) ## ])
) ##)
ideal_joined, eqns_joined = Engine.realize(ctx_joined) ##ideal_joined, eqns_joined = Engine.realize(ctx_joined)
freedom = Engine.dimension(ideal_joined) ##freedom = Engine.dimension(ideal_joined)
println("$(length(points)) points on $(length(spheres)) spheres: $freedom degrees of freedom") ##println("$(length(points)) points on $(length(spheres)) spheres: $freedom degrees of freedom")
# --- test rational cut --- # --- test rational cut ---
coordring = base_ring(ideal_joined) coordring = base_ring(small_ideal_tan_sph)
vbls = Variable.(symbols(coordring)) vbls = Variable.(symbols(coordring))
# test a random witness set # test a random witness set
system = CompiledSystem(System(eqns_joined, variables = vbls)) system = CompiledSystem(System(small_eqns_tan_sph, variables = vbls))
norm2 = vec -> real(dot(conj.(vec), vec)) norm2 = vec -> real(dot(conj.(vec), vec))
Random.seed!(6071) Random.seed!(6071)
n_planes = 3 n_planes = 36
samples = [] samples = []
for _ in 1:n_planes for _ in 1:n_planes
real_solns = solution.(Engine.Numerical.real_samples(system, freedom)) real_solns = solution.(Engine.Numerical.real_samples(system, freedom))
@ -94,21 +101,21 @@ for _ in 1:n_planes
end end
end end
println("$(length(samples)) sample solutions:") println("$(length(samples)) sample solutions:")
for soln in samples ##for soln in samples
## display([vbls round.(soln, digits = 6)]) ## [verbose] ## ## display([vbls round.(soln, digits = 6)]) ## [verbose]
k_sq = abs2(soln[1]) ## k_sq = abs2(soln[1])
if abs2(soln[end-2]) > 1e-12 ## if abs2(soln[end-2]) > 1e-12
if k_sq < 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))") ## println(" center at infinity: z coordinates $(round(soln[end], digits = 6)) and $(round(soln[end-1], digits = 6))")
else ## else
sum_sq = soln[4]^2 + soln[7]^2 + soln[end-2]^2 / k_sq ## 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))") ## println(" center on z axis: r² = $(round(1/k_sq, digits = 6)), x² + y² + h² = $(round(sum_sq, digits = 6))")
end ## end
else ## else
sum_sq = sum(soln[[4, 7, 10]] .^ 2) ## 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))") ## println(" center at origin: r² = $(round(1/k_sq, digits = 6)); x² + y² + z² = $(round(sum_sq, digits = 6))")
end ## end
end ##end
# show a sample solution # show a sample solution
function show_solution(ctx, vals) function show_solution(ctx, vals)