算額（その487）

算額（その487）

宮城県丸森町小斎日向 鹿島神社 大正年間

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

算額の破損のため，図以外の情報は殆どない。
外円の中に弦と斜線があり，区切られた領域に甲円が 4 個，乙円が 2 個入っている。

外円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (0, r0 - r1), (x1, a - r1), (0, r1 - r0)
乙円の半径と中心座標を r2, (x2, a + r2)
弦と y 軸の交点の y 座標を a
斜線と円の交点の y 座標を -b
とおき，以下の連立方程式を解く。

条件式が5個で，変数が b, r0, r1, x1, r2, x2 の 6 個なので，どれか一つを既知とすれば解を求めて図形を書くことができる。r0 = 1 としても一般性を損なわない（r0 が異なる図形はすべて相似）。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r0::positive, 
     r1::positive, x1::positive, r2::positive, x2::positive;
r0 = 1
a = r0 - 2r1
eq1 = distance(0, a, sqrt(r0^2 - b^2), -b, 0, r1 - r0) - r1^2
eq2 = distance(0, a, sqrt(r0^2 - b^2), -b, x1, a - r1) - r1^2
eq3 = x1^2 + (a - r1)^2 - (r0 - r1)^2
eq4 = x2^2 + (a + r2)^2 - (r0 - r2)^2
eq5 = x2^2 + (r0 - r1 - a - r2)^2 - (r1 + r2)^2;

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (b, r1, x1, r2, x2) = u
   return [
       -r1^2 + (1 - b^2)*(3*r1 - 2)^2*(b - 2*r1 + 1)^2/(-4*b*r1 + 2*b + 4*r1^2 - 4*r1 + 2)^2 + (r1 - ((b^2 - 1)*(3*r1 - 2)/2 + (r1 - 1)*(-2*b*r1 + b + 2*r1^2 - 2*r1 + 1))/(-2*b*r1 + b + 2*r1^2 - 2*r1 + 1) - 1)^2,  # eq1
       -r1^2 + (x1 - (-b^2*x1 + b*r1*sqrt(1 - b^2) - 2*r1^2*sqrt(1 - b^2) + r1*sqrt(1 - b^2) + x1)/(2*(-2*b*r1 + b + 2*r1^2 - 2*r1 + 1)))^2 + (-3*r1 - (sqrt(1 - b^2)*(-b*x1 + 2*r1*x1 + r1*sqrt(1 - b^2) - x1)/2 + (1 - 3*r1)*(-2*b*r1 + b + 2*r1^2 - 2*r1 + 1))/(-2*b*r1 + b + 2*r1^2 - 2*r1 + 1) + 1)^2,  # eq2
       x1^2 + (1 - 3*r1)^2 - (1 - r1)^2,  # eq3
       x2^2 - (1 - r2)^2 + (-2*r1 + r2 + 1)^2,  # eq4
       x2^2 + (r1 - r2)^2 - (r1 + r2)^2,  # eq5
   ]
end;

iniv = BigFloat[0.85, 0.38, 0.61, 0.24, 0.60]
res = nls(H, ini=iniv)

   (BigFloat[0.847036874527671344254215992490411686652141768831937314871633149986214879073071, 0.3787321874818335132405467689419389110250823757243052813412166596148529968650207, 0.6061552534203894328837729024835844965150345459568913206939271172607303312643575, 0.2352941176470588235294117647058823529411764705882352941176470588235292092020814, 0.5970375394498355116094971924511155185714872080275178420141772660186381316118884], true)

外円の半径を 1 としたときの，図を描くためのパラーメータは以下の通り。

a = 0.242536;  b = 0.847037;  r0 = 1;  r1 = 0.378732;  x1 = 0.606155;  r2 = 0.235294;  x2 = 0.597038

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (b, r1, x1, r2, x2) = res[1] # [55.0, 20, 35, 12, 35]
   r0 = 1
   a = r0 - 2r1
   @printf("a = %g;  b = %g;  r0 = %g;  r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g\n",
       a, b, r0, r1, x1, r2, x2)
   plot()
   circle(0, 0, r0, :black)
   circle(0, r0 - r1, r1, :blue)
   circle(0, r1 - r0, r1, :blue)
   circle(x1, a - r1, r1, :blue)
   circle(-x1, a - r1, r1, :blue)
   circle(x2, a + r2, r2)
   circle(-x2, a + r2, r2)
   segment(-sqrt(r0^2 - a^2), a, sqrt(r0^2 - a^2), a, :green)
   segment(0, a, sqrt(r0^2 - b^2), -b, :red)
   segment(0, a, -sqrt(r0^2 - b^2), -b, :red)
   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, a, " a")
       point(0, -b, " -b", :black, :left, :vcenter)
       point(sqrt(r0^2 - b^2), -b, "(√(r0^2-b^2),b)", :black, :left, :top)
       point(0, r1 - r0, " r1-r0", :blue, :left, :vcenter)
       point(x1, a - r1, " 甲円:r1,(x1,a-r1)", :blue, :center, :top, delta=-delta/2)
       point(x2, a + r2, "乙円:r2\n(x2,a+r2)", :red, :center, :vcenter)
       point(0, r0 - r1, " r0-r1", :blue, :left, :vcenter)
       point(0, r0, " r0", :black, :left, :bottom, delta=delta)
   else
       plot!(showaxis=false)
   end
end;