算額（その1543）
千葉県富津市寺尾 六所神社 明治4年(1871)
山口正義：やまぶき 第 57 号
https://yamabukiwasan.sakura.ne.jp/ymbk57.pdf
山口正義：やまぶき 第 59 号
https://yamabukiwasan.sakura.ne.jp/ymbk59.pdf
キーワード：球9個，3次元
#Julia, #SymPy, #算額, #和算
外球の中に甲球 1 個，乙球 7 個を容れる。外球の直径が 12 寸，甲球の直径が 8 寸のとき，乙球の直径はいかほどか。
外球の半径と中心座標を R, (0, 0, 0)
甲球の半径と中心座標を r1, (0, 0, R - r1)
乙球の半径と中心座標を r2, (x2, y2, z2)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# # julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
@syms R, r1, r2, x2, y2, z2
eq1 = x2^2 + (R - r1 - z2)^2 - (r1 + r2)^2
eq2 = x2^2 + z2^2 - (R - r2)^2
eq3 = x2*sind(Sym(360)/14) - r2
res = solve([eq1, eq2, eq3], (r2, x2, z2))[2] # 2 of 2
(2*R*r1*(-R*cos(2*pi/7) + R - r1 + r1*cos(2*pi/7))/(R^2 - 2*R*r1*cos(2*pi/7) + r1^2), 2*R*r1*(-R*cos(2*pi/7) + R - r1 + r1*cos(2*pi/7))/((R^2 - 2*R*r1*cos(2*pi/7) + r1^2)*sin(pi/7)), R*(R^2 - 2*R*r1 - r1^2 + 2*r1^2*cos(2*pi/7))/(R^2 - 2*R*r1*cos(2*pi/7) + r1^2))
# 2r2
2res[1](R => 12/2, r1 => 8/2).evalf() |> println
3.27511575021090
外球，甲球の直径がそれぞれ 12 寸，8 寸のとき，乙球の直径は 3.27511575021090 である。
答・術は 4.7071 寸としているが，山口も 3.275 有奇として異を唱えている。
function draw(R, r1, more=false)
pyplot(size=(500, 250), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, x2, z2) = (2*R*r1*(-R*cos(2*pi/7) + R - r1 + r1*cos(2*pi/7))/(R^2 - 2*R*r1*cos(2*pi/7) + r1^2), 2*R*r1*(-R*cos(2*pi/7) + R - r1 + r1*cos(2*pi/7))/((R^2 - 2*R*r1*cos(2*pi/7) + r1^2)*sin(pi/7)), R*(R^2 - 2*R*r1 - r1^2 + 2*r1^2*cos(2*pi/7))/(R^2 - 2*R*r1*cos(2*pi/7) + r1^2))
p1 = plot()
circle(0, 0, R, :cyan)
circle(0, 0, x2, :blue)
rotate(x2, 0, r2, angle=360/7)
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, R, "R", :black, :center, :bottom, delta=delta/2)
point(x2, 0, "乙球:r2,(x2,0,z2)", :red, :right, delta=-delta/2, deltax=15delta)
end
p2 = plot()
circle(0, 0, R, :cyan)
circle(0, R- r1, r1, :blue)
circle(x2, z2, r2)
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, R, "R", :black, :center, :bottom, delta=delta/2)
point(0, R - r1, "甲球:(0,0, R-r1)", :blue, :center, delta=-delta/2)
point(x2, z2, "乙球:r2,(x2,0,z2)", :red, :right, delta=-delta/2, deltax=15delta)
end
plot(p1, p2)
end;
draw(12/2, 8/2, true)
