算額（その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;

算額（その486）

算額（その486）

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

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

算額の破損のため，図以外の情報は殆どない。
等円 2 個，甲円 1 個，乙円と丙円がそれぞれ 4 個が図のように配置されている。

等円の半径と中心座標を r1, (r1, r1)
甲円の半径と中心座標を r2, (0, y2)
乙円の半径と中心座標を r3, (x3, r1 + r3), (x3, r1 - r3)
丙円の半径と中心座標を r4, (r4, r4), (3r4, r4)
とおき，以下の連立方程式を解く。
未知数が r1, r2, y2, r3, x3, r4 の 6 個で，条件式が 3 個しかないので，未知数のうち 3 個を既知としなければ解は求まらない。
推測であるが，等円，甲円，乙円の直径が与えられていたのであろう。図を描くために y2, x3 を求めるが，丙円の直径を求めるだけであれば，甲円，乙円の情報は不要である。
条件式はそれぞれが独立なので，eq1, eq2, eq3 を連立方程式ではなく単一の方程式として解くことができる。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, y2::positive, 
     r3::positive, x3::positive, r4::positive;

eq1 = r1^2 + (y2 - r1)^2 - (r1 + r2)^2
eq2 = (r1 - x3)^2 + r3^2 - (r1 - r3)^2
eq3 = (r1 - 3r4)^2 + (r1 - r4)^2 - (r1 + r4)^2;

eq1 は等円と甲円についての方程式で，甲円の大きさは任意であり，その大きさに応じて中心座標 y2 が変わる。中心座標 y2 は以下のように求まる。2番目のものが適解である。

res1 = solve(eq1, y2)
res1 |> println

   Sym[r1 - sqrt(r2)*sqrt(2*r1 + r2), r1 + sqrt(r2)*sqrt(2*r1 + r2)]

eq2 は等円と乙円についての方程式で，乙円の直径は等円の直径の半分以下であれば任意であり，その大きさに応じて中心座標 x3 が変わる。中心座標 x3 は以下のように求まる。最初のものが適解である。

res2 = solve(eq2, x3)
res2 |> println

   Sym[-sqrt(r1)*sqrt(r1 - 2*r3) + r1, sqrt(r1)*sqrt(r1 - 2*r3) + r1]

eq3 は等円と丙円が外接するときの方程式で，等円の大きさに応じて丙円の大きさが決まる。丙円の大きさは等円の大きさの 1/9 である。

res3 = solve(eq3, r4)
res3 |> println

   Sym[r1/9, r1]

以上をまとめると，甲円と乙円は制限はあるものの任意に決めることができる。等円の大きさが決まれば丙円の大きさも決まる。

例えば，等円，甲円，乙円の半径を 90，70，30 とすれば，
y2 = 222.28756555322954
x3 = 38.03847577293368
r3 = 10.0
である。

res1[2] |> println
res2[1] |> println
res3[1] |> println

   r1 + sqrt(r2)*sqrt(2*r1 + r2)
   -sqrt(r1)*sqrt(r1 - 2*r3) + r1
   r1/9

(r1, r2, r3) = (90, 70, 30)
# y2
r1 + sqrt(r2)*sqrt(2*r1 + r2) |> println
# x3
-sqrt(r1)*sqrt(r1 - 2*r3) + r1 |> println
# r4
r1/9 |> println

   222.28756555322954
   38.03847577293368
   10.0

using Plots

function draw(r1, r2, r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r4, y2, x3) = (r1/9, r1 + sqrt(r2)*sqrt(2*r1 + r2), -sqrt(r1)*sqrt(r1 - 2*r3) + r1)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  r4 = %g;  y2 = %g;  x3 = %g\n", r1, r2, r3, r4, y2, x3)
   plot()
   circle(r1, r1, r1, :blue)
   circle(-r1, r1, r1, :blue)
   circle(0, y2, r2)
   circle(x3, r1 + r3, r3, :green)
   circle(-x3, r1 + r3, r3, :green)
   circle(x3, r1 - r3, r3, :green)
   circle(-x3, r1 - r3, r3, :green)
   circle(r4, r4, r4, :orange)
   circle(-r4, r4, r4, :orange)
   circle(3r4, r4, r4, :orange)
   circle(-3r4, r4, r4, :orange)
   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(r1, r1, " 等円:r1,(r1,r1)", :blue, :left, :vcenter)
       point(0, y2, " 甲円:r2,(0,y2)", :red, :left, :vcenter)
       point(x3, r1+r3, " 乙円:r3,(x3,r1+r3)", :green, :left, :vcenter)
       point(x3, r1-r3, " 乙円:r3,(x3,r1-r3)", :green, :left, :vcenter)
       point(3r4, r4, " 丙円:r4,(3r4,r4)", :black, :left, :vcenter)
   else
       plot!(showaxis=false)
   end
end;