算額（その1080）

算額（その1080）

九十八 岩手県江刺市 雨宝堂 現中善観音堂 文政10年(1827)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html
キーワード：円5個，半円4個，長方形，斜線2本

長方形の中に 2 本の斜線と甲円（半円），乙円，丙円，丁円を容れる。丁円の直径が 1 寸のとき，乙円の直径はいかほどか。

原点を，中央の乙円の中心に置く。」
長方形の短辺，長辺の長さを 2a, 2r2 + 2r1
甲円の半径と中心座標を r1, (r2 + r1)
乙円の半径と中心座標を r2, (a, 0)
丙円の半径と中心座標を r3,（0, r2 + r3)
丁円の半径と中心座標を r4, (x4, y4)
とおき，以下の連立方程式を解く。

eq1, eq2 は他の条件と重複するので除外する。

残念ながら，SymPy で一度に解くことができないので，今回は既知の変数を前もって定義するという対応を取る。

右下がりの斜線と長方形の右辺を下に延長し交点を α とする。右の甲円，乙円の中心 の y 座標を β，γとすると，αβ，αγ を一辺とする直角三角形は相似で，相似比は 1:2 である。したがって，乙円と甲円の半径 r2, r1 は 1:2 である。

また，右下がりの斜線と y 軸の交点の y 座標を δ，中央下の丙円，乙円の中心の y 座標を ε，ζ とすると,δε，δζ を一辺とする直角三角形は相似で，相似比は 1:2 である。したがって，丙円と乙円の半径 r3, r2 は 1:2 である。
これを合わせると，r3:r2:r1 = 4:2:1 なので，r2 = 2r3, r1 = 4r3 という条件を加え，eq5, eq6 は除外する。

これで，もともとは 6 元連立方程式だったものが 4 元連立方程式になった。

include("julia-source.txt")

using SymPy

@syms a::positive, r1::positive, r2::positive,
     r3::positive, x3::positive, y3::positive,
     r4::positive, x4::positive, y4::positive
r2 = 2r3
r1 = 4r3
#eq1 = dist2(0, r1 + r2, a/2, 0, a, 0, r2)
#eq2 = dist2(0, r1 + r2, a/2, 0, 0, 0, r2)
eq3 = dist2(0, r1 + r2, a/2, 0, x4, y4, r4)
eq4 = dist2(0, r1 + r2, a/2, 0, 0, r2 + r3, r3)
#eq5 = dist2(0, r1 + r2, a/2, 0, a, r1 + r2, r1)
eq6 = (a - x4)^2 + (r1 + r2 - y4)^2 - (r1 + r4)^2 |> expand
eq7 = (a - x4)^2 + y4^2 - (r2 + r4)^2 |> expand;
#eq8 = r3/(r1 - r3) - r2/(r1 + r2)
#res = solve([eq3, eq4, eq5, eq6, eq7, eq8], (a, r1, r2, r3, x4, y4))

println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
#println(eq5, ",  # eq5")
println(eq6, ",  # eq6")
println(eq7, ",  # eq7")
#println(eq8, ",  # eq8")

   36*a^2*r3^2 - 12*a^2*r3*y4 - a^2*r4^2 + a^2*y4^2 - 144*a*r3^2*x4 + 24*a*r3*x4*y4 - 144*r3^2*r4^2 + 144*r3^2*x4^2,  # eq3
   8*r3^2*(a^2 - 18*r3^2),  # eq4
   a^2 - 2*a*x4 + 20*r3^2 - 8*r3*r4 - 12*r3*y4 - r4^2 + x4^2 + y4^2,  # eq6
   a^2 - 2*a*x4 - 4*r3^2 - 4*r3*r4 - r4^2 + x4^2 + y4^2,  # eq7

最も式が簡単な eq4 を解いて a を求め，残りの 3 個の方程式に代入することで，更に次元を 1 つ減らすことができる。

ans_a = solve(eq4, a)[1]
ans_a |> println

   3*sqrt(2)*r3

eq3 = eq3(a => ans_a) |> simplify |> (x -> x/18r3^2)
eq6 = eq6(a => ans_a)
eq7 = eq7(a => ans_a) |> simplify
eq3 |> println
eq6 |> println
eq7 |> println

   36*r3^2 - 24*sqrt(2)*r3*x4 - 12*r3*y4 - 9*r4^2 + 8*x4^2 + 4*sqrt(2)*x4*y4 + y4^2
   38*r3^2 - 8*r3*r4 - 6*sqrt(2)*r3*x4 - 12*r3*y4 - r4^2 + x4^2 + y4^2
   14*r3^2 - 4*r3*r4 - 6*sqrt(2)*r3*x4 - r4^2 + x4^2 + y4^2

これで，SymPy でも連立方程式を解くことができるようになる。r3, x4, y4 が，与えられる r4 を含む式で表される。

res = solve([eq3, eq6, eq7], (r3, x4, y4))[1];
res[1] |> simplify |> println
res[2] |> simplify |> println
res[3] |> simplify |> println

乙円の半径 r2 は 丙円の半径 r3 の 2倍，更に r3 は 丁円の半径 r4 の (2√2 + 3)/4 倍である。つまり，r2 = r2*(2√2 + 3)/2 である。
丁円の直径が 1 寸のとき丙円の直径は 2.914213562373095 寸である。

その他のパラメータは以下のとおりである。

r4 = 0.5;  r3 = 0.728553;  x4 = 1.61959;  y4 = 1.29044;  a = 3.09099;  r2 = 1.45711;  r1 = 2.91421

パラメータは以下のようにして逐次求めることができる。

r4 = 0.5
r3 = r4*(2*sqrt(2) + 3)/4
x4 = r4*(12 + 19*sqrt(2))/12
y4 = r4*(7 + 6*sqrt(2))/6
a = 3*sqrt(2)*r3
r2 = 2r3
r1 = 4r3
(a, r1, r2, r3, x4, y4)

   (3.0909902576697323, 2.914213562373095, 1.4571067811865475, 0.7285533905932737, 1.6195857368787003, 1.290440114519881)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r4 = 0.5
   #(a, r1, r2, r3, x4, y4) = res[1]
   r3 = r4*(2*sqrt(2) + 3)/4
   x4 = r4*(12 + 19*sqrt(2))/12
   y4 = r4*(7 + 6*sqrt(2))/6
   a = 3*sqrt(2)*r3
   r2 = 2r3
   r1 = 4r3
   @printf("丁円の直径が %g のとき，乙円の直径は %g である。\n", 2r4, 2r2)
   @printf("r4 = %g;  r3 = %g;  x4 = %g;  y4 = %g;  a = %g;  r2 = %g;  r1 = %g\n",
       r4, r3, x4, y4, a, r2, r1)
   plot([a, a, -a, -a, a], [-r2, r2 + 2r1, r2 + 2r1, -r2, -r2], color=:green, lw=0.5)
   circle(0, 0, r2, :blue)
   circle(a, 0, r2, :blue, beginangle=90, endangle=270)
   circle(-a, 0, r2, :blue, beginangle=-90, endangle=90)
   circle(0, r2 + r3, r3, :magenta)
   circle(0, 2r1 + r2 - r3, r3, :magenta)
   circle(a, r2 + r1, r1, beginangle=90, endangle=270)
   circle(-a, r2 + r1, r1, beginangle=-90, endangle=90)
   circle2(x4, y4, r4, :green)
   x01 = a*(r1 + 2r2)/2(r1 + r2)
   x02 = a*r1/2(r1 + r2)
   segment(a, -r2 - r1, -x02, r2 + 2r1, :orange)
   segment(a, -r2, a, -r2 - r1, :orange)
   segment(-x01, -r2, x02, r2 + 2r1, :orange)
   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(a, -r2 - r1, " α", :gray50, :left, :vcenter)
       point(a, 0, " β", :gray40, :left, :vcenter)
       point(a, r2 + r1, " γ", :gray40, :left, :vcenter)
       point(0, 0, "ζ ", :gray40, :right, :vcenter, delta=delta)
       point(0, r2 + r3, "ε ", :gray40, :right, :vcenter)
       point(0, r2 + r1, "δ ", :gray40, :right, :vcenter)
       point(a, r2 + r1, "甲円:r1\n(a,r2+r1) ", :red, :right, :vcenter)
       point(a, 0, "乙円:r2 \n(a,0) ", :blue, :right, :vcenter)
       point(0, 0, "乙円:r2,(0,0) ", :blue, :center, delta=-delta/2)
       point(0, r2 + r3, " 丙円:r3,(0,r2+r3) ", :black, :left, :vcenter)
       point(0, 2r1 + r2 - r3, " 丙円:r3,(0,2r1+r2-r3) ", :black, :left, :vcenter)
       point(x4, y4, "丁円:r4,(x4,y4)", :black, :left, delta=-delta/2, deltax=-3delta)
   end
end;