算額（その1308）
百三十四 群馬県富岡市一ノ宮 貫前神社 明治20年(1887)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，直角三角形
#Julia, #SymPy, #算額, #和算
直角三角形に内接する大円は等円の直径の計算には無関係なお飾りである。大円を除けば，算額（その928）の第一象限の図と同じになる。
直角三角形の 3 辺を，「鈎」，「股」，「弦」
大円の半径と中心座標を r0, (r0, r0)
等円の半径と中心座標を r1, (r1, y1), (x1, r1), (x, y); x = (r1 + x1)/2, y = (y1 + r1)/2
とおき，以下の連立方程式を解く。
1. 大円の半径 r0 は以下の式で求められる。
include("julia-source.txt");
using SymPy
@syms r0::positive, r1::positive,
x1::positive, y1::positive,
x::positive, y::positive,
鈎::positive, 股::positive, 弦::positive;
eq0 = 鈎 + 股 - sqrt(鈎^2 + 股^2) ⩵ 2r0
ans_r0 = solve(eq0, r0)[1]
ans_r0 |> println
ans_r0(鈎 => 5.4, 股 => 7.2) |> println
股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2
1.80000000000000
2. r1,x1,y1 は SymPy の能力的に，同時には解けないので，逐次解いていく。
eq1 = dist2(0, 鈎, 股, 0, r1, y1, r1)
eq2 = dist2(0, 鈎, 股, 0, x1, r1, r1)
eq3 = (x1 - r1)^2 + (y1 - r1)^2 - (4r1)^2;
まず，eq1, eq2 から x1, y1 を求める。
(ans_x1, ans_y1) = solve([eq1, eq2], (x1, y1))[1] # 1 of 4
(-r1*sqrt(股^2 + 鈎^2)/鈎 - 股*(r1 - 鈎)/鈎, -r1*sqrt(股^2 + 鈎^2)/股 - 鈎*(r1 - 股)/股)
eq3 の x1, y1 に代入し，r1 を求める。
eq13 = eq3(x1 => ans_x1, y1 => ans_y1)
ans_r1 = solve(eq13, r1)[2] # 2 of 2
ans_r1 |> println
股*鈎*(股^3 + 股^2*鈎 + 股^2*sqrt(股^2 + 鈎^2) + 股*鈎^2 - 4*股*鈎*sqrt(股^2 + 鈎^2) + 鈎^3 + 鈎^2*sqrt(股^2 + 鈎^2))/(2*(股^4 + 股^3*鈎 + 股^3*sqrt(股^2 + 鈎^2) - 6*股^2*鈎^2 + 股^2*鈎*sqrt(股^2 + 鈎^2) + 股*鈎^3 + 股*鈎^2*sqrt(股^2 + 鈎^2) + 鈎^4 + 鈎^3*sqrt(股^2 + 鈎^2)))
鈎，股が 5.4 寸, 7.2 寸のとき，等円の直径は 2 寸である。
ans_r1(鈎 => 5.4, 股 => 7.2) * 2 |> println
2.00000000000000
sqrt(股^2 + 鈎^2) が何回も出てくるので，「弦」で置き換える。
ans_r1 = ans_r1(sqrt(股^2 + 鈎^2) => 弦) |> simplify
ans_r1 |> println
股*鈎*(弦*股^2 - 4*弦*股*鈎 + 弦*鈎^2 + 股^3 + 股^2*鈎 + 股*鈎^2 + 鈎^3)/(2*(弦*股^3 + 弦*股^2*鈎 + 弦*股*鈎^2 + 弦*鈎^3 + 股^4 + 股^3*鈎 - 6*股^2*鈎^2 + 股*鈎^3 + 鈎^4))
すべてのパラメータは以下のとおりである。
鈎 = 5.4; 股 = 7.2; 弦 = 9; r0 = 1.8; r1 = 1; x1 = 4.2; y1 = 3.4; x = 2.6; y = 2.2
function draw(鈎, 股, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
弦 = sqrt(鈎^2 + 股^2)
r0 = 股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2
r1 = 股*鈎*(弦*股^2 - 4*弦*股*鈎 + 弦*鈎^2 + 股^3 + 股^2*鈎 + 股*鈎^2 + 鈎^3)/(2*(弦*股^3 + 弦*股^2*鈎 + 弦*股*鈎^2 + 弦*鈎^3 + 股^4 + 股^3*鈎 - 6*股^2*鈎^2 + 股*鈎^3 + 鈎^4))
(x1, y1) = (-r1*sqrt(股^2 + 鈎^2)/鈎 - 股*(r1 - 鈎)/鈎, -r1*sqrt(股^2 + 鈎^2)/股 - 鈎*(r1 - 股)/股)
(x, y) = ((r1 + x1)/2, (y1 + r1)/2)
@printf("鈎，股が %g, %g のとき，等円の直径は %g である。\n", 鈎, 股, 2r1)
@printf("鈎 = %g; 股 = %g; 弦 = %g; r0 = %g; r1 = %g; x1 = %g; y1 = %g; x = %g; y = %g\n", 鈎, 股, 弦, r0, r1, x1, y1, x, y)
plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:green, lw=0.5)
circle(r0, r0, r0, :blue)
circle(r1, y1, r1)
circle(x, y, r1)
circle(x1, r1, r1)
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(0, 鈎, "鈎", :green, :left, :bottom, delta=delta/2, deltax=delta/2)
point(股, 0, "股", :green, :left, :bottom, delta=delta/2, deltax=delta/2)
point(r0, r0, "大円:r0,(r0,r0)", :blue, :center, delta=-delta)
point(r1, y1, "等円:r1,(r1,y1)", :red, :center, delta=-delta)
point(x, y, "等円:r1,(x,r)", :red, :center, delta=-delta)
point(x1, r1, "等円:r1,(x1,r1)", :red, :center, delta=-delta)
end
end;
draw(5.4, 7.2, true)
#SymPy, #算額, #和算
