算額（その905）

算額（その905）

一〇〇 桶川市小針領家 氷川諏訪神社 明治30年(1897)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

球の中に「蕎麦の実」が入っている。球の直径が 3 寸のとき，蕎麦の実の一辺の長さはいかほどか。

注：「蕎麦の実」というのは，正四面体のことである。

球の半径と中心座標を r, (0, 0, 0) とする。
正四面体の一辺の長さを a とする。

頂点を A, B, C, D として 3 次元座標を割り当てる。
A: (x1, 0, z1)
B: (x2, y2, z1)
C: (x2, -y2, z1)
D: (0, 0, z4)
r = z4 および a = 2y2 である。
4 点は球の表面上にあり，どの2点間の距離も等しく a である。
以下の連立方程式を解けば，4点が決まる（a も決まる）。

include("julia-source.txt");

using SymPy
@syms a::positive, r::positive,
     x1::positive, z1::negative, x2::negative,
     y2::positive, x3::positive, z4::positive

a = 2y2
eq1 = x1^2 + z1^2 - r^2
eq2 = x2^2 + y2^2 + z1^2 - r^2
eq3 = (x1 - x2)^2 + y2^2 - a^2
eq4 = x1^2 + (r - z1)^2 - a^2
eq5 = x2^2 + y2^2 + (r - z1)^2 - a^2
solve([eq1, eq2, eq3, eq4], (x1, z1, x2, y2))

   1-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (2*sqrt(2)*r/3, -r/3, -sqrt(2)*r/3, sqrt(6)*r/3)

正四面体の一辺の長さは 2r*√6/3 = 球の直径 × √6/3 である。
球の直径が 3 寸のとき，正四面体の一辺の長さは √6 = 2.449489742783178 である。

r = 3/2
x1 = r*2√2/3
z1 = -r/3
x2 = -r*√2/3
y2 = r*√6/3
a = 2y2  # = 2r*√6/3

   2.449489742783178

術は「2 を 3 で割り，その平方根を求め，球の直径を掛ける」 sqrt(2/3) * 2r = √6/3 * 2r である。

算額（その904）

算額（その904）

一〇〇 桶川市小針領家 氷川諏訪神社 明治30年(1897)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

外円内に甲乙丙の 3 円と斜線（弦）が入っている。外円，甲円の直径がそれぞれ 3 寸，2 寸のとき，乙円の直径はいかほどか。

注：外円，甲円，丙円は一直線上にあり，弦は外円と甲円の接点を通る。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (0, r3 - R)
弦と外円の交点座標を (x, y)
とおき，以下の連立方程式を解く。
SymPy の性能上，一度に解くことができないので，逐次的に解いてゆく。

乙円の半径は eq1 だけで求めることができる。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive,
     r2::positive, x2::positive, y2::positive,
     r3::positive, x::positive, y::positive
r3 = R - r1
eq1 = r3/(2R - r3) - (R - 2r2)/R
ans_r2 = solve(eq1, r2)[1]
ans_r2 |> println
ans_r2(R => 3//2, r1 => 2/2)*2 |> println

   R*r1/(R + r1)
   1.20000000000000

乙円の半径は 外円，甲円の直径に基づいて計算できる。術の通り，外円と甲円の直径の積を外円と甲円の直径の和で割ればよい。

外円，甲円の直径がそれぞれ 3 寸，2 寸のとき，乙円の直径は 1 寸 2 分である。

算額の問題を解くということだけならば，ここまでで終了。

図を描くために，その他のパラメータを求める。
次に，r2 が決定されたという条件下で，弦の交点座標を求める。

r2 = R*r1/(R + r1)
eq2 = x^2 + y^2 - R^2 |> expand
eq3 = (R + y)*y2 - x2*x |> expand;
res2 = solve([eq2, eq3], (x, y))[1]

   (y2*(R + (R*x2^2 - R*y2^2)/(x2^2 + y2^2))/x2, (R*x2^2 - R*y2^2)/(x2^2 + y2^2))

最後に，r2, x, y が決定されたという条件下で，乙円の中心座標 x2, y2 を求める。

eq4 = x2^2 + y2^2 - (R - r2)^2 |> expand
eq5 = dist2(0, R, x, -y, x2, y2, r2) |> expand;

eq14 = eq4(x => res2[1], y => res2[2])
ans_x2 = solve(eq14, x2)[1]
ans_x2 |> println

   sqrt(R^4 - R^2*y2^2 - 2*R*r1*y2^2 - r1^2*y2^2)/(R + r1)

eq15 = eq5(x => res2[1], y => res2[2]) |> simplify |> numerator
eq25 = eq15(x2 => ans_x2);
ans_y2 = solve(eq25, y2)[2]
ans_y2 |> println

   R^2*(R - r1)/(R^2 + 2*R*r1 + r1^2)

y2 の値も，外円の半径と甲円の半径から求めることができる。
その後は，逆順にたどると全てのパラメータの値が求まる。

(R, r1) = (1.5, 1)
r2 = R*r1/(R + r1)
y2 = R^2*(R - r1)/(R^2 + 2*R*r1 + r1^2)
x2 = sqrt(R^4 - R^2*y2^2 - 2*R*r1*y2^2 - r1^2*y2^2)/(R + r1)
x = y2*(R + (R*x2^2 - R*y2^2)/(x2^2 + y2^2))/x2
y = (R*x2^2 - R*y2^2)/(x2^2 + y2^2)
(x2, y2, x, y) |> println

   (0.8818163074019442, 0.18, 0.5878775382679626, 1.3800000000000001)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1) = (1.5, 1)
   r2 = R*r1/(R + r1)
   y2 = R^2*(R - r1)/(R^2 + 2*R*r1 + r1^2)
   x2 = sqrt(R^4 - R^2*y2^2 - 2*R*r1*y2^2 - r1^2*y2^2)/(R + r1)
   x = y2*(R + (R*x2^2 - R*y2^2)/(x2^2 + y2^2))/x2
   y = (R*x2^2 - R*y2^2)/(x2^2 + y2^2)
   r3 = R - r1
   @printf("r2 = %g;  x2 = %g;  y2 = %g;  r3 = %g;  x = %g;  y = %g\n", r2, x2, y2, r3, x, y)
   plot()
   circle(0, 0, R)
   circle(0, R - r1, r1, :blue)
   circle(0, r3 - R, r3, :green)
   circle(x2, y2, r2, :magenta)
   segment(0, R, x, -y, :orange)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=1)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R, " R", :red, :left, :bottom, delta=delta/2)
       point(x, -y, "(x,y)", :red, :left, delta=-delta/2)
       point(0, R - r1, "甲円:r1 \n(0,R-r1) ", :blue, :right, :vcenter)
       point(x2, y2, "乙円:r2 \n(x2,y2) ", :magenta, :right, :vcenter)
       point(0, r3 - R, "丙円:r3,(0,r3-R) ", :green, :center, delta=-delta/2)
       point(0, r3 - R)
   end
end;