算額（その1425）
八十九 陸前高田市小友町 常膳寺観音堂 天保13年(1842)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円6個，外円，弦3本
#Julia, #SymPy, #算額, #和算
全円の中に 3 本の弦を引き，区画された領域に甲円 3 個，乙円 2 個を容れる。乙円の直径が 1 寸のとき，甲円の直径はいかほどか。
全円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1), (x1, R - 3r1)
乙円の半径と中心座標を r2, (x2, R - 2r1 + r2)
とおき，以下の連立方程式を解く。
注1：左右にある乙円と甲円の x 座標値は同じになる。x2 = x1
注2：eq1, eq2 は甲円，乙円が全円に内接するという条件を記述したものである。
注3：コメントアウトした eq3, eq4は斜線と甲円，乙円の中心の距離を記述したものであるが，連立方程式を解くことができないほど複雑になる。図に灰色で示した 2 つの円は甲円と乙円と合同なものである。この位置にある甲円，乙円についての式が有効な eq3, eq4 である。これを使えば，連立方程式は容易に解くことができる
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms R::positive, r1::positive, r2::positive, x1::positive, x2, x0::positive
x2 = x1
eq1 = x1^2 + (R - 3r1)^2 - (R - r1)^2
eq2 = x2^2 + (R - 2r1 + r2)^2 - (R - r2)^2
# eq3 = dist2(0, -R, x0, sqrt(R^2 - x0^2), x1, R - 3r1, r1)
# eq4 = dist2(0, -R, x0, sqrt(R^2 - x0^2), x2, R - 2r1 + r2, r2)
eq3 = 2R*r1 - sqrt(x0^2 + (R - sqrt(R^2 - x0^2))^2)*(2R - r1)
eq4 = r2*(2R - r1) - r1*R
res = solve([eq1, eq2, eq3, eq4], (R, r1, x1, x0))[1];
# R
res[1] |> println
r2*(2 + sqrt(5))
# r1
res[2] |> println
r2*(1 + sqrt(5))/2
# x1
res[3] |> println
r2*sqrt(2 + 2*sqrt(5))
# x0
res[4] |> sympy.sqrtdenest |> println
2*r2*sqrt(-8 + 4*sqrt(5))
甲円の半径 r1 は，乙円の半径 r2 の (1 + √5)/2 倍である。
乙円の直径が 1 寸のとき，甲円の直径は (1 + √5)/2 = 1.618033988749895 寸である。
全てのパラメータは以下のとおりである
r2 = 0.5; R = 2.11803; r1 = 0.809017; x1 = 1.27202; x0 = 0.971737; y0 = 1.88197; y = 0.5
function draw(r2, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
R = r2*(√5 + 2)
r1 = r2*(√5 + 1)/2
x1 = r2*sqrt(2√5 + 2)
x0 = 2r2*sqrt(4√5 - 8)
x2 = x1
y0 = sqrt(R^2 - x0^2)
y = R - 2r1
@printf("乙円の直径が %g のとき，甲円の直径は %g である。\n", 2r2, 2r1)
@printf("r2 = %g; R = %g; r1 = %g; x1 = %g; x0 = %g; y0 = %g; y = %g\n", r2, R, r1, x1, x0, y0, y)
plot([x0, 0, -x0], [y0, -R, y0], color=:green, lw=0.5)
circle(0, 0, R)
circle(0, R - r1, r1)
circle2(x1, y - r1, r1)
circle2(x2, y + r2, r2, :blue)
circle(0, y - r2, r2, :gray80)
circle(0, r1 - R, r1, :gray80)
segment(sqrt(R^2 - y^2), y, -sqrt(R^2 - y^2), y, :magenta)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
vline!([0], color=:gray80, lw=0.5)
hline!([0], color=:gray80, lw=0.5)
point(0, R, "R", :magenta, :center, :bottom, delta=delta/2)
point(x1, y + r2, "乙円:r2\n(x2,y+r2)", :blue, :center, delta=-delta/2)
point(0, R - r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta/2)
point(x1, y - r1, "甲円:r1,(x1,y-r1)", :red, :center, delta=-delta/2)
point(x0, sqrt(R^2 - x0^2), "(x0,sqrt(R^2-x0^2) ", :green, :left, :bottom, delta=delta/2)
point(0, y, " y", :magenta, :left, :bottom, delta=delta/2)
end
end;
draw(1/2, true)
