算額（その470）

算額（その470）

宮城県角田市横倉 愛宕神社 明治15年(1882)1月
http://www.wasan.jp/miyagi/yokokuraatago.html

徳竹亜紀子，谷垣美保: 2021年度の算額調査，仙台高等専門学校名取キャンパス 研究紀要，第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf

甲円内に弦を挟んで乙円 1 個，丙円 1 個，丁円 2 個が入っている。
乙円，丙円，丁円の直径がそれぞれ 4 寸，3 寸，1 寸のとき，甲円の直径を求めよ。

甲円の半径と中心座標を r1, (0, 0)
乙円の半径と中心座標を r2, (x2, 2r4 - r1 + r2)
丙円の半径と中心座標を r3, (x3, 2r4 - r1 + r3)
丁円の半径と中心座標を r4, (0, r4 - r1), (x4, 3r4 - r1)
とし，以下のように記号を定め方程式を解く。

include("julia-source.txt")

using SymPy
@syms r1::positive, r2::positive, x2::positive,
     r3::positive, x3::negative,
     r4::positive, x4::negative;
eq1 = x2^2 + (2r4 - r1 + r2)^2 - (r1 - r2)^2
eq2 = x3^2 + (2r4 - r1 + r3)^2 - (r1 - r3)^2
eq3 = (x2 - x4)^2 + (r2 - r4)^2 - (r2 + r4)^2
eq4 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2
solve([eq1, eq2, eq3, eq4], (r1, x2, x3, x4))

   2-element Vector{NTuple{4, Sym}}:
    (r2 + r3 + r4, 2*sqrt(r4)*(-sqrt(r2)*sqrt(r3) + r3)*sqrt(2*sqrt(r2)*sqrt(r3) + r2 + r3)/(r2 - r3), 2*sqrt(r4)*(sqrt(r2)*sqrt(r3) - r2)*sqrt(2*sqrt(r2)*sqrt(r3) + r2 + r3)/(r2 - r3), -2*sqrt(r4)*sqrt(2*sqrt(r2)*sqrt(r3) + r2 + r3))
    (r2 + r3 + r4, 2*sqrt(r4)*(sqrt(r2)*sqrt(r3) + r3)*sqrt(-2*sqrt(r2)*sqrt(r3) + r2 + r3)/(r2 - r3), -2*sqrt(r4)*(sqrt(r2)*sqrt(r3) + r2)*sqrt(-2*sqrt(r2)*sqrt(r3) + r2 + r3)/(r2 - r3), -2*sqrt(r4)*sqrt(-2*sqrt(r2)*sqrt(r3) + r2 + r3))

2 組の解が得られるが，2 組目のものが適解である。

   r1 = 4;  x2 = 1.73205;  x3 = -2;  x4= -0.267949
   甲円の直径 = 8

甲円の半径は r2 + r3 + r4 すなわち，乙円，丙円，丁円の半径の和である。
乙円，丙円，丁円の直径がそれぞれ 4 寸，3 寸，1 寸ならば，甲円の直径は 4 + 3 + 1 = 8 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, r4) = (4, 3, 1) .// 2
   (r1, x2, x3, x4) = (r2 + r3 + r4, 2*sqrt(r4)*(sqrt(r2)*sqrt(r3) + r3)*sqrt(-2*sqrt(r2)*sqrt(r3) + r2 + r3)/(r2 - r3), -2*sqrt(r4)*(sqrt(r2)*sqrt(r3) + r2)*sqrt(-2*sqrt(r2)*sqrt(r3) + r2 + r3)/(r2 - r3), -2*sqrt(r4)*sqrt(-2*sqrt(r2)*sqrt(r3) + r2 + r3))
   @printf("r1 = %g;  x2 = %g;  x3 = %g;  x4= %g\n", r1, x2, x3, x4)
   @printf("甲円の直径 = %g\n", 2r1)
   plot()
   circle(0, 0, r1, :red)
   circle(x2, 2r4 - r1 + r2, r2, :magenta)
   circle(x3, 2r4 - r1 + r3, r3, :blue)
   circle(x4, 3r4 - r1, r4, :green)
   circle(0, r4 - r1, r4, :green)
   segment(-sqrt(r1^2 - (2r4 - r1)^2), 2r4 - r1, sqrt(r1^2 - (2r4 - r1)^2), 2r4 - r1)
   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, r1, " r1", :red, :left, :bottom, delta=delta/2)
       point(x2, 2r4 - r1 + r2, "乙円:r2,(x2,2r4-r1+r2)", :magenta, :center, :top, delta=-delta)
       point(x3, 2r4 - r1 + r3, "丙円:r3,(x3,2r4-r1+r3)", :blue, :center, :top, delta=-delta)
       point(0, r4 - r1, " 丁円:r4\n (0,r4-r1)", :black, :left, :vcenter)
       point(x4, 3r4 - r1, "丁円:r4,(x4,3r4-r1) ", :black, :right, :vcenter)
       point(0.2r1, 0.5r1, " 甲円:r1,(0,0)", :red, :left, :vcenter)
       point(0, 2r4 - r1, " 2r4-r1", :black, :left, :bottom, delta=delta)
   end
end;