算額（その1099）

算額（その1099）

六十六 岩手県花泉町金沢字大柳 金沢八幡宮 明治29年(1896)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円3個，外円，直角三角形2個

全円の中に交差する（合同な）直角三角形を 2 個と，大円，小円を 1 個ずつ容れる。大円と小円の直径が与えられたとき，全円の直径を求めよ。

内接する直角三角形なので，斜辺は全円の直径である。一つの直角三角形の斜辺が x 軸上になるように図を描く。
全円の半径と中心座標を　R, (0, 0)
大円の直径と中心座標を r1, (x1, r1); x1 < 0
小円の半径と中心座標を r2, (0, R - r2)
直角三角形の直径上ではない直角の頂点の座標を (x0, y0)
とおき，以下の方程式を解く。

なお，証明はしないが x1 = r2 - r1 である（図の水色の垂直線を参照）ので，方程式の解は簡単に求まる。

include("julia-source.txt")

using SymPy
@syms R::positive, x0::positive, y0::positive,
     r1::positive, x1::negative, r2::positive
y0 = sqrt(R^2 - x0^2)
x1 = r2 - r1  # 未証明
eq1 = x1^2 + r1^2 - (R - r1)^2
eq2 = dist2(-x0, y0, x0, -y0, -x1, r1, r1)
eq3 = dist2(-R, 0, x0, y0, 0, R - r2, r2);

最初に eq1 を解くことで全円の半径 R が求まる。「術」と同じ式になる。

ans_R = solve(eq1, R)[2]
ans_R |> println

   r1 + sqrt(2*r1^2 - 2*r1*r2 + r2^2)

次いで，eq2 を解いて x0 を求める。
なお，eq3 を解いても求まるが，本質的におなじ式であるが eq3 から求める式は SymPy では簡約化できず，長く複雑なものになる。

ans_x0 = solve(eq2, x0)[3] |> factor
ans_x0 |> println

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

大円，小円の直径が 20, 10 のとき，全円の直径は 42.3607 である。

   r1 = 10;  r2 = 5;  R = 21.1803;  x0 = 12.7082;  y0 = 16.9443;  x1 = -5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (10, 5)
   R = r1 + sqrt(2*r1^2 - 2*r1*r2 + r2^2)
   x0 = R*r2*(2*r1 - r2)/(2*r1^2 - 2*r1*r2 + r2^2)
   y0 = sqrt(R^2 - x0^2)
   x1 = r2 - r1  # 未証明
   @printf("大円，小円の直径が %g, %g のとき，全円の直径は %g である。\n", 2r1, 2r2, 2R)
   @printf("r1 = %g;  r2 = %g;  R = %g;  x0 = %g;  y0 = %g;  x1 = %g\n", r1, r2, R, x0, y0, x1)
   plot()
   circle(0, 0, R, :blue)
   circle(x1, -r1, r1)
   circle(0, R - r2, r2, :green)
   plot!([x0, -x0, R, x0], [-y0, y0, 0, -y0], color=:magenta, lw=0.5)
   plot!([R, -R, x0, R], [0, 0, y0, 0], color=:magenta, lw=0.5)
   vline!([r2], lw=0.5)
   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(x1, -r1, "大円:r1,(x1,-r1)", :red, :center, delta=-delta/2)
       point(0, R - r2, "小円:r2,(0,R-r2)", :green, :center, delta=-delta/2)
       point(r2, 0, "r2", :green, :center, delta=-delta/2)
       point(x0, y0, "(x0,y0)  ", :magenta, :right, :bottom, delta=-delta/2)
       point(x0, -y0, "(x0,-y0)  ", :magenta, :right, :bottom, delta=-delta/2)
       point(-x0, y0, "  (-x0,y0)", :magenta, :left, :bottom, delta=-delta/2)
   end
end;