算額（その817）
宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf
楕円内に等円 2 個が入っている。楕円の短径が 1 のとき，原点を点対称とする円と楕円の 2 接点間の距離を求めよ。
条件が「楕円の短径」のみなので，解けるかどうか不安であったが，「楕円の長径」は無関係であることがわかった。
楕円の長半径と短半径，中心座標を a, b, (0, 0)
円の半径と中心座標を r, (r, 0)
円と，楕円の接点を (x0，y0)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, b::positive, x0::positive, y0::positive, r::positive
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0/(a^2*y0) + (x0 - r)/y0
eq3 = (x0 - r)^2 + y0^2 - r^2
res = solve([eq1, eq2, eq3], (x0, y0, r))
2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
(a*b/sqrt(a^2 - b^2), -b*sqrt((a^2 - 2*b^2)/(a - b))/sqrt(a + b), b*sqrt(a^2 - b^2)/a)
(a*b/sqrt(a^2 - b^2), b*sqrt((a^2 - 2*b^2)/(a - b))/sqrt(a + b), b*sqrt(a^2 - b^2)/a)
2 組の解が得られるが，符号の違いのみなのでどちらでもよい。
接点間の距離 diag は以下のようになり，点対称な接点間の距離は，楕円の短径のみで決まる。常に，楕円の短径の √2 倍である。
すなわち，問のように短径が 1 寸ならば，接点間の距離は 1.41421356... である。
@syms d
diag = 2sqrt(res[2][1]^2 + res[2][2]^2)
diag |> println
2*sqrt(a^2*b^2/(a^2 - b^2) + b^2*((a^2 - 2*b^2)/(a - b))/(a + b))
長い式は，以下のように簡約化される。
apart(diag, d) |> simplify |> println
2*sqrt(2)*b
b に特定の値を代入して数値解を求める。
apart(diag, d)(b => 1/2) |> simplify |> println
1.41421356237310
その他のパラメータは以下のとおり。
b = 1.234; x0 = 1.35383; y0 = 1.1012; r = 1.12477
function draw(more)
pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b) = (3, 1.234)
(x0, y0, r) = (a*b/sqrt(a^2 - b^2), b*sqrt((a^2 - 2*b^2)/(a - b))/sqrt(a + b), b*sqrt(a^2 - b^2)/a)
println("点対称な接点間の距離 = $(2sqrt(x0^2 + y0^2)) = $(2√2b)")
@printf("b = %g; x0 = %g; y0 = %g; r = %g\n", b, x0, y0, r)
plot()
ellipse(0, 0, a, b, color=:red)
circle2(r, 0, r, :blue)
segment(x0, y0, -x0, -y0, :green)
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, b, " b", :red, :left, :bottom, delta=delta/2)
point(x0, y0, "(x0,y0)", :red, :left, :bottom, delta=delta)
point(-x0, -y0, "(-x0,-y0)", :red, :left, :bottom, delta=delta)
#plot!(xlims=(-0.1,0.1), ylims=(-0.1,0.1))
end
end;
