算額（その469）

算額（その469）

宮城県石巻市雄勝町 葉山神社 明治17年(1884)

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

半円内に等円が 4 個入っている。黒い部分の面積（黒積）はいかほどか。

include("julia-source.txt")

using SymPy

@syms r0::positive, r::positive, x1::negative, x2::positive,
     x3::negative, y3::positive;
@syms r0::positive, r::positive, x1::negative, x2::positive, x3::negative, y3::positive
eq1 = x1^2 + r^2 - (r0 - r)^2
eq2 = (x1 - x3)^2 + (y3 - r)^2 - 4r^2
eq3 = x3^2 + y3^2 - (r0 - r)^2
eq4 = (x2 + 2r)^2 + r^2 - (r0 - r)^2
eq5 = x1 + x2 - 2x3

res = solve([eq1, eq2, eq3, eq4, eq5], (r0, x1, x2, x3, y3))

   2-element Vector{NTuple{5, Sym}}:
    (r + sqrt(2)*r*sqrt(sqrt(2) + 2), -r*(1 + sqrt(2)), r*(-1 + sqrt(2)), -r, r + sqrt(2)*r)
    (-sqrt(2)*r*sqrt(sqrt(2) + 2) + r, -r*(1 + sqrt(2)), r*(-1 + sqrt(2)), -r, r + sqrt(2)*r)

2 組の解が求まるが，最初のものが適解である。

res[1]

   (r + sqrt(2)*r*sqrt(sqrt(2) + 2), -r*(1 + sqrt(2)), r*(-1 + sqrt(2)), -r, r + sqrt(2)*r)

等円の半径が 1/2 のとき，
r0 = 1.80656;  x1 = -1.20711;  x2 = 0.207107  x3 = -0.5;  y3 = 1.20711
である。

黒積は以下のようにして求めることができる。

扇形 OAB の面積は半径 r0 の円の 1/4 なので，π*r0^2/4 = π*r*(1 + sqrt(2√2 + 4))/4

@syms r

r0 = r*(1 + sqrt(2√Sym(2) + 4))
gray = PI*r0^2/4;

水色の三角形2個の面積は (2r*r/2)*2 = 2r^2

cyan = 2r^2;

赤い扇形の中心角は 157.5 度で，2 つ分の面積は (πr^2 * 157.5/360)*2 = π*r^2*(7/8)

red = PI*r^2*(7//8);

黄色い扇形の中心角は 135 度で，面積は πr^2*135/360 = π*r^2*(3/8)

yellow = PI*r^2*(3//8);

灰色の部分の面積は gray - cyan - red - yellow

S = gray - cyan - red - yellow;
S |> expand |> simplify

   r^2*(-2 + sqrt(2)*pi/2 + pi*sqrt(2*sqrt(2) + 4)/2)

等円の半径 r = 1/2 とすると，黒積は 1.08153252024683 である。

@syms r
S(r => 1/2).evalf() |> println

   1.08153252024683

術では
sqrt(1/2) を「位」とおけば，
((sqrt(位 + 1) + 位)*PI/4 - 1/2)^(2r)^2 で黒積が求められる。

位 = sqrt(1/2)
円積率 = pi/4
等円の直径 = 1
((sqrt(位 + 1) + 位)*円積率 - 1/2)*等円の直径^2

   1.0815325202468267

using Plots

function circlef2(ox, oy, r, color=:red; beginangle=0, endangle=360, by=0.5, n=0, alpha=0.4)
   n != 0 && (by = (endangle - beginangle) / n)
   θ = beginangle:by:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   append!(x, [0, x[1]])
   append!(y, [0, y[1]])
   plot!(ox .+ x, oy .+ y, linecolor=:black, linewidth=0.5, seriestype=:shape, fillcolor=color, alpha=alpha)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   (r0, x1, x2, x3, y3) = (
       r*(1 + sqrt(2√2 + 4)),
       r*(-√2 - 1),
       r*(√2 - 1),
       -r,
       r*(1 + √2))
   @printf("r0 = %g;  x1 = %g;  x2 = %g  x3 = %g;  y3 = %g\n", r0, x1, x2, x3, y3)
   plot()
   circle(0, 0, r0, :blue, beginangle=0, endangle=180)
   circlef2(0, 0, r0, :black, beginangle=22.5, endangle=112.5)
   circle(x1, r, r)
   circle(x2, r, r)
   circlef2(x2, r, r, :orange, beginangle=0, endangle=135, alpha=0.8)
   circle(x2 + 2r, r, r)
   circlef2(x2 + 2r, r, r, beginangle=22.5, endangle=180, alpha=0.8)
   circle(x3, y3, r)
   circlef2(x3, y3, r, beginangle=-45, endangle=112.5, alpha=0.8)
   segment(0, 0, r0*x3/sqrt(x3^2 + y3^2), r0*y3/sqrt(x3^2 + y3^2))
   segment(0, 0, r0*y3/sqrt(x3^2 + y3^2), -r0*x3/sqrt(x3^2 + y3^2))
   segment(x2, r, x3, y3)
   segment(x2, r, x2 + 2r, r)
   segment(0, 0, x2, r)
   plot!([0, x2 + 2r, x2, 0], [0, r, r, 0], linewidth=0.5, seriestype=:shape, fillcolor=:cyan, alpha=0.8)
   plot!([0, x2, x3, 0], [0, r, y3, 0], linewidth=0.5, seriestype=:shape, fillcolor=:cyan, alpha=0.8)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
       plot!(ylims=(-0.1r0, r0))
       point(r0*x3/sqrt(x3^2 + y3^2), r0*y3/sqrt(x3^2 + y3^2), "A", :blue, :right, :bottom, delta=delta)
       point(r0*y3/sqrt(x3^2 + y3^2), -r0*x3/sqrt(x3^2 + y3^2), " B", :blue, :left, :vcenter)
       point(0, 0, "O ", :blue, :right, :top, delta=-delta)
       point(x1, r, "(x1,r)", :red)
       point(x2, r, "(x2,r)", :white, :left, :bottom, delta=delta)
       point(x2 + 2r, r, "(x2+2r,r)", :white, :right, :bottom, delta=delta)
       point(x3, y3, " (x3,y3)", :white, :left, :vcenter)
   else
      plot!(showaxis=false)
   end
end;