算額（その470）

算額（その470）

宮城県角田市横倉 愛宕神社 明治15年(1882)1月
http://www.wasan.jp/miyagi/yokokuraatago.html

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

甲円内に弦を挟んで乙円 1 個，丙円 1 個，丁円 2 個が入っている。
乙円，丙円，丁円の直径がそれぞれ 4 寸，3 寸，1 寸のとき，甲円の直径を求めよ。

甲円の半径と中心座標を r1, (0, 0)
乙円の半径と中心座標を r2, (x2, 2r4 - r1 + r2)
丙円の半径と中心座標を r3, (x3, 2r4 - r1 + r3)
丁円の半径と中心座標を r4, (0, r4 - r1), (x4, 3r4 - r1)
とし，以下のように記号を定め方程式を解く。

include("julia-source.txt")

using SymPy
@syms r1::positive, r2::positive, x2::positive,
     r3::positive, x3::negative,
     r4::positive, x4::negative;
eq1 = x2^2 + (2r4 - r1 + r2)^2 - (r1 - r2)^2
eq2 = x3^2 + (2r4 - r1 + r3)^2 - (r1 - r3)^2
eq3 = (x2 - x4)^2 + (r2 - r4)^2 - (r2 + r4)^2
eq4 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2
solve([eq1, eq2, eq3, eq4], (r1, x2, x3, x4))

   2-element Vector{NTuple{4, Sym}}:
    (r2 + r3 + r4, 2*sqrt(r4)*(-sqrt(r2)*sqrt(r3) + r3)*sqrt(2*sqrt(r2)*sqrt(r3) + r2 + r3)/(r2 - r3), 2*sqrt(r4)*(sqrt(r2)*sqrt(r3) - r2)*sqrt(2*sqrt(r2)*sqrt(r3) + r2 + r3)/(r2 - r3), -2*sqrt(r4)*sqrt(2*sqrt(r2)*sqrt(r3) + r2 + r3))
    (r2 + r3 + r4, 2*sqrt(r4)*(sqrt(r2)*sqrt(r3) + r3)*sqrt(-2*sqrt(r2)*sqrt(r3) + r2 + r3)/(r2 - r3), -2*sqrt(r4)*(sqrt(r2)*sqrt(r3) + r2)*sqrt(-2*sqrt(r2)*sqrt(r3) + r2 + r3)/(r2 - r3), -2*sqrt(r4)*sqrt(-2*sqrt(r2)*sqrt(r3) + r2 + r3))

2 組の解が得られるが，2 組目のものが適解である。

   r1 = 4;  x2 = 1.73205;  x3 = -2;  x4= -0.267949
   甲円の直径 = 8

甲円の半径は r2 + r3 + r4 すなわち，乙円，丙円，丁円の半径の和である。
乙円，丙円，丁円の直径がそれぞれ 4 寸，3 寸，1 寸ならば，甲円の直径は 4 + 3 + 1 = 8 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, r4) = (4, 3, 1) .// 2
   (r1, x2, x3, x4) = (r2 + r3 + r4, 2*sqrt(r4)*(sqrt(r2)*sqrt(r3) + r3)*sqrt(-2*sqrt(r2)*sqrt(r3) + r2 + r3)/(r2 - r3), -2*sqrt(r4)*(sqrt(r2)*sqrt(r3) + r2)*sqrt(-2*sqrt(r2)*sqrt(r3) + r2 + r3)/(r2 - r3), -2*sqrt(r4)*sqrt(-2*sqrt(r2)*sqrt(r3) + r2 + r3))
   @printf("r1 = %g;  x2 = %g;  x3 = %g;  x4= %g\n", r1, x2, x3, x4)
   @printf("甲円の直径 = %g\n", 2r1)
   plot()
   circle(0, 0, r1, :red)
   circle(x2, 2r4 - r1 + r2, r2, :magenta)
   circle(x3, 2r4 - r1 + r3, r3, :blue)
   circle(x4, 3r4 - r1, r4, :green)
   circle(0, r4 - r1, r4, :green)
   segment(-sqrt(r1^2 - (2r4 - r1)^2), 2r4 - r1, sqrt(r1^2 - (2r4 - r1)^2), 2r4 - r1)
   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, r1, " r1", :red, :left, :bottom, delta=delta/2)
       point(x2, 2r4 - r1 + r2, "乙円:r2,(x2,2r4-r1+r2)", :magenta, :center, :top, delta=-delta)
       point(x3, 2r4 - r1 + r3, "丙円:r3,(x3,2r4-r1+r3)", :blue, :center, :top, delta=-delta)
       point(0, r4 - r1, " 丁円:r4\n (0,r4-r1)", :black, :left, :vcenter)
       point(x4, 3r4 - r1, "丁円:r4,(x4,3r4-r1) ", :black, :right, :vcenter)
       point(0.2r1, 0.5r1, " 甲円:r1,(0,0)", :red, :left, :vcenter)
       point(0, 2r4 - r1, " 2r4-r1", :black, :left, :bottom, delta=delta)
   end
end;

 

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