算額（その464）

算額（その464）

算額（その463）「埼玉県秩父市大宮 秩父神社 明治20年」を解いていて，条件を間違えてしまってできた図が案外きれいだったので記録しておこう。

山口正義(2015): やまぶき, 第27号
https://yamabukiwasan.sakura.ne.jp/ymbk27.pdf

大円内に甲円 2 個，乙円 4 個，丙円　4 個，丁円 2 個が入っている。甲円の直径が 1 寸のとき，乙円の直径はいかほどか。

大円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (0, r0 - r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, r3)
乙円の半径と中心座標を r4, (0, r4)
として以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms r0::positive, r1::positive, r2::positive, x2::negative,
     r3::positive, x3::negative, r4::positive;
#@syms r0, r1, r2, x2, r3, x3, r4;

eq1 = x2^2 + r2^2 - (r0 - r2)^2
eq2 = x2^2 + (r0 - r1 - r2)^2 - (r1 + r2)^2
eq3 = x3^2 + (r0 - r1 - r3)^2 - (r1 + r3)^2
eq4 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq5 = x3^2 + (r3 - r4)^2 - (r3 + r4)^2
eq6 = 2r4 + 2r1 - r0
res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (r0, r2, x2, r3, x3, r4))

   1-element Vector{NTuple{6, Sym}}:
    (4*r1*(sqrt(2) + 3)/7, r1*(5 + 4*sqrt(2))/14, -2*r1*sqrt(7*sqrt(2) + 21)/7, 2*r1*(1 + 5*sqrt(2))/49, 2*r1*sqrt(sqrt(2) + 3)*(-4*sqrt(7) + sqrt(14))/49, r1*(-1 + 2*sqrt(2))/7)

乙円の直径は甲円の直径の (5 + 4*sqrt(2))/14 倍である。

すなわち甲円の直径が 1 寸のとき，乙円の直径は 0.7612038749637414 = 7分6厘1毛あまりである。

r0 = 1.2612;  r1 = 0.5;  r2 = 0.380602;  x2 = -0.794104;  r3 = 0.164716;  x3 = -0.293341;  r4 = 0.130602
甲円の直径 = 1;  乙円の直径 = 0.761204

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1/2
   (r0, r2, x2, r3, x3, r4) = (4*r1*(sqrt(2) + 3)/7, r1*(5 + 4*sqrt(2))/14, -2*r1*sqrt(7*sqrt(2) + 21)/7, 2*r1*(1 + 5*sqrt(2))/49, 2*r1*sqrt(sqrt(2) + 3)*(-4*sqrt(7) + sqrt(14))/49, r1*(-1 + 2*sqrt(2))/7)
   @printf("r0 = %g;  r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  r4 = %g\n", r0, r1, r2, x2, r3, x3, r4)
   @printf("甲円の直径 = %g;  乙円の直径 = %g\n", 2r1, 2r2)
   plot()
   circle(0, 0, r0)
   circle(0, r0 - r1, r1, :blue)
   circle(0, r1 - r0, r1, :blue)
   circle4(x2, r2, r2, :orange)
   circle4(x3, r3, r3, :green)
   circle(0, r4, r4)
   circle(0, -r4, r4)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r0 - r1, " 甲円:r1\n (0,r0-r1)", :blue, :left, :vcenter)
       point(x2, r2, " 乙円:r2\n (x2,r2)", :orange, :left, :vcenter)
       point(x3, r3, "丙円:r3(r3+r4,0) ", :black, :right, :vcenter)
       point(0, r4, " 丁円:r4(0,r4)", :black, :left, :vcenter)
       point(r0, 0, "r0 ", :red, :right, :bottom, delta=delta/4)
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;