算額（その1661） - 算額あれこれ

長野県上田市別所温泉北向観音堂文化4年(1807)

中村信弥「改訂増補長野県の算額」県内の算額(P.84)
http://www.wasan.jp/zoho/zoho.html

岐阜県大垣市西外側町大垣八幡神社（大垣八幡宮）天保年間(1831~1845)

http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

キーワード：3次元，球，外球
#Julia #SymPy #算額 #和算 #数学

大球の中に甲球 1 個，乙球 4 個，丙球 1 個を容れる。甲球の直径が 6 寸，丙球の直径が 2 寸のとき，乙球の直径はいかほどか。

大球の半径と中心座標を \(R,\ (0,\ 0,\ 0)\)
甲球の半径と中心座標を \(r_1,\ (0,\ 0,\ R - r_1)\)
乙球の半径と中心座標を \(r_2,\ (x_2,\ 0,\ z_2);\ x_2 = \sqrt{2}r_2\)
丙球の半径と中心座標を \(r_3,\ (0,\ 0,\ r_3 -R)\)
とおき，以下の連立方程式を解く。

include("julia-source.txt"); # julia-source.txt ソース

using SymPy
@syms r1::positive, r2::positive, z2::negative,
r3::positive, R::positive
x2 = √Sym(2)r2
eq1 = x2^2 + (R - r1 - z2)^2 - (r1 + r2)^2 # 甲球と乙球が外接する
eq2 = x2^2 + (z2 - r3 + R)^2 - (r2 + r3)^2 # 乙球と丙球が外接する
eq3 = x2^2 + z2^2 - (R - r2)^2 # 乙球は大球に内接する
res = solve([eq1, eq2, eq3], (r2, z2, R))[2]

(2*r1*r3/(r1 + r3), (-r1^2 - r1*sqrt(r1^2 + 6*r1*r3 + r3^2) + r3^2 + r3*sqrt(r1^2 + 6*r1*r3 + r3^2))/(2*(r1 + r3)), r1/2 + r3/2 + sqrt(r1^2 + 6*r1*r3 + r3^2)/2)

# r2 乙球の半径
res[1] |> println
res[1](r1 => 6/2, r3 => 2/2) |> println

2*r1*r3/(r1 + r3)
1.50000000000000

乙球の半径は \(r_2 = 2r_1 r_3/(r_1 + r_3)\) である。
甲球の直径が 6 寸，丙球の直径が 2 寸のとき，乙球の直径は 3 寸である。

描画関数プログラムのソースを見る

function draw(r1, r3, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, z2, R) = (2*r1*r3/(r1 + r3), (-r1^2 - r1*sqrt(r1^2 + 6*r1*r3 + r3^2) + r3^2 + r3*sqrt(r1^2 + 6*r1*r3 + r3^2))/(2*(r1 + r3)), r1/2 + r3/2 + sqrt(r1^2 + 6*r1*r3 + r3^2)/2)
x2 = √2r2
@printf("r1 = %g; r3 = %g; r2 = %g; z2 = %g; R = %g\n", r1, r3, r2, z2, R)
p1 = plot(xlabel="x-axis", ylabel="z-axis")
circle(0, 0, R, :cyan)
circle(0, R - r1, r1, :blue)
circle2(x2, z2, r2)
circle(0, z2, r2)
circle(0, r3 - R, r3, :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, R - r1, "甲円:r1\n(0,0,R-r1)", :blue, :center, delta=-delta)
point(x2, z2, "乙円:r2\n(x2,0,z2)", :red, :center, delta=-delta)
point(0, r3 - R, "丙円:r3,(0,0,r3-R)", :green, :center, delta=-delta)
end
p2 = plot(xlabel="x-axis", ylabel="y-axis")
circle(0, 0, R, :cyan)
circle(0, 0, r1, :blue)
rotate(x2, 0, r2, :red, angle=90)
circle(0, 0, r3, :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)
end
plot(p1, p2)
end;