算額（その867）
奈良県橿原市木原町 耳成山口神社（天神社） 嘉永7年(1854)
http://www.wasan.jp/nara/miminasi.html
(3) 奈良県大和郡山市小泉町 耳成山口神社 嘉永7年(1854)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
牧下英世：数学教育を通して取り組んだ総合的な学習とその実証的な研究-算額を―算額を用いた課題学習とそのフィールドワークの実践から―，2002筑波大学附属駒場論集第42集
https://tsukuba.repo.nii.ac.jp/record/6466/files/13.pdf
外円内に甲円，乙円，丙円，丁円が互いに内接・外接して入っている。甲円，丁円の直径がそれぞれ 5 尺 2 寸，7 寸のとき，外円の直径はいかほどか。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (0, r2 - R)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (x4, y4)
とおき，以下の連立方程式を解く。
eq0 は算法助術の公式55 であるが，この問題では r2 = R - r1 なので eq1 のように簡単になる。
なお，デカルトの円定理による条件式 eq00 を使うこともできるが SymPy では解けない（数値解を求めることはできる）。
include("julia-source.txt");
using SymPy
@syms R::positive, r1::positive, r2::positive, r3::positive, r4::positive,
x3::positive, y3::negative, x4::positive, y4::negative,
d1::positive, d2::positive, d3::positive, d4::positive
# (k1, k2, k3, k4) = 1.0 ./ (r1, r2, r3, -R)
# eq00 = 1/(-(k1 + k2 + k3) + 2sqrt(k1*k2 + k2*k3 + k3*k1)) - R
(r1, r4) = (52, 7) .// 2
r2 = R - r1
(d1, d2, d3, d4) = (R, r1, r2, r3)
eq0 = d1^2*d2^2*d3^2 + d1^2*d4^2*(d2 - d3)^2 + 2*d1*d2^2*d3^2*d4 - 2*d1*d2*d3*d4*(d1 - d4)*(d2 + d3) + d2^2*d3^2*d4^2
eq1 = -R^2*r1 + R^2*r3 + R*r1^2 - R*r1*r3 + r1^2*r3
eq2 = x3^2 + y3^2 - (R - r3)^2 |> expand
eq3 = x4^2 + y4^2 - (R - r4)^2 |> expand
eq4 = x3^2 + (r2 - R - y3)^2 - (r2 + r3)^2 |> expand
eq5 = x4^2 + (r2 - R - y4)^2 - (r2 + r4)^2 |> expand
eq6 = (x3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2 |> expand
res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (R, r3, x3, y3, x4, y4))
1-element Vector{NTuple{6, Sym{PyCall.PyObject}}}:
(13 + 13*sqrt(785)/15, 728/73, 1456/73, -13 - 221*sqrt(785)/1095, 14, -19*sqrt(785)/30 - 13)
甲円，丁円の直径がそれぞれ 5 尺 2 寸，7 寸のとき，外円の直径は 2(13 + 13√785/15) = 74.56427585055592 寸である。
2(13 + 13√785/15)
74.56427585055592
その他のパラメータは以下のとおりである。
r1 = 26; r4 = 3.5; R = 37.2821; r2 = 11.2821; r3 = 9.9726; x3 = 19.9452; y3 = -18.6547; x4 = 14; y4 = -30.7446
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r4) = (26, 3.5)
(R, r3, x3, y3, x4, y4) = 13 .* (1 + √785/15, 56/73, 112/73, -1 - 17√785/1095, 14/13, -1 - 19√785/390)
r2 = R - r1
@printf("r1 = %g; r4 = %g; R = %g; r2 = %g; r3 = %g; x3 = %g; y3 = %g; x4 = %g; y4 = %g\n", r1, r4, R, r2, r3, x3, y3, x4, y4)
plot()
circle(0, 0, R)
circle(0, R - r1, r1, :green)
circle(0, r2 - R, r2, :magenta)
circle2(x3, y3, r3, :blue)
circle2(x4, y4, r4, :brown)
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, R - r1, "甲円:r1,(0,R-r1)", :green, :center, delta=-delta/2)
point(0, r2 - R, "乙円:r2,(0,r2-R)", :magenta, :center, delta=-delta/2)
point(x3, y3, "丙円:r3,(x3,y3)", :blue, :center, delta=-delta/2)
point(x4, y4, " 丁円:r4,(x4,y4)", :black, :left, :vcenter)
point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
end
end;
