算額（その349）
岐阜県大垣市赤坂町 金生山明星輪寺 元治2年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF
台形内に青の半円と，赤，白，黄，浅青，緑の合計 7 個の円を入れる。赤円の直径が与えられたときに，緑円の直径を求めよ。
台形の下底（図の右側の辺），上底の長さをそれぞれ a, b とする。青円の中心を座標原点とする。。
青円の半径と中心座標 r1, (0, 0)
赤円の半径と中心座標 r2, (r2, r2)
黃円の半径と中心座標 r4, (0, r1 - r4)
緑円の半径と中心座標 r6, (x6, y6)
白円の半径と中心座標 r3, (r1 - r3, y3)
浅青円の半径と中心座標 r5, (r5 - r1, y5)
r2 を既知とすれば，残りの未知数は 11 変数なので，11 元連立方程式になるが，条件は 10 個しかない。a, b は確定できないので，b も与えないと解けない。いずれにせよ solve() では一度には解けないので，白円と浅青円を除いてまず解を求め，しかる後に得られた解を既知として描画に必要なパラメータを求める。
include("julia-source.txt");
using SymPy
@syms r1::positive, r2::positive, r3::positive, y3::positive, r4::positive,
r5::positive, y5::positive, r6::positive, x6::positive, y6::positive;
r2 = 1
eq1 = 2r2^2 - (r1 - r2)^2
eq3 = r2^2 + (r1 - r4 - r2)^2 - (r2 + r4)^2
eq5 = (x6 - r2)^2 + (y6 - r2)^2 - (r2 + r6)^2
eq6 = x6^2 + (y6 - r1 + r4)^2 - (r4 + r6)^2
eq7 = x6^2 + y6^2 - (r1 - r6)^2
res = solve([eq1, eq3, eq5, eq6, eq7], (r1, r4, r6, x6, y6))
1-element Vector{NTuple{5, Sym}}:
(1 + sqrt(2), -1 + sqrt(2), 1/5, 3/5, 1 + 4*sqrt(2)/5)
赤円の半径を 1 としたとき，緑円の半径は 1/5 である（赤円の直径の 1/5 である）。
残る 6 変数 r3, y3, r5, y5, a, b に対して，条件式は 5 個しかない。
図を見ればわかるが，a と b は互いに従属しており，a が小さくなれば b は大きくなる。
ここでは，b = 2 として解を求める。
using SymPy
@syms r1::positive, r2::positive, r3::positive, y3::positive, r4::positive,
r5::positive, y5::positive, r6::positive, x6::positive, y6::positive,
a::positive, b::positive;
r2 = 1
b = 2
(r1, r4, r6, x6, y6) = (1 + sqrt(Sym(2)), -1 + sqrt(Sym(2)), 1//5, 3//5, 1 + 4*sqrt(Sym(2))/5)
eq2 = (r1 - r3)^2 + y3^2 - (r1 + r3)^2
eq4 = (r1 - r5)^2 + y5^2 - (r1 + r5)^2
eq8 = r3/(a - y3) - r1/a
eq9 = r5/(b - y5) - r1/b
eq10 = 4r1^2 + (a - b)^2 - (a + b)^2
res = solve([eq2, eq4, eq8, eq9, eq10], (r3, y3, r5, y5, a))
1-element Vector{NTuple{5, Sym}}:
(-66932*sqrt(236*sqrt(2) + 334) - 94628*sqrt(118*sqrt(2) + 167) - 7 + 9*sqrt(2) + 66912*sqrt(2)*sqrt(118*sqrt(2) + 167) + 47328*sqrt(2)*sqrt(236*sqrt(2) + 334), -4 + 2*sqrt(118*sqrt(2) + 167)/(2*sqrt(2) + 3), (-sqrt(2)/2 - 1/2)*(-3 - 2*sqrt(2) + sqrt(20*sqrt(2) + 29)) + 1 + sqrt(2), -3 - 2*sqrt(2) + sqrt(20*sqrt(2) + 29), sqrt(2) + 3/2)
二段階に分けて求めたパラメータで図を描き，解が正しいことを確認する。
r1 = 2.41421; r2 = 1; r4 = 0.414214
r6 = 0.2; x6 = 0.6; y6 = 2.13137
r3 = 0.533631; y3 = 2.27006
r5 = 0.313594; y5 = 1.74021
a = 2.91421; b = 2
using Plots
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 1
(r1, r4, r6, x6, y6) = (1 + sqrt(2), -1 + sqrt(2), 1/5, 3/5, 1 + 4*sqrt(2)/5)
@printf("r1 = %g; r2 = %g; r4 = %g\n", r1, r2, r4)
@printf("r6 = %g; x6 = %g; y6 = %g\n", r6, x6, y6)
b = 2
(r3, y3, r5, y5, a) = (-66932*sqrt(236*sqrt(2) + 334) - 94628*sqrt(118*sqrt(2) + 167) - 7 + 9*sqrt(2) + 66912*sqrt(2)*sqrt(118*sqrt(2) + 167) + 47328*sqrt(2)*sqrt(236*sqrt(2) + 334), -4 + 2*sqrt(118*sqrt(2) + 167)/(2*sqrt(2) + 3), (-sqrt(2)/2 - 1/2)*(-3 - 2*sqrt(2) + sqrt(20*sqrt(2) + 29)) + 1 + sqrt(2), -3 - 2*sqrt(2) + sqrt(20*sqrt(2) + 29), sqrt(2) + 3/2)
@printf("r3 = %g; y3 = %g\n", r3, y3)
@printf("r5 = %g; y5 = %g\n", r5, y5)
@printf("a = %g; b = %g\n", a, b)
plot([r1, r1, -r1, -r1, r1], [0, a, b, 0, 0], color=:black, lw=0.5)
circlef(0, 0, r1, :dodgerblue, beginangle=0, endangle=180)
circlef(r2, r2, r2, :firebrick1)
circlef(-r2, r2, r2, :firebrick1)
circlef(0, r1 - r4, r4, :yellow)
circlef(x6, y6, r6, :limegreen)
circlef(-x6, y6, r6, :limegreen)
circle(r1 - r3, y3, r3, :gray)
circlef(r5 -r1, y5, r5, :skyblue1)
if more
point(r1, 0, " r1", :black, :left, :bottom)
point(r1, a, "(r1,a)", :black, :right, :bottom)
point(-r1, b, "(-r1,b)\n", :black, :left, :bottom)
point(r2, r2, "(r2,r2)", :black)
point(0, r1 - r4, " r1-r4", :black)
point(x6, y6, "(x6,y6)", :black)
point(r1 - r3, y3, "(r1-r3,y3)", :black)
point(r5 - r1, y5, "(r5-r1,y5)", :black)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
plot!(xlims=(-2.5, 2.7))
else
plot!(showaxis=false)
end
end;
