算額（その710）
八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg
キーワード：円6個，楕円2個，正方形
#Julia, #SymPy, #算額, #和算
一辺の長さが 2.5 寸の正方形の中に楕円 2 個，楕円の中に大円 2 個，小円 1 個ずつが入っている。小円の直径が 0.4 寸のとき，大円の直径はいかほどか。
正方形の一辺の長さを 2a
大円の半径と中心座標を r1, (x1, a/2)
小円の半径と中心座標を r2, (0, a - r2)
楕円と大円の接点の座標を (x0, y0)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms a::positive, b::positive,
r1::positive, x1::positive,
r2::positive,
x0::positive, y0::positive
b = a/2
eq1 = x1^2 + (a/2 - r2)^2 - (r1 + r2)^2
eq2 = -b^2*x0/(a^2*(y0 - a/2)) + (x0 - x1)/(y0 - a/2)
eq3 = (x0 - x1)^2 + (y0 - a/2)^2 - r1^2
eq4 = x0^2/a^2 + (y0 -a/2)^2/b^2 - 1
res = solve([eq1, eq2, eq3, eq4], (r1, x1, x0, y0))
2-element Vector{NTuple{4, Sym}}:
(a/2 - r2/2, sqrt(3)*sqrt(r2)*sqrt(2*a - r2)/2, 2*sqrt(3)*sqrt(r2)*sqrt(2*a - r2)/3, a/2 - sqrt(9*a^2 - 24*a*r2 + 12*r2^2)/6)
(a/2 - r2/2, sqrt(3)*sqrt(r2)*sqrt(2*a - r2)/2, 2*sqrt(3)*sqrt(r2)*sqrt(2*a - r2)/3, a/2 + sqrt(9*a^2 - 24*a*r2 + 12*r2^2)/6)
2 組の解が得られるが，y0 の値が異なる（楕円の長軸に対して対称）だけで，両方とも適解である。
大円の半径は (a - r2)/2なので，直径は a - r2 である。
正方形の一辺の長さが 2.5 寸，小円の直径が 0.4 寸のとき，大円の直径は (2a - 2r2)/2 = (2.5 - 0.4)/2 = 1.05 である。
術は，「正方形の一辺の長さから小円の直径を引き，半分にする」 (2.5 - 0.4)/2 = 1.05
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, r2) = (25, 4) .// 20
b = a/2
(r1, x1, x0, y0) = (a/2 - r2/2, sqrt(3)*sqrt(r2)*sqrt(2*a - r2)/2, 2*sqrt(3)*sqrt(r2)*sqrt(2*a - r2)/3, a/2 + sqrt(9*a^2 - 24*a*r2 + 12*r2^2)/6)
plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:blue, lw=0.5)
circle4(x1, a/2, r1)
circle(0, a - r2, r2, :orange)
circle(0, r2 - a, r2, :orange)
ellipse(0, a/2, a, b)
ellipse(0, -a/2, a, b)
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(x0, y0, " (x0,y0)", :green, :left, :bottom, delta=delta/2)
point(0, a, " a", :blue, :left, :bottom, delta=delta/2)
point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
point(x1, a/2, " 大円:r1,(x1,a/2)", :red, :center, delta=-delta/2)
point(0, a - r2, " 小円:r2,(0,a-r2)", :black, :center, :bottom, delta=delta/2)
end
end;
