算額（その1651）
三十一 一関市舞川 観福寺内地蔵堂前額 明治34年(1901)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
今有如図 03085
https://w.atwiki.jp/sangaku/pages/157.html
キーワード：円1個，外円，楕円2個，斜線2本
#Julia, #SymPy, #算額, #和算
全円の中に等楕円 2 個，斜線 2 本を容れる。上側の楕円は短径で外円と一点で接する（外円は楕円の曲率円である）。楕円の長径，短径が与えられたとき，斜の長さを求める術を述べよ。
外円の半径と中心座標を R, (0, 0)
楕円の長半径，短半径と中心座標を a, b, (0, R - b), (0, R - 3b)
下側の楕円と斜線の接点座標を (x0, y0)
斜線と全円の交点座標を (x, -sqrt(R^2 - x^2))
とおき，以下の連立方程式を解く。
斜 = sqrt(x^2 + (R + sqrt(R^2 - x^2))^2) も，連立方程式に組み入れておく。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms R::positive, a::positive, b::positive, x0::positive, y0::positive, x::positive, 斜::positive
eq1 = R - a^2/b
eq2 = x0^2/a^2 + (y0 - (R - 3b))^2/b^2 - 1
eq3 = -b^2*x0/(a^2*(y0 - (R - 3b)))- (y0 - R)/x0
eq4 = (y0 - R)/x0 + (R + sqrt(R^2 - x^2))/x
eq5 = 斜 - sqrt(x^2 + (R + sqrt(R^2 - x^2))^2)
res = solve([eq1, eq2, eq3, eq4, eq5], (斜, R, x0, y0, x))[1]
(a^2*sqrt(2*a^2 + 16*b^2 + 2*sqrt(a^4 - 16*a^2*b^2 + 64*b^4))/(b*sqrt(a^2 + 8*b^2)), a^2/b, 2*sqrt(2)*a/3, a^2/b - 8*b/3, 4*sqrt(2)*a^3/(a^2 + 8*b^2))
# 斜
res[1]
2sqrt(a^4 - 16a^2*b^2 + 64b^4) = -2(a^2 - 8b^2) を代入すると，SymPy が自動的に簡約化してくれる。
斜 = (a^2*sqrt(2*a^2 + 16*b^2 - 2*(a^2 - 8*b^2))/(b*sqrt(a^2 + 8*b^2)))
たとえば，長半径 = 5，短半径 = 2 のとき，斜線の長さは 18.7317162316339 である。
斜(a => 5, b => 2).evalf() |> println
18.7317162316339
術は，以下のようになっており，同じ答えを得る。
長径 = 10
短径 = 4
法 = sqrt(長径^2 + 8短径^2)
斜 = √8長径^2/法
18.73171623163388
長径・短径を与えるか長半径・短半径を与えるかの違いだけで，術と同じ解が得られた。
function draw(a, b, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(斜, R, x0, y0, x) = (a^2*sqrt(2*a^2 + 16*b^2 + 2*sqrt(a^4 - 16*a^2*b^2 + 64*b^4))/(b*sqrt(a^2 + 8*b^2)), a^2/b, 2*sqrt(2)*a/3, a^2/b - 8*b/3, 4*sqrt(2)*a^3/(a^2 + 8*b^2))
@printf("a = %g; b = %g; 斜 = %g; R = %g; x0 = %g; y0 = %g; x = %g; ", a, b, 斜, R, x0, y0, x)
plot()
circle(0, 0, R)
ellipse(0, R - b, a, b, color=:blue)
ellipse(0, R - 3b, a, b, color=:blue)
segment(0, R, x, -sqrt(R^2 - x^2))
segment(0, R, -x, -sqrt(R^2 - x^2))
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)", :blue, :left, :vcenter)
point(x, -sqrt(R^2 - x^2), "(x,-sqrt(R^2-x^2)) ", :red, :right, :vcenter)
point(0, R - 3b, " R-3b", :blue, :left, :vcenter)
point(0, R, "R", :red, :center, :bottom, delta=delta/2)
end
end;
draw(5, 2, true)
