算額（その791）
寛政十一年己未四月 丸山良玄門人 北越中宿邑 米持杢左衛門富房
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf
右鈎，左鈎，容円の直径が 9 寸，5 寸，3 寸のとき，雙股はいかほどか。
右鈎，左鈎，雙股 を R, L, a
容円の半径と中心座標を r, (x, r)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms x::positive, a::positive, d,
r::positive, L::positive, R::positive
eq1 = dist(0, L, a, 0, x, r) - r^2
eq2 = dist(a, R, 0, 0, x, r) - r^2;
eq1 = numerator(apart(eq1, d))
eq2 = numerator(apart(eq2, d))
eq1 |> println
eq2 |> println
L^2*a^2 - 2*L^2*a*x - L^2*r^2 + L^2*x^2 - 2*L*a^2*r + 2*L*a*r*x
-R^2*r^2 + R^2*x^2 - 2*R*a*r*x
2 組の解が得られるが，最初のものが適解である。
右鈎，左鈎，容円の直径が 9 寸，5 寸，3 寸のとき，雙股は 112 寸である。
res = solve([eq1, eq2], (a, x));
res[1][1] |> println
res[1][1](L => 5, R => 9, r => 3//2) |> println
2*r*sqrt((L - 2*r)/(L*R^2 - 4*L*R*r + 4*L*r^2 - 2*R^2*r + 4*R*r^2))*(-L*R + L*r + R*r)/(-L + 2*r)
12
容円の中心座標は (9/2, 3/2) である。
res[1][2] |> println
res[1][2](L => 5, R => 9, r => 3//2) |> println
R*r*sqrt((L - 2*r)/(L*R^2 - 4*L*R*r + 4*L*r^2 - 2*R^2*r + 4*R*r^2))
9/2
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(L, R, r) = (5, 9, 3/2)
(a, x) = (2*r*sqrt((L - 2*r)/(L*R^2 - 4*L*R*r + 4*L*r^2 - 2*R^2*r + 4*R*r^2))*(-L*R + L*r + R*r)/(-L + 2*r), R*r*sqrt((L - 2*r)/(L*R^2 - 4*L*R*r + 4*L*r^2 - 2*R^2*r + 4*R*r^2)))
@printf("右鈎 = %g; 左鈎 = %g; 容円直径 = %g\n", R, L, 2r)
@printf("雙股 = %g; x = %g\n", a, x)
plot([0, 0, a, 0, a, a, 0], [0, L, 0, 0, R, 0], color=:blue, lw=0.5)
circle(x, r, r)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(a, R, "(a,R)", :blue, :right, :bottom, delta=delta/2)
point(0, L, " R", :blue, :left, :bottom, delta=delta/2)
point(x, r, "容円:r,(x,r)", :red, :center, delta=-delta/2)
point(0, L/2, " 左鈎", :blue, :left, :vcenter, mark=false)
point(a, R/2, "右鈎 ", :blue, :right, :vcenter, mark=false)
point(a/2, 0, "雙股", :blue, :center, :bottom, delta=delta/2, mark=false)
end
end;
