算額（その661）

算額（その661）

茨城県東海村 村松虚空蔵（村松虚空蔵尊） 文化15年(1818)
中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

正三角形内に全円が内接している。2本の斜線を引き，区画領域に甲円，乙円を入れる。甲円の直径が 6 寸のとき，乙円の直径を求めよ。

正三角形の一辺の長さを 2a とし，斜線と底辺の交点座標を (b, 0) とする。
全円の半径と中心座標を r0, (0, r0)
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (0, 2r0 + r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, b::positive, r0::positive,
     r1::positive, x1::positive, r2::positive
eq0 = r0 - a/√(Sym(3))
eq1 = (a - x1)/√(Sym(3)) - r1
eq2 = r2/(√(Sym(3))a - 2r0 - r2) - b/sqrt(b^2 + 3a^2)
eq3 = x1^2 + (r0 - r1)^2 - (r0 + r1)^2
eq4 = dist(0, √(Sym(3))a, b, 0, x1, r1) - r1^2
res = solve([eq0, eq1, eq2, eq3, eq4], (a, b, r0, r2, x1))

   2-element Vector{NTuple{5, Sym}}:
    (3*sqrt(3)*r1, 11*sqrt(3)*r1/7, 3*r1, 33*r1/49, 2*sqrt(3)*r1)
    (3*sqrt(3)*r1, 3*sqrt(3)*r1, 3*r1, r1, 2*sqrt(3)*r1)

2 組の解が得られるが，最初のものが適解である。

乙円の半径は甲円の半径の 33/49 倍である。「術」でも「置甲円径三十三因四十九除而得乙円径」とあり，一致する。
甲円の直径が 6 のとき，乙円の直径は 6*33/49 = 4.040816326530612 となる。「答」には，「乙円径四寸九分寸之二」とあるが，「乙円径四寸四九分寸之二（4 寸と 2/49）」の誤りだそうだ。

4 + 2/49, (4*49 + 2)/49, 198/49

   (4.040816326530612, 4.040816326530612, 4.040816326530612)

その他のパラメータは以下の通り。

 r1 = 3;  a = 15.5885;  b = 8.16538;  r0 = 9;  r2 = 2.02041;  x1 = 10.3923

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 6//2
   (a, b, r0, r2, x1) = (3*sqrt(3)*r1, 11*sqrt(3)*r1/7, 3*r1, 33*r1/49, 2*sqrt(3)*r1)
   @printf("乙円の直径 = %g;  r1 = %g;  a = %g;  b = %g;  r0 = %g;  r2 = %g;  x1 = %g\n", 2r2, r1, a, b, r0, r2, x1)
   plot([a, 0, -a, 0], [0, √3a, 0, 0], color=:green, lw=0.5)
   plot!([b, 0, -b], [0, √3a, 0], color=:orange, lw=0.5)
   circle(0, r0, r0)
   circle(x1, r1, r1, :blue)
   circle(-x1, r1, r1, :blue)
   circle(0, 2r0 + r2, r2, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(0, r0, " 全円:r0\n (0,r0)", :red, :left, :vcenter)
       point(x1, r1, " 甲円:r1,(x1,r1)", :black, :center, delta=-delta/2)
       point(0, 2r0 + r2, " 乙円:r2,(2r0+r2)", :black, :left, :vcenter)
       point(0, √3a, " √3a", :green, :left, :bottom)
       point(a, 0, "a", :green, :left, :bottom, delta=delta/2)
   end
end;