算額（その342）

算額（その342）

和算で遊ぼう！！ 「三春まちなか寺子屋」2017レポート
https://miharukoma.com/wp-content/uploads/2018/01/%E4%B8%89%E6%98%A5%E3%81%BE%E3%81%A1%E3%81%AA%E3%81%8B%E5%AF%BA%E5%AD%90%E5%B1%8B2017%E3%83%AC%E3%83%9D%E3%83%BC%E3%83%88.pdf
8月 直毘神社、沫なぎ神社

扇面（円弧）内に，図のように相交わる大円の一部（弧）と小円があり，その隙間に内接する全円を入れる。大円の直径が 12寸，小円の直径が 4寸，扇長（扇の半径）が 4寸5分のとき，全円の直径を求めよ。

扇面の円弧と大円の円弧の交点を (x, y) とする。
大円の半径，中心座標を r0, (0, y0)
扇の半径，中心座標を r2, (0, 0)
小円の円の半径，中心座標を r3, (0, r2 - r3)
全円の円の半径，中心座標を r4, (0, r2 - r4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms y0::positive, x::positive, y::positive, r0::positive, r2::positive, r3::positive, r4::positive

(r0, r2, r3) = (6, 9//2, 2)
eq1 = x^2 + (y - y0)^2 - r0^2
eq2 = x^2 + y^2 - r2^2
eq3 = r2 - 2r4 - y0 + r0
eq4 = r0 - 2r2 + 2r4;

res = solve([eq1, eq2, eq3, eq4], (r4, y0, x, y))

   2-element Vector{NTuple{4, Sym}}:
    (-(r0 - 2*r2)/2, 2*r0 - r2, -sqrt(3)*r0*sqrt(-(r0 - 2*r2)*(3*r0 - 2*r2))/(2*(2*r0 - r2)), (3*r0^2 - 4*r0*r2 + 2*r2^2)/(2*(2*r0 - r2)))
    (-(r0 - 2*r2)/2, 2*r0 - r2, sqrt(3)*r0*sqrt(-(r0 - 2*r2)*(3*r0 - 2*r2))/(2*(2*r0 - r2)), (3*r0^2 - 4*r0*r2 + 2*r2^2)/(2*(2*r0 - r2)))

2 番目の対が適解である。また，若干簡約化できる。
これによれば，全円の半径は「扇長 - 大円の半径/2 = 4.5 - 6/2 = 1.5」なので，直径は 3 である。
全円の直径は「2扇長 - 大円の半径 = 2*4.5 - 6 = 3」でもよい。

術では，全円の直径は「(扇長＋大円の直径) - (扇長×大円の直径)/小円の直径」としている。「小円の直径 = 2*扇長*(4/9)」なので
全円の直径は簡約化され，「扇長 - 大円の直径/8 = 4.5 - 12/8 = 3」となる。

   r4 = 1.5;  y0 = 7.5;  x = 3.6;  y = 2.7

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r2, r3) = (6, 9//2, 2)
   r4 = -r0/2 + r2
   y0 = 2*r0 - r2
   x = r0*sqrt(-9*r0^2 + 24*r0*r2 - 12*r2^2)/(2*(2*r0 - r2))
   y = (3*r0^2/2 - 2*r0*r2 + r2^2)/(2*r0 - r2)
   @printf("r4 = %g;  y0 = %g;  x = %g;  y = %g\n", r4, y0, x, y)
   θ = round(Int64, atand(y, x))
   θ2 = round(Int64, atand(y0 - y, x))
   plot()
   circle(0, y0, r0, beginangle=180+θ2, endangle=360-θ2)
   circle(0, 0, r2, :black, beginangle=θ, endangle=180-θ)
   segment(0, 0, x, y)
   segment(0, 0, -x, y)
   circle(0, r2 - r4, r4, :green)
   circle(0, r2 - r3, r3, :blue)
   circle(0, 0, r2 - 2r4, :black, beginangle=θ, endangle=180-θ)
   if more
       point(x, y, "(x,y)")
       point(0, y0, " y0", :red)
       point(0, r2, " r2", :red, :left, :bottom)
       point(0, r0-r2, " r0-r2", :red, :left, :bottom)
       point(0, r2 - r4, " r2-r4", :green)
       point(0, r2 - 2r4, " r2-2r4", :green)
       point(0, r2 - r3, " r2-r3", :blue)
       point(0, r2-2r3, " r2-2r3", :blue)
       #point(0, y0 - r0, " y0-r0", :blue)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;