算額（その1458）
七〇 加須市大字外野 棘脱地蔵堂
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，外円，半円1個，楕円2個
#Julia, #SymPy, #算額, #和算
外円の中に外円と同じ直径の円弧を描き，甲円 1 個，乙円 2 個，等楕円 2 個を容れる。甲円の直径は外円の直径の半分であり，楕円と外円および楕円と弧は 1 点で接する。乙円の直径が 1 寸のとき，楕円の長径の最大値はいかほどか。
注：円弧は外円の中心を通る。また，外円（および円弧）は楕円の曲率円である。
外円の半径と中心座標を R, (0, 0)
弧の半径と中心座標を R, (0, -R)
甲円の半径と中心座標を r1, (0, R - r1); R = 2r1
乙円の半径と中心座標を r2, (x2, y2)
楕円の長径と短径を 2a, 2b; b = r1/2; 曲率円であるから a = sqrt(b*R)
とおき，以下の連立方程式を解く。
R = 2r1, b = r1/2 のとき， R = a^2/b ゆえ，a = sqrt(2r1*r1/2) = r1 である。
よって，「楕円の長径 = 甲円の直径」である。
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,
x2::positive, y2::positive,
a::positive, b::positive
R = 2r1
b = r1/2
eq1 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + (R + y2)^2 - (R + r2)^2
eq3 = x2^2 + y2^2 - (R - r2)^2
res = solve([eq1, eq2, eq3], (r1, x2, y2))[1]
(5*r2/3, 4*sqrt(3)*r2/3, r2/3)
r1 = 5r2/3である。
よって，楕円の長径は 2r1 = 10r2/3 で，r2 = 1/2 寸のとき 2r1 = 5/3 ≒ 1.67 寸。
以下は，図を描くためにパラメータを求める。
R = 2r1 = 10r2/3
b = 5r2/6 である。
蛇足
外円が楕円の曲率円であるから R = a^2/b より，
a = sqrt(b*R) = sqrt(5r2/6 * 10r2/3) = sqrt(25r2^2/9) = r2*5/3
楕円の長径 = 2a = 2r2*5/3 = 乙円の直径 * 5/3 ≒ 乙円の直径 * 1.67
乙円の直径が 1 寸のとき，楕円の長径は 5/3 寸 ≒ 1.67 寸である。甲円の直径と等しい。
function draw(r2, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, x2, y2) = (5*r2/3, 4*sqrt(3)*r2/3, r2/3)
R = 2r1
b = r1/2 # R/4
a = sqrt(b*R) # 曲率円: R = a^2/b
@printf("乙円の直径が %g のとき，楕円の長径は %g である。\n", 2r2, 2a)
@printf("r2 = %g; r1 = %g; x2 = %g; y2 = %g; R = %g; a = %g; b = %g\n",
r2, r1, x2, y2, R, a, b)
plot()
circle(0, 0, R, :magenta)
circle(0, -R, R, beginangle=30, endangle=150)
circle(0, R - r1, r1, :green)
circle2(x2, y2, r2)
ellipse(0, -R/4, a, R/4, color=:blue)
ellipse(0, -3R/4, a, R/4, color=:blue)
if more == true
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, "R", :green, :center, :bottom, delta=delta)
point(0, R - r1, "甲円:r1,(0,R-r1)", :green, :center, delta=-delta)
point(x2, y2, "乙円:r2,(x2,y2)", :red, :center, delta=-delta)
point(0, -R/4, "楕円:a,b,(0,-R/4)", :blue, :center, delta=-delta)
point(0, -3R/4, "楕円:a,b,(0,-3R/4)", :blue, :center, delta=-delta)
#plot!(xlims=(-0.2, 0.2), ylims=(-0.1, 0.05))
end
end;
draw(1/2, true)
