算額（その517）

算額（その517）

千葉胤秀: 『算法新書』巻の五，文政13年(1830)
一関市博物館>>和算に挑戦>>平成23年度出題問題&解答例>>平成23年度出題問題(3)[上級問題]
https://www.city.ichinoseki.iwate.jp/museum/wasan/h23/hard.html

外円内に5個の円が内接・外接している。外円の直径は 6 寸，甲円の直径が 2 寸のとき，丙円の直径を求めよ。

外円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (0, r0 - r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (0, r3 - 0)
丁円の半径と中心座標を r4, (0, r0 - 2r1 - r4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0, r1, r2, x2, y2, r3, r4

eq1 = x2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + (y2 - r3 + r0)^2 - (r2 + r3)^2
eq3 = x2^2 + y2^2 - (r0 - r2)^2
eq4 = x2^2 + (r0 - 2r1 - r4 - y2)^2 - (r2 + r4)^2
eq5 = r1 + r3 + r4 - r0

res = solve([eq1, eq2, eq3, eq4, eq5], (r2, x2, y2, r3, r4))

   3-element Vector{NTuple{5, Sym}}:
    (0, 0, r0, r0, -r1)
    (2*r0*r1*(r0 - r1)/(r0^2 + r1^2), -2*sqrt(2)*r0*r1*(r0 - r1)/(r0^2 + r1^2), r0*(r0^2 - 2*r0*r1 - r1^2)/(r0^2 + r1^2), r0*(r0 - r1)/(r0 + r1), (r0*r1 - r1^2)/(r0 + r1))
    (2*r0*r1*(r0 - r1)/(r0^2 + r1^2), 2*sqrt(2)*r0*r1*(r0 - r1)/(r0^2 + r1^2), r0*(r0^2 - 2*r0*r1 - r1^2)/(r0^2 + r1^2), r0*(r0 - r1)/(r0 + r1), (r0*r1 - r1^2)/(r0 + r1))

3 組の解が得られるが，3 番目の解が適解である。

乙円の半径(r2)は r0*(r0 - r1)/(r0 + r1) で，外円の直径(2r0)が 6 寸，甲円の直径(2r1)が 2 寸のとき，丙円の直径は 3.0 である。

r0 = 3
r1 = 1
2(r0*(r0 - r1)/(r0 + r1))

   3.0

   r0 = 3;  r1 = 1;  r2 = 1.2;  x2 = 1.69706;  y2 = 0.6;  r3 = 1.5;  r4 = 0.5

「甲円の直径（半径）の逆数」と「丙円の直径（半径）の逆数」の和は，「丙円の逆数の2倍」と等しい。
今回の図では左右にある丙円は同じ大きさであるが大きさが異なる場合も「丙円1の逆数と丙円2の逆数の和」に等しい。

1/r1 + 1/res[3][4] |> simplify |> println

   (1.0*r0*(r0 - r1) + r0 + r1)/(r0*(r0 - r1))

2/res[3][1] |> simplify |>  println

   (r0^2 + r1^2)/(r0*r1*(r0 - r1))

1/r1 + 1/r3 = 1.66667;  2/r2 = 1.66667

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1) = (6, 2) .// 2
   (r2, x2, y2, r3, r4) = (
       2*r0*r1*(r0 - r1)/(r0^2 + r1^2),
       2*sqrt(2)*r0*r1*(r0 - r1)/(r0^2 + r1^2),
       r0*(r0^2 - 2*r0*r1 - r1^2)/(r0^2 + r1^2),
       r0*(r0 - r1)/(r0 + r1),
       (r0*r1 - r1^2)/(r0 + r1))
   @printf("r0 = %g;  r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  r3 = %g;  r4 = %g\n", r0, r1, r2, x2, y2, r3, r4)
   @printf("1/r1 + 1/r3 = %g;  2/r2 = %g\n", 1/r1 + 1/r3, 2/r2)
   plot()
   circle(0, 0, r0)
   circle(0, r0 - r1, r1, :blue)
   circle(x2, y2, r2, :green)
   circle(-x2, y2, r2, :green)
   circle(0, r3 - r0, r3, :magenta)
   circle(0, r0 - 2r1 - r4, r4, :orange)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r0, " r0", :red, :left, :bottom, delta=delta/2)
       point(0, r0 - r1, " 甲円:r1\n (0,r0-r1)", :blue, :left, :vcenter)
       point(x2, y2, " 乙円:r2,(x2,y2)", :green, :center, :top, delta=-delta/2)
       point(0, r3 - r0, " 丙円:r3\n (0,r3-r0)", :magenta, :left, :vcenter)
       point(0, r0 - 2r1 - r4, " 丁円:r4\n (0,r0-2r1-r4)", :black, :left, :bottom, delta=delta)
   end
end;