算額（その524）

算額（その524）

和算図形問題あれこれ - 令和4年10月の問題-No.2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html

外円の中に甲円と乙円が入っている。甲円と乙円の直径が与えられたとき，外円の直径を求めよ。

外円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (0, r0 - r1)
乙円の半径と中心座標を r2, (x2, r3 - r0 - r2)
外円の下部にある円弧の半径と中心座標を r3, (0, -r0)
として，以下の連立方程式を解く。
乙円の半径を求めるためには eq1, eq2, eq3 を解けばよい。
eq4, eq5 は外円と円弧の交点座標を求めるためのものである。

include("julia-source.txt");

using SymPy
@syms r0::positive, r1::positive, r2::positive, x2::positive,
   r3::positive, xa::positive, ya::negative;

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

   4-element Vector{NTuple{5, Sym}}:
    (r1*(r1 - r2)/(r1 - 2*r2), -2*sqrt(2)*sqrt(r1)*r2*sqrt(1/(r1 - 2*r2)), 2*r1*r2/(r1 - 2*r2), -2*r1^(3/2)*r2*sqrt(1/(r1 - 2*r2))/(r1 - r2), -r1*(r1^2 - 2*r1*r2 - r2^2)/((r1 - 2*r2)*(r1 - r2)))
    (r1*(r1 - r2)/(r1 - 2*r2), -2*sqrt(2)*sqrt(r1)*r2*sqrt(1/(r1 - 2*r2)), 2*r1*r2/(r1 - 2*r2), 2*r1^(3/2)*r2*sqrt(1/(r1 - 2*r2))/(r1 - r2), -r1*(r1^2 - 2*r1*r2 - r2^2)/((r1 - 2*r2)*(r1 - r2)))
    (r1*(r1 - r2)/(r1 - 2*r2), 2*sqrt(2)*sqrt(r1)*r2*sqrt(1/(r1 - 2*r2)), 2*r1*r2/(r1 - 2*r2), -2*r1^(3/2)*r2*sqrt(1/(r1 - 2*r2))/(r1 - r2), -r1*(r1^2 - 2*r1*r2 - r2^2)/((r1 - 2*r2)*(r1 - r2)))
    (r1*(r1 - r2)/(r1 - 2*r2), 2*sqrt(2)*sqrt(r1)*r2*sqrt(1/(r1 - 2*r2)), 2*r1*r2/(r1 - 2*r2), 2*r1^(3/2)*r2*sqrt(1/(r1 - 2*r2))/(r1 - r2), -r1*(r1^2 - 2*r1*r2 - r2^2)/((r1 - 2*r2)*(r1 - r2)))

4 組の解が得られる。符号を無視すればどれでもよいが，図のような位置関係ならば 4 番目の解が適解である。

外円の半径は r1*(r1 - r2)/(r1 - 2r2) である。

甲円，乙円の直径が7, 2 の場合，外円の直径は 35/3 ≒ 11.6667 である。

   r0 = 5.83333;  r1 = 3.5;  r2 = 1;  x2 = 4.32049;  r3 = 4.66667;  xa = 4.27707;  ya = -3.96667
   外円の直径 = 11.6667

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (7, 2) .// 2
   (r0, x2, r3, xa, ya) = (r1*(r1 - r2)/(r1 - 2*r2), 2*sqrt(2)*sqrt(r1)*r2*sqrt(1/(r1 - 2*r2)), 2*r1*r2/(r1 - 2*r2), 2*r1^(3/2)*r2*sqrt(1/(r1 - 2*r2))/(r1 - r2), -r1*(r1^2 - 2*r1*r2 - r2^2)/((r1 - 2*r2)*(r1 - r2)))
   @printf("r0 = %g;  r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  xa = %g;  ya = %g\n", r0, r1, r2, x2, r3, xa, ya)
   @printf("外円の直径 = %g\n", 2r0)
   plot()
   circle(0, 0, r0)
   circle(0, r0 - r1, r1, :blue)
   circle(x2, r3 - r0 - r2, r2, :green)
   circle(-x2, r3 - r0 - r2, r2, :green)
   x = sqrt(r0^2 - (r3 - r0)^2)
   segment(-x, r3 - r0, x, r3 - r0, :orange)
   θ = atand((ya + r0)/xa)
   circle(0, -r0, r3, :magenta, beginangle=θ, endangle=180 - θ)    
   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, r3 - r0, " r3-r0", :black, :left, :bottom, delta=delta)
       point(0, r0 - r1, " 甲円:r1\n (0,r0-r1)", :blue, :left, :vcenter)
       point(x2, r3 - r0 - r2, "乙円:r2(x2,r3-r0-r2) ", :black, :right, :vcenter)
       point(xa, ya, " (xa,ya)", :red, :left, :vcenter)
   end
end;