算額（その966）
一七 大里郡岡部村岡 稲荷社 文化13年(1816)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
直角三角形の中に大小の正三角形と甲円，乙円が入っている。乙円の直径が 153 寸のとき，甲円の直径はいかほどか。
大きい正三角形の一辺の長さを a
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (r2, y2)
とする。
乙円，甲円が入っている二等辺三角形の相似比が 1:√3 であることは関係式を少し書き下せばわかり，甲円の直径も乙円の直径の √3 倍であることが導ける。
ここでは，図を各パラメータを全て求めるため，以下の連立方程式を解く。
SymPy の能力的に一度に解けないので，逐次解いてゆく。
include("julia-source.txt");
using SymPy
@syms a::positive, r1::positive, x1::positive,
r2::positive, y2::positive
@syms a, r1, x1, r2, y2
eq1 = dist2(0, 0, a/2, √Sym(3)a/2, r2, y2, r2)
eq2 = dist2(a, 0, a/2, √Sym(3)a/2, x1, r1, r1)
eq3 = dist2(0, a/√Sym(3), a/2, √Sym(3)a/2, r2, y2, r2)
eq4 = dist2(2a, 0, a/2, √Sym(3)a/2, x1, r1, r1);
# res = solve([eq1, eq2, eq3, eq4], (r1, x1, y2, a))
y2 を求める。
ans_y2 = solve(eq1, y2)[2]
ans_y2 |> println
r2*(sqrt(3) + 2)
eq3 に y2 を代入し，a を求める。
eq3 = eq3(y2 => ans_y2) |> simplify;
ans_a = solve(eq3, a)[2]
ans_a |> println
2*r2*(sqrt(3) + 2)
eq2, eq4 に a を代入し，eq2 から x1 を求める。
eq2 = eq2(a => ans_a) |> simplify;
eq4 = eq4(a => ans_a) |> simplify;
ans_x1 = solve(eq2, x1)[2]
ans_x1 |> println
sqrt(3)*r1/3 + 2*sqrt(3)*r2 + 4*r2
eq4 に x1 を代入すれば r1, r2 だけを含む式になるので, r1 を求める（r2 を含む式になる）。
eq4 = eq4(x1 => ans_x1) |> simplify
eq4 |> println
r1^2/3 - 8*sqrt(3)*r1*r2/3 - 4*r1*r2 + 4*sqrt(3)*r2^2 + 7*r2^2
ans_r1 = solve(eq4, r1)[2] |> sympy.sqrtdenest |> simplify
ans_r1 |> println
sqrt(3)*r2
既知の r2 により r1 が求まったので，順次遡って以下のパラメータを得る。
y2 = r2*(√3 + 2)
a = 2y2
r1 = √3r2
x1 = √3r1/3 + a
甲円の半径 r1 は乙円の半径 r2 の √3 倍である。
乙円の直径が 153 寸のとき，甲円の直径は 153√3 = 265.00377355803823 寸である。
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = 153/2
y2 = r2*(√3 + 2)
a = 2y2
r1 = √3r2
x1 = √3r1/3 + a
@printf("乙円の直径が %g のとき，甲円の直径は %g である。\n", 2r2, 2r1)
@printf("r2 = %g; r1 = %g; x1 = %g; y2 = %g; a = %g\n", r2, r1, x1, y2, a)
plot([0, 2a, 0, 0], [0, 0, 2a/√3, 0], color=:blue, lw=0.5)
circle(x1, r1, r1)
circle(r2, y2, r2, :green)
segment(0, a/√3, a/2, √3a/2, :magenta)
segment(0, 0, a/2, √3a/2, :magenta)
segment(a, 0, a/2, √3a/2, :magenta)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:black, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(x1, r1, "甲円:r1\n(x1,r1)", :red, :center, delta=-delta)
point(r2, y2, " 乙円:r2,(r2,y2)", :green, :left, :vcenter)
point(a, 0, "a ", :magenta, :right, :bottom, delta=delta)
point(a/2, √3a/2, "(a/2,√3a/2)", :magenta, :left, :bottom, delta=delta)
end
end;
