算額（その494）

算額（その494）

北越長岡蔵王権現社（金峰神社） 寛政7年(1795)
涌田和芳，山田章，島津大，若山想思，市野梨保子，矢野敦大，山崎瑠聖: 長岡金峰神社の紛失算額，長岡工業高等専門学校研究紀要，第 55 巻，p.39-44，2019
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_55/55_39wakuta.pdf

長径，短径が a, b の楕円内に半径 r の円が 3 個入っている。左右の円は楕円と 2 点で接し中央の円と外接している。中央の円は楕円に内接している。

長径，短径が与えられたとき，円の半径を求めよ。

注: 和算では楕円の長径，短径は今でいう長径，短径の 2 倍である。同じく，円の径も直径を意味する。
ここでは他の記事と同じで，長径，短径は今の意味，円も半径で扱う。

eq2 は，山本賀前著『算法助術』の84番目の公式である。
山本賀前: 算法助術，天保12年(1841)，東北大学DB.
インターネットでのダウンロード・閲覧は
中村信弥編著:和算の図形公式-『算法助術』-，2008.
http://www.wasan.jp/kosiki/kosiki3.pdf

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r::positive, x::positive;

eq1 = x^2 + (b - r)^2 - 4r^2
eq2 = (b^2 - r^2)*(a^2 - b^2)/b^2 - x^2
res = solve([eq1, eq2], (r, x))

   1-element Vector{Tuple{Sym, Sym}}:
    ((-2*a^2*b^2 + 2*b^4 + 4*b^2*(a - b)*(a + b)*(a^2 + b^2)/(a^2 + 2*b^2))/(2*b*(a - b)*(a + b)), 2*b*sqrt(a - b)*sqrt(a + b)*sqrt(a^2 + b^2)/(a^2 + 2*b^2))

円の半径は以下のように簡約化される。これは，術で述べられている式と同じである。

res[1][1] |> simplify |> println

   a^2*b/(a^2 + 2*b^2)

右側の円の中心座標（x座標）も以下のように簡約化される。

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

   2*b*sqrt(a^4 - b^4)/(a^2 + 2*b^2)

長径と短径が 18, 4.5 のとき，円の半径は 4 である。

(a, b) = (36, 9) .// 2
a^2*b/(a^2 + 2*b^2) |> println

   4//1

   r = 4;  x = 7.98436;  a = 18;  b = 4.5

算額の図は，長径が短径の 2 倍程度なので上の図のような感じであるが，長径と短径が 18, 4.5 のときは，かなり横長なものになる。

using Plots

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x) = (a^2*b/(a^2 + 2b^2), 2b*sqrt(a^4 - b^4)/(a^2 + 2b^2))
   @printf("r = %g;  x = %g;  a = %g;  b = %g\n", r, x, a, b)
   plot()
   circle(0, b - r, r, :red)
   circle(x, 0, r, :red)
   circle(-x, 0, r, :red)
   ellipse(0, 0, a, b, color=:green)
   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(x, 0, "x", :red, :center, :bottom, delta=delta)
       point(a, 0, " a", :green, :left, :bottom, delta=delta)
       point(0, b - r, " b-r", :red, :left, :vcenter)
       point(0, b, " b", :green, :left, :bottom, delta=delta)
   else
       plot!(showaxis=false)
   end
end;