算額（その403）
東京都あきる野市二宮 二宮神社 寛政6年(1794)
山口正義：やまぶき，第15号
https://yamabukiwasan.sakura.ne.jp/ymbk15.html
鈎股弦の中に内接する全円と 5 個の等円がある。全円の直径が 44寸，弦が 110 寸のとき，等円の直径はいかほどか。
全円の半径と中心座標を r1, (r1, r1)
等円の半径を r2 として以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms 鈎::positive, 股::positive, 弦::positive, r1::positive,
r2::positive, x::positive, y::positive;
eq1 = 鈎 + 股 - 弦 ~ 2r1
eq2 = 鈎^2 + 股^2 ~ 弦^2
eq3 = 弦 ~ 股 - x + 鈎 - y + 8r2
eq4 = r1/(鈎 - r1) ~ r2/(鈎 - y)
eq5 = r1/(股 - r1) ~ r2/(股 - x)
res = solve((eq1, eq2, eq3, eq4, eq5), (鈎, 股, r2, x, y))
2-element Vector{NTuple{5, Sym}}:
(r1 + 弦/2 - sqrt(-4*r1^2 - 4*r1*弦 + 弦^2)/2, r1 + 弦/2 + sqrt(-4*r1^2 - 4*r1*弦 + 弦^2)/2, r1*弦/(8*r1 + 弦), r1*(8*r1 + 5*弦 + 4*sqrt(-4*r1^2 - 4*r1*弦 + 弦^2))/(8*r1 + 弦), r1*(8*r1 + 5*弦 - 4*sqrt(-4*r1^2 - 4*r1*弦 + 弦^2))/(8*r1 + 弦))
(r1 + 弦/2 + sqrt(-4*r1^2 - 4*r1*弦 + 弦^2)/2, r1 + 弦/2 - sqrt(-4*r1^2 - 4*r1*弦 + 弦^2)/2, r1*弦/(8*r1 + 弦), r1*(8*r1 + 5*弦 - 4*sqrt(-4*r1^2 - 4*r1*弦 + 弦^2))/(8*r1 + 弦), r1*(8*r1 + 5*弦 + 4*sqrt(-4*r1^2 - 4*r1*弦 + 弦^2))/(8*r1 + 弦))
最初の解が適解である。
等円の直径は 2r1*弦/(8*r1 + 弦) である。
r1 = 44//2, 弦 = 110 のとき，等円の直径は 16.923077 である。
鈎 = 66; 股 = 88; r2 = 8.46154; x = 62.6154; y = 49.0769
等円の直径 = 16.923077
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, 弦) = (44//2, 110)
(鈎, 股, r2, x, y) = (r1 + 弦/2 - sqrt(-4*r1^2 - 4*r1*弦 + 弦^2)/2, r1 + 弦/2 + sqrt(-4*r1^2 - 4*r1*弦 + 弦^2)/2, r1*弦/(8*r1 + 弦), r1*(8*r1 + 5*弦 + 4*sqrt(-4*r1^2 - 4*r1*弦 + 弦^2))/(8*r1 + 弦), r1*(8*r1 + 5*弦 - 4*sqrt(-4*r1^2 - 4*r1*弦 + 弦^2))/(8*r1 + 弦))
@printf("鈎 = %g; 股 = %g; r2 = %g; x = %g; y = %g\n",鈎, 股, r2, x, y)
@printf("等円の直径 = %.8g\n", 2r2)
plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
circle(r1, r1, r1)
#circle(x, r2, r2, :blue)
#circle(r2, y, r2, :blue)
for i = 0:4
circle(x - i*2r2*股/弦, r2 + i*2r2*鈎/弦, r2, :blue)
end
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(r1, r1, "全円:r1,(r1,r1)", :red, :center, delta=-delta)
point(x, r2, "(x,r2)", :blue)
point(r2, y, "(r2,y)", :blue)
point(0, 鈎, "鈎 ", :black, :right, :vcenter)
point(股, 0, "股", :black, :left, :bottom, delta=delta)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
plot!(xlims=(-5, 90))
else
plot!(showaxis=false)
end
end;
