算額（その348）

算額（その348）

岐阜県大垣市赤坂町 金生山明星輪寺 元治2年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

長方形内に青の半円 2 個，赤円 1 個，緑円 2 個，浅青円 1 個，白円 2 個，および 薄紅円 2 個，無名円 2 個を入れる。緑円は青，赤，浅青のいずれにも外接し，緑，浅青円は長方形の辺にも接している。

注：算額では，緑円と赤円，浅青円と赤円が接していないような，誤解を与えかねない図になっている。

長方形の左下隅を原点とし，以下のように記号を定める。
青円の半径と中心座標 r1, (r1, 0)
赤円の半径と中心座標 r2, (2r1 - r2, y2)
緑円の半径と中心座標 r3, (r3, y3)
浅青円の半径と中心座標 r5, (r5, y2)
白円の半径と中心座標 r4, (2r1 - r4, y4)
薄紅円の半径と中心座標 r6, (x6, y6)
無名円の半径と中心座標 r7, (x7, y7)
長方形の幅は 2r1，高さは 2y2 である。

r1 を既知とすれば，残りの未知数は 13 変数なので，13 元連立方程式になるが，solve() では一度には解けないので，薄紅円と無名円を除いてまず解を求め，しかる後に得られた解を既知として，薄紅円と無名円のパラメータを求める。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, y2::positive, r3::positive, y3::positive,
     r4::positive, y4::positive, r5::positive, r6::positive, x6::positive,
     y6::positive, r7::positive, x7::positive, y7::positive;

eq1 = (r1 - r2)^2 + y2^2 - (r1 + r2)^2
eq2 = (r1 - r3)^2 + y3^2 - (r1 + r3)^2
eq3 = (r1 - r4)^2 + y4^2 - (r1 + r4)^2
eq5 = (2r1 - r2 - r3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq6 = (r2 - r4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq9 = (r3 - r5)^2 + (y2 - y3)^2 - (r3 + r5)^2
eq13 = 2r5 + 2r2 - 2r1
res1 = solve([eq1, eq2, eq3, eq5, eq6, eq9, eq13], (r2, y2, r3, y3, r4, y4, r5))
solve([eq1, eq2, eq3, eq5, eq6, eq9, eq13, eq4, eq7, eq8, eq10, eq11, eq12], (r2, y2, r3, y3, r4, y4, r5, r6, x6, y6, r7, x7, y7))

   2-element Vector{NTuple{7, Sym}}:
    (3*r1/4, sqrt(3)*r1, r1/3, 2*sqrt(3)*r1/3, 3*r1*(7 - 4*sqrt(3)), 2*r1*(-3 + 2*sqrt(3)), r1/4)
    (3*r1/4, sqrt(3)*r1, r1/3, 2*sqrt(3)*r1/3, 3*r1*(4*sqrt(3) + 7), 2*r1*(3 + 2*sqrt(3)), r1/4)

最初の組のものが適解である。
浅青円の半径は r1/4 すなわち，青円の半径の 1/4 である。

得られた r2, y2, r3, y3, r4, y4, r5 を既知として，次に，r6, x6, y6, r7, x7, y7 に関する連立方程式を解く。

using SymPy

@syms r1::positive, r2::positive, y2::positive, r3::positive, y3::positive,
     r4::positive, y4::positive, r5::positive, r6::positive, x6::positive,
     y6::positive, r7::positive, x7::positive, y7::positive;

(r2, y2, r3, y3, r4, y4, r5) = (3*r1/4, sqrt(Sym(3))*r1, r1/3, 2*sqrt(Sym(3))*r1/3, 3*r1*(7 - 4*sqrt(Sym(3))), 2*r1*(-3 + 2*sqrt(Sym(3))), r1/4)
eq4 = (x6 - r1)^2 + y6^2 - (r1 + r6)^2
eq7 = (2r1 - r2 - x6)^2 + (y2 - y6)^2 - (r2 + r6)^2
eq8 = (2r1 - r2 - x7)^2 + (y2 - y7)^2 - (r2 + r7)^2
eq10 = (x6 - r3)^2 + (y6 - y3)^2 - (r3 + r6)^2
eq11 = (x7 - r3)^2 + (y7 - y3)^2 - (r3 + r7)^2
eq12 = (x7 - r5)^2 + (y2 - y7)^2 - (r5 + r7)^2
res2 = solve([eq4, eq7, eq8, eq10, eq11, eq12], (r6, x6, y6, r7, x7, y7))

   1-element Vector{NTuple{6, Sym}}:
    (3*r1*(-8 + 5*sqrt(3))/22, r1*(32 - 9*sqrt(3))/22, 3*r1/11 + 5*sqrt(3)*r1/11, 3*r1*(-25 + 16*sqrt(3))/143, r1*(109 - 24*sqrt(3))/143, 6*r1/143 + 122*sqrt(3)*r1/143)

二段階に分けて求めたパラメータで図を描き，解が正しいことを確認する。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1
   (r2, y2, r3, y3, r4, y4, r5) = (3*r1/4, sqrt(3)*r1, r1/3, 2*sqrt(3)*r1/3, 3*r1*(7 - 4*sqrt(3)), 2*r1*(-3 + 2*sqrt(3)), r1/4)
   (r6, x6, y6, r7, x7, y7)= (3*r1*(-8 + 5*sqrt(3))/22, r1*(32 - 9*sqrt(3))/22, 3*r1/11 + 5*sqrt(3)*r1/11, 3*r1*(-25 + 16*sqrt(3))/143, r1*(109 - 24*sqrt(3))/143, 6*r1/143 + 122*sqrt(3)*r1/143)
   plot([0, 2r1, 2r1, 0, 0], [0, 0, 2y2, 2y2, 0], color=:black, lw=0.5)
   circlef(r1, 0, r1, :dodgerblue, beginangle=0, endangle=180)
   circlef(r1, 2y2, r1, :dodgerblue, beginangle=180, endangle=360)
   circlef(2r1 - r2, y2, r2, :red)
   circlef(r3, y3, r3, :green)
   circlef(r3, 2y2 - y3, r3, :green)
   circle(2r1 - r4, y4, r4, :orange)
   circle(2r1 - r4, 2y2 - y4, r4, :orange)
   circlef(r5, y2, r5, :skyblue1)
   circlef(x6, y6, r6, :pink)
   circlef(x6, 2y2 - y6, r6, :pink)
   circlef(x7, y7, r7, :orange)
   circlef(x7, 2y2 - y7, r7, :orange)
   if more
       point(r1, 0, "r1", :black, :center, :bottom)
       point(2r1 - r2, y2, "(2r1-r2,y2)", :black, :center, :bottom)
       point(r3, y3, "(r3,y3)", :black, :center, :bottom)
       point(2r1 - r4, y4, "(2r1-r4,y4)   ", :black, :center, :bottom)
       point(r5, y2, "(r5,y2)", :black, :center, :bottom)
       point(x6, y6, "  (x6,y6)", :black, :left, :bottom)
       point(x7, y7, "  (x7,y7)", :black, :left, :vcenter)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;