算額（その441）
東都下谷稲荷社 寛政元年己酉3月(1624)
黒川喜兵衛盛榮
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
円弧の中に甲円 1 個，乙円 2 個，丙円 2 個があり，それぞれの円は隣同士は概説し，円弧を構成する外円には内接し，弦に接している。
甲円，丙円の直径がそれぞれ 12 寸，3 寸のとき，乙円の直径はいかほどか。
以下のように記号を定め方程式を解く。
外円の中心を原点とする。
外円の半径，中心座標を r0, (0, 0)
甲円の半径，中心座標を r1, (0, a + r1)
乙円の半径，中心座標を r2, (x2, a + r2)
丙円の半径，中心座標を r3, (x3, a + r3)
弦が y 軸と交わる座標を (0, a)
として，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms r0::positive, a::positive,
r1::positive, r2::positive, x2::positive,
r3::positive, x3::positive;
eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq3 = x2^2 + (a + r2)^2 - (r0 -r2)^2
eq4 = x3^2 + (a + r3)^2 - (r0 -r3)^2
eq5 = a + 2r1 - r0
res = solve([eq1, eq2, eq3, eq4, eq5], (a, r0, r2, x2, x3))
2-element Vector{NTuple{5, Sym}}:
((-2*r1^(3/2)*sqrt(r3) - r1^2 + 3*r1*r3)/(r1 - r3), (-2*r1^(3/2)*sqrt(r3) + r1^2 + r1*r3)/(r1 - r3), 2*(-r1^(7/2)*sqrt(r3) + 2*r1^(5/2)*r3^(3/2) - r1^(3/2)*r3^(5/2) - r1^3*r3 + 2*r1^2*r3^2 - r1*r3^3)/(r1^3 - 3*r1^2*r3 + 3*r1*r3^2 - r3^3), -2*sqrt(2)*r1^(3/2)*sqrt(r3)/sqrt(-r1^(3/2)*sqrt(r3) + r1*r3), sqrt(-8*r1^(3/2)*sqrt(r3) + 8*r1*r3))
((2*r1^(3/2)*sqrt(r3) - r1^2 + 3*r1*r3)/(r1 - r3), (2*r1^(3/2)*sqrt(r3) + r1^2 + r1*r3)/(r1 - r3), 2*(r1^(7/2)*sqrt(r3) - 2*r1^(5/2)*r3^(3/2) + r1^(3/2)*r3^(5/2) - r1^3*r3 + 2*r1^2*r3^2 - r1*r3^3)/(r1^3 - 3*r1^2*r3 + 3*r1*r3^2 - r3^3), 2*sqrt(2)*r1^(3/2)*sqrt(r3)/sqrt(r1^(3/2)*sqrt(r3) + r1*r3), sqrt(8*r1^(3/2)*sqrt(r3) + 8*r1*r3))
2 番めのものが適解である。
names = ("a", "r0", "r2", "x2", "x3")
for (i, name) in enumerate(names)
@printf("%s = %g; ", name, res[2][i](r1 => 6, r3 => 1.5).evalf())
end
a = 6; r0 = 18; r2 = 4; x2 = 9.79796; x3 = 14.6969;
甲円，丙円の直径が 12 寸，3 寸のとき，乙円の直径は 8 寸である。
using Plots
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r3) =(12, 3) ./ 2
(a, r0, r2, x2, x3) = ((2*r1^(3/2)*sqrt(r3) - r1^2 + 3*r1*r3)/(r1 - r3), (2*r1^(3/2)*sqrt(r3) + r1^2 + r1*r3)/(r1 - r3), 2*(r1^(7/2)*sqrt(r3) - 2*r1^(5/2)*r3^(3/2) + r1^(3/2)*r3^(5/2) - r1^3*r3 + 2*r1^2*r3^2 - r1*r3^3)/(r1^3 - 3*r1^2*r3 + 3*r1*r3^2 - r3^3), 2*sqrt(2)*r1^(3/2)*sqrt(r3)/sqrt(r1^(3/2)*sqrt(r3) + r1*r3), sqrt(8*r1^(3/2)*sqrt(r3) + 8*r1*r3))
@printf("r1 = %g; r3 = %g\n", r1, r3)
@printf("a = %g; r0 = %g; r2 = %g; x2 = %g; x3 = %g\n", a, r0, r2, x2, x3)
@printf("乙円の直径 = %g\n", 2r2)
plot()
circle(0, 0, r0, :green)
circle(0, a + r1, r1, :magenta)
circle(x2, a + r2, r2, :blue)
circle(-x2, a + r2, r2, :blue)
circle(x3, a + r3, r3, :red)
circle(-x3, a + r3, r3, :red)
b = sqrt(r0^2 - a^2)
segment(-b, a, b, a, :orange)
if more == true
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(0, r0, " r0", :green, :left, :bottom, delta=delta/2)
point(0, a, " a", :orange, :left, :top, delta=-delta/2)
point(0, a + r1, "甲円:r1,(0,a+r1)", :magenta, :center, delta=-delta/2)
point(x2, a + r2, "乙円:r2,(x2,a+r2)", :black, :center, delta=-delta/2)
point(x3, a + r3, "丙円:r3,(x3,a+r3)", :black, :center, delta=-delta/2)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
plot!(ylims=(-3, r0 + 0.5))
else
plot!(showaxis=false)
end
end;
