算額（その1624）
長野県上水内牟礼村牟礼 渋薬師堂嘉永2年(1849) 大久保善賢氏保管
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
キーワード：球5個，回転楕円体，3次元
#Julia, #Julia, #SymPy, #算額, #和算, #数学
回転楕円体の中に，甲球 2 個，乙球 3 個を容れる。回転楕円体の長径と短径が与えられたとき，甲球の直径を得る術を述べよ。
回転楕円体の長径，短径を 2a, 2b
甲球の半径と中心座標を r1, (0, 0, z1)
乙球の半径と中心座標を r2, (b - r2, 0, 0)
甲球と回転楕円体の x-z 平面上の交点座標を (x0, 0, z0)
とおき，以下の連立方程式を解く。
まず，eq1 を解き r2 を求める。
次いで， eq3, eq4, eq5 を解き，z1, x0, z0 を求める。
最後に，eq2 を解き r1 を求める。
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,
z0::positive, r1::positive, z1::positive,
r2::positive
eq1 = (2r2/√Sym(3) + r2) - b
eq2 = z1^2 + (b - r2)^2 - (r1 + r2)^2 |> expand
eq3 = x0^2/a^2 + z0^2/b^2 - 1 |> expand
eq4 = (x0 - z1)^2 + z0^2 - r1^2 |> expand
eq5 = -b^2*x0/(a^2*z0) + (x0 - z1)/z0 |> expand;
# r2 乙球の半径
ans_r2 = solve(eq1, r2)[1] |> simplify
ans_r2 |> println
b*(-3 + 2*sqrt(3))
# z1, x0, z0
res = solve([eq3, eq4, eq5], (z1, x0, z0))[4] # 4 of 4
(sqrt((a*b^2 - a*r1^2 + b^3 - b*r1^2)/(a - b))*(a - b)/b, a^2*sqrt((b^2 - r1^2)/(a - b))/(b*sqrt(a + b)), sqrt((a^2*r1^2 - b^4)/(a - b))/sqrt(a + b))
eq2 に，ここまでで得られた r2, x1, x0, z0 を代入し，r1 を求める。
eq12 = eq2(r2 => ans_r2, z1 => res[1], x0 => res[2], z0 => res[3]) |> simplify;
# r1 甲球の半径
ans_r1 = solve(eq12, r1)[1]
ans_r1 |> println
b*(a^2 - 4*sqrt(3)*b^2 + 6*b^2)/a^2
甲球の半径は b*(a^2 - 4√3b^2 + 6b^2)/a^2 である。
長径が 4，短径が 2 のとき，甲球の直径は 1.53589838486225 である。
2*ans_r1(a => 4/2, b => 2/2).evalf()|> println
1.53589838486225
術は以下のようになっている。
@syms 短径, 長径
甲球径 = 2(短径/2 - (2√3 -3)*短径^3/長径^2)
甲球径(短径 => 2, 長径 => 4) |> println
1.53589838486225
すべてのパラメータは以下の通りである。
回転楕円体の長径が 4，短径が 2 のとき，甲球の直径は 1.5359，乙球の直径は 0.928203 である。
a = 2; b = 1; r1 = 0.767949; r2 = 0.464102; z1 = 1.1094; x0 = 1.4792; z0 = 0.673049
function draw1(a, b, r1, z1, r2, x0, z0, more)
p1 = plot()
circle(0, 0, b, :green)
rotate(0, b - r2, r2, :blue)
circle(0, 0, r1)
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(0, b, "b", :green, :center, :bottom, delta=delta/2)
point(0, b - r2, "乙球:r2(0,b-r2,0)", :blue, :center, delta=-delta/2)
point(0, r1, "r1", :red, :center, :bottom, delta=delta/2)
end
end;
function draw2(a, b, r1, z1, r2, x0, z0, more)
plot()
circle2(z1, 0, r1)
circle(0, b - r2, r2, :blue)
circle(0, -(b - r2)/2, r2, :blue)
ellipse(0, 0, a, b, color=:green)
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(0, b - r2, "乙球:r2\n(0,b-r2,0)", :blue, :center, delta=-delta)
point(a, 0, "a", :green, :left, :bottom, delta=delta)
point(0, b, "b", :green, :center, :bottom, delta=delta)
point(z1, 0, "甲球:r1,(0, 0, z1)", :red, :center, delta=-delta)
point(x0, z0, "(x0,z0)", :green, :left, :bottom, delta=delta)
end
end;
function draw(a, b, func, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = b*(a^2 - 4*sqrt(3)*b^2 + 6*b^2)/a^2
r2 = b*(-3 + 2*sqrt(3))
z1 = sqrt((a*b^2 - a*r1^2 + b^3 - b*r1^2)/(a - b))*(a - b)/b
x0 = a^2*sqrt((b^2 - r1^2)/(a - b))/(b*sqrt(a + b))
z0 = sqrt((a^2*r1^2 - b^4)/(a - b))/sqrt(a + b)
@printf("回転楕円体の長径が %g，短径が %g のとき，甲球の直径は %g，乙球の直径は %g である。\n", 2a, 2b, 2r1, 2r2)
@printf("a = %g; b = %g; r1 = %g; r2 = %g; z1 = %g; x0 = %g; z0 = %g\n", a, b, r1, r2, z1, x0, z0)
func(a, b, r1, z1, r2, x0, z0, more)
end;
draw(4/2, 2/2, draw1, true)
draw(4/2, 2/2, draw2, true)
