算額（その818）

算額（その818）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

長方形と，長方形の長辺を底辺として共有する二等辺三角形の中に，甲円，乙円，丙円を入れる。
丙円の直径が 1 寸のとき，乙円の直径はいかほどか。

甲円の半径と中心座標を r1, (r1, r1)
乙円の半径と中心座標を r2, (0, 2r1 + r2)
丙円の半径と中心座標を r3, (0, r3), (x3, y3)
長方形の長辺と短辺を 2r1, r1
二等辺三角形の高さを h
とおき，以下の連立方程式を解く。

乙円についての方程式 eq5 は，他の方程式とは独立なので，まずは eq1〜eq4 の連立方程式を解いて
h, r1, x3, y3 を求める。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms h::positive, r1::positive, r2::positive,
     r3::positive, x3::positive, y3::positive, d
r3 = 1//2
eq1 = r1^2 + (r1 - r3)^2 - (r1 + r3)^2
eq2 = (x3 - r1)^2 + (y3 - r1)^2 - (r1 - r3)^2
eq3 = (x3 - r1)*2r1 - h*(y3 - r1)
eq4 = dist(0, h, 2r1, 0, x3, y3) - r3^2
eq4 = numerator(apart(eq4, d));
(h, r1, x3, y3) = solve([eq1, eq2, eq3, eq4], (h, r1, x3, y3))[2]

   (2*sqrt(2*sqrt(7) + 8)/3 + 2*sqrt(7)*sqrt(2*sqrt(7) + 8)/3, 2, 3*sqrt(2*sqrt(7) + 8)/8 + 2, 3*sqrt(7)/8 + 13/8)

h, r1, x3, y3 が求められた前提で，次の方程式を解く。

eq5 = r2/(h - 2r1 - r2) - 2r1/sqrt(h^2 + 4r1^2)
res2 = solve(eq5, r2)[1]
res2 |> println

   -144/(36 + 3*sqrt(144 + (2*sqrt(2*sqrt(7) + 8) + 2*sqrt(7)*sqrt(2*sqrt(7) + 8))^2)) + 24*sqrt(2*sqrt(7) + 8)/(36 + 3*sqrt(144 + (2*sqrt(2*sqrt(7) + 8) + 2*sqrt(7)*sqrt(2*sqrt(7) + 8))^2)) + 24*sqrt(7)*sqrt(2*sqrt(7) + 8)/(36 + 3*sqrt(144 + (2*sqrt(2*sqrt(7) + 8) + 2*sqrt(7)*sqrt(2*sqrt(7) + 8))^2))

長い式は簡約化できる。

res2 |> factor |> println

   4*(-6 + sqrt(2)*sqrt(sqrt(7) + 4) + sqrt(14)*sqrt(sqrt(7) + 4))/(6 + sqrt(32*sqrt(7) + 128))

乙円の直径は 2.83398951148328 である。

2res2.evalf() |> println

   2.83398951148328

その他のパラメータは以下のとおりである。

   h = 8.861;  r1 = 2;  r2 = 1.41699;  x3 = 3.36716;  y3 = 2.61716

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1/2
   (h, r1, x3, y3) = (2*sqrt(2*sqrt(7) + 8)/3 + 2*sqrt(7)*sqrt(2*sqrt(7) + 8)/3, 2, 3*sqrt(2*sqrt(7) + 8)/8 + 2, 3*sqrt(7)/8 + 13/8)
   r2 = 4*(-6 + sqrt(2)*sqrt(sqrt(7) + 4) + sqrt(14)*sqrt(sqrt(7) + 4))/(6 + sqrt(32*sqrt(7) + 128))
   @printf("乙円の直径 = %g\n", 2r2)
   @printf("h = %g;  r1 = %g;  r2 = %g;  x3 = %g;  y3 = %g\n", h, r1, r2, x3, y3)
   plot([2r1, 2r1, -2r1, -2r1, 2r1], [0, 2r1, 2r1, 0, 0], color=:black, lw=0.5)
   circle2(r1, r1, r1)
   circle(0, 2r1 + r2, r2, :blue)
   circle(0, r3, r3, :green)
   circle2(x3, y3, r3, :green)
   plot!([-2r1, 0, 2r1], [0, h, 0], color=:orange, lw=0.5)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(x3, y3, "丙円:r3,(x3,y3) ", :black, :right, :vcenter)
       point(0, r3, " 丙円:r3,(0,r3)", :black, :left, :vcenter)
       point(r1, r1, "甲円:r1,(r1,r1)", :red, :center, delta=-delta/2)
       point(0, 2r1 + r2, "乙円:r2,(0,2r1+r2)", :blue, :center, delta=-delta/2)
       point(0, h, " h", :orange, :left, :vcenter)
       point(2r1, 2r1, "(2r1,2r1)", :black, :right, :bottom, delta=delta/2)
   end
end;