算額（その1617）
山形県鶴岡市大山 椙尾神社 文政元年8月
http://www.wasan.jp/yamagata/sugio.html
キーワード：円2個，円弧，直角三角形
#Julia, #SymPy, #算額, #和算, #数学
外円の一部の円弧（弓形）の中に楕円を容れる。楕円の長径が 6 寸，短径が 3 寸，矢が 3.1 寸のとき外円の直径はいかほどか。
外円の半径と中心座標を R, (0, 0)
楕円の半径と中心座標を a, b, (0, R - 矢 + b)
円弧と楕円の接点座標を (x0, y0)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms R::positive, x0::positive, y0::positive, a::positive, b::positive, 矢::positive
eq1 = x0^2/a^2 + (y0 - (R - 矢 + b))^2/b^2 - 1
eq2 = x0^2 + y0^2 - R^2
eq3 = -b^2*x0/(a^2*(y0 - (R - 矢 + b))) + x0/y0;
res = solve([eq2, eq3], (x0, y0))[2];
# x0
@syms d
ans_x0 = res[1] |> simplify |> x -> apart(x, d) |> simplify
ans_x0 |> println
sqrt(-2*R^2*a^2*b^2 + R^2*b^4 - 2*R*a^4*b + 2*R*a^4*矢 - a^4*b^2 + 2*a^4*b*矢 - a^4*矢^2)/(a^2 - b^2)
# y0
ans_y0 = res[2]
ans_y0 |> println
a^2*(R + b - 矢)/(a^2 - b^2)
eq11 = eq1(x0 => ans_x0, y0 => ans_y0) |> simplify |> numerator;
# R
ans_R = solve(eq11, R)[1] # 1 of 2
ans_R |> println
a*(a*(-b + 矢) - sqrt(矢)*sqrt(-2*a^2*b + a^2*矢 + 2*b^3 - b^2*矢))/b^2
長径が 6，短径が 3，矢が 3.1 のとき，外円の直径は 8.94254 である。
2ans_R(a => 6/2, b => 3/2, 矢 => 3.1).evalf() |> println
8.94253969560281
function draw(a, b, 矢, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
R = a*(a*(-b + 矢) - sqrt(矢)*sqrt(-2*a^2*b + a^2*矢 + 2*b^3 - b^2*矢))/b^2
x0 = sqrt(-(R*b^2 + a^2*b - a^2*矢)*(2*R*a^2 - R*b^2 + a^2*b - a^2*矢))/((a - b)*(a + b))
y0 = a^2*(R + b - 矢)/(a^2 - b^2)
y = R - 矢
x = sqrt(R^2 - y^2)
θ = atand(y, x)
@printf("長径が %g，短径が %g，矢が %g のとき，外円の直径は %g である。\n", 2a, 2b, 矢, 2R)
@printf("x0 = %g; y0 = %g\n", x0, y0)
plot()
circle(0, 0, R, :blue, beginangle=θ, endangle=180-θ, n=1000)
ellipse(0, R - 矢 + b, a, b, color=:red)
segment(-x, y, x, y, :blue)
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)
dimension_line(0, R, 0, y, "矢 ", :blue, :right, dx=-3delta)
dimension_line(0, y + b, 0, y + 2b, " b", :red, :left)
dimension_line(0, y + b, a, y + b, " a", :red, :center, delta=-2delta)
point(0, R, "R", :blue, :center, :bottom, delta=2delta)
point(0, y, "R-矢+b", :red, :center, delta=-2delta)
point(x0, y0, "(x0,y0)", :red, :left, :bottom)
point(0, y + b, "", :red)
point(a, y + b, "", :red)
point(0, y + 2b, "", :red)
end
end;
draw(6/2, 3/2, 3.1, true)
