算額（その489）

算額（その489）

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

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

算額の破損のため，図以外の情報は殆どない。
外円を 3 本の弦で 5 個の領域に区切り，各領域に 甲円 2 個，乙円 3 個が入っている。

外円の半径と中心座標を r0, (0, 0); 一般性を失わずに r0 = 1 と設定できる。
甲円の半径と中心座標を r1, (x1, y1)
中円の半径と中心座標を r2, (0, r0 - r2), (0, 3r2 - r0), (0, r2 - r0)
右上がりの斜線と外円の交点座標を (sqrt(r0^2 - b^2), b)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms b::positive, r0::positive, r1::positive, 
     x1::positive, y1::positive, r2::positive;

r0 = 1
eq1 = x1^2 + y1^2 - (r0 - r1)^2
eq2 = distance(sqrt(r0^2 -(2r2 - r0)^2), 2r2 - r0, -sqrt(r0^2 - b^2), b, x1, y1) - r1^2
eq3 = distance(-sqrt(r0^2 -(2r2 - r0)^2), 2r2 - r0, sqrt(r0^2 - b^2), b, x1, y1) - r1^2
eq4 = distance(sqrt(r0^2 -(2r2 - r0)^2), 2r2 - r0, -sqrt(r0^2 - b^2), b, 0, 3r2 - r0) - r2^2
eq5 = distance(sqrt(r0^2 -(2r2 - r0)^2), 2r2 - r0, -sqrt(r0^2 - b^2), b, 0, r0 - r2) - r2^2;

# res = solve([eq1, eq2, eq3, eq4, eq5], (b, r1, x1, y1, r2))

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, y1, r2) = u
   x = sqrt(1 - b^2)
   y = sqrt(1 - r2)
   z = sqrt(r2)
   return [
       x1^2 + y1^2 - (1 - r1)^2,  # eq1
       -r1^2 + (x1 - (x1*(-2*b^2*r2 + b^2 + 2*b*z*x*y + 4*b*r2^2 - 4*b*r2 + 2*b - 4*r2^(3/2)*x*y + 2*z*x*y - 2*r2 + 1) - (b - 2*r2 + 1)*(-2*b^2*z*y + b^2*x1 + 4*b*r2^(3/2)*y + 2*b*z*y1*y - 2*b*z*y - 4*b*r2*x1 - 2*b*r2*x + 2*b*x1 + b*y1*x + b*x - 4*r2^(3/2)*y1*y + 2*z*y1*y + 4*r2^2*x1 + 4*r2^2*x - 4*r2*x1 - 2*r2*y1*x - 4*r2*x + x1 + y1*x + x)/2)/(-2*b^2*r2 + b^2 + 2*b*z*x*y + 4*b*r2^2 - 4*b*r2 + 2*b - 4*r2^(3/2)*x*y + 2*z*x*y - 2*r2 + 1))^2 + (y1 - (y1*(-2*b^2*r2 + b^2 + 2*b*z*x*y + 4*b*r2^2 - 4*b*r2 + 2*b - 4*r2^(3/2)*x*y + 2*z*x*y - 2*r2 + 1) - (x + sqrt(1 - (2*r2 - 1)^2))*(-2*b^2*z*y + b^2*x1 + 4*b*r2^(3/2)*y + 2*b*z*y1*y - 2*b*z*y - 4*b*r2*x1 - 2*b*r2*x + 2*b*x1 + b*y1*x + b*x - 4*r2^(3/2)*y1*y + 2*z*y1*y + 4*r2^2*x1 + 4*r2^2*x - 4*r2*x1 - 2*r2*y1*x - 4*r2*x + x1 + y1*x + x)/2)/(-2*b^2*r2 + b^2 + 2*b*z*x*y + 4*b*r2^2 - 4*b*r2 + 2*b - 4*r2^(3/2)*x*y + 2*z*x*y - 2*r2 + 1))^2,  # eq2
       -r1^2 + (x1 - (x1*(-2*b^2*r2 + b^2 + 2*b*z*x*y + 4*b*r2^2 - 4*b*r2 + 2*b - 4*r2^(3/2)*x*y + 2*z*x*y - 2*r2 + 1) - (b - 2*r2 + 1)*(2*b^2*z*y + b^2*x1 - 4*b*r2^(3/2)*y - 2*b*z*y1*y + 2*b*z*y - 4*b*r2*x1 + 2*b*r2*x + 2*b*x1 - b*y1*x - b*x + 4*r2^(3/2)*y1*y - 2*z*y1*y + 4*r2^2*x1 - 4*r2^2*x - 4*r2*x1 + 2*r2*y1*x + 4*r2*x + x1 - y1*x - x)/2)/(-2*b^2*r2 + b^2 + 2*b*z*x*y + 4*b*r2^2 - 4*b*r2 + 2*b - 4*r2^(3/2)*x*y + 2*z*x*y - 2*r2 + 1))^2 + (y1 - (y1*(-2*b^2*r2 + b^2 + 2*b*z*x*y + 4*b*r2^2 - 4*b*r2 + 2*b - 4*r2^(3/2)*x*y + 2*z*x*y - 2*r2 + 1) + (x + sqrt(1 - (2*r2 - 1)^2))*(2*b^2*z*y + b^2*x1 - 4*b*r2^(3/2)*y - 2*b*z*y1*y + 2*b*z*y - 4*b*r2*x1 + 2*b*r2*x + 2*b*x1 - b*y1*x - b*x + 4*r2^(3/2)*y1*y - 2*z*y1*y + 4*r2^2*x1 - 4*r2^2*x - 4*r2*x1 + 2*r2*y1*x + 4*r2*x + x1 - y1*x - x)/2)/(-2*b^2*r2 + b^2 + 2*b*z*x*y + 4*b*r2^2 - 4*b*r2 + 2*b - 4*r2^(3/2)*x*y + 2*z*x*y - 2*r2 + 1))^2,  # eq3
       -r2^2 + (b - 2*r2 + 1)^2*(2*b^2*z*y - 10*b*r2^(3/2)*y + 4*b*z*y - b*r2*x + 12*r2^(5/2)*y - 10*r2^(3/2)*y + 2*z*y + 2*r2^2*x - r2*x)^2/(4*(-2*b^2*r2 + b^2 + 2*b*z*x*y + 4*b*r2^2 - 4*b*r2 + 2*b - 4*r2^(3/2)*x*y + 2*z*x*y - 2*r2 + 1)^2) + (3*r2 - 1 - (b^3*r2 + 2*b^2*z*x*y - 18*b^2*r2^2 + 15*b^2*r2 - 2*b^2 + 44*b*r2^3 - 60*b*r2^2 + 27*b*r2 - 4*b - 8*r2^(5/2)*x*y + 8*r2^(3/2)*x*y - 2*z*x*y - 24*r2^4 + 44*r2^3 - 34*r2^2 + 13*r2 - 2)/(2*(-2*b^2*r2 + b^2 + 2*b*z*x*y + 4*b*r2^2 - 4*b*r2 + 2*b - 4*r2^(3/2)*x*y + 2*z*x*y - 2*r2 + 1)))^2,  # eq4
       -r2^2 + (-r2 + 1 - (b^3*z*y + 5*b^2*r2^(3/2)*y - 3*b^2*z*y + 3*b^2*r2*x/2 - b^2*x + 2*b*r2^(3/2)*y - 3*b*z*y + 2*b*r2^2*x - b*r2*x - b*x + 4*r2^(7/2)*y - 8*r2^(5/2)*y + r2^(3/2)*y + z*y + 6*r2^3*x - 10*r2^2*x + 7*r2*x/2)/(2*b^2*z*y + 4*b*r2^(3/2)*y - 2*b*z*y + 2*b*r2*x - b*x - 4*z*y + 4*r2^2*x - 4*r2*x - x))^2 + (3*b^3*r2/2 - b^3 - b^2*z*x*y - b^2*r2^2 + 3*b^2*r2/2 - b^2 - 2*b*r2^(3/2)*x*y + 2*b*z*x*y - 2*b*r2^3 + 2*b*r2^2 - 3*b*r2/2 + b + 8*r2^(5/2)*x*y - 10*r2^(3/2)*x*y + 3*z*x*y - 4*r2^4 + 10*r2^3 - 5*r2^2 - 3*r2/2 + 1)^2/(2*b^2*z*y + 4*b*r2^(3/2)*y - 2*b*z*y + 2*b*r2*x - b*x - 4*z*y + 4*r2^2*x - 4*r2*x - x)^2,  # eq5
   ]
end;
r0 = 1
iniv = BigFloat[0.77, 0.32, 0.55, 0.3, 0.31]
res = nls(H, ini=iniv)

   (BigFloat[0.7784615384615384615384615384615384615384615384615384615384616011678188086041049, 0.3461538461538461538461538461538461538461538461538461538461544534842427177340537, 0.576923076923076923076923076923076923076923076923076923076925033822024574673162, 0.3076923076923076923076923076923076923076923076923076923076934100784216705502225, 0.3076923076923076923076923076923076923076923076923076923076922943346533297705356], true)

   b = 0.778462;  r1 = 0.346154;  x1 = 0.576923;  y1 = 0.307692;  r2 = 0.307692

using Plots

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

算額（その488）

算額（その488）

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

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

算額の破損のため，図以外の情報は殆どない。
扇の中に正三角形と大円，中円，小円が入っている。

外円の半径と中心座標を r0, (0, 0)
大円の半径と中心座標を r1, (0, 0); r1 = r0*sin(pi/6)
中円の半径と中心座標を r2, (0, r0 - r2); r2 = r0(1 - sin(pi/6))/2
小円の半径と中心座標を r3, (x3, y3)
下部の円弧の半径と中心座標を r4, (0, -r0); r4 = r0 - r1
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive, 
     r3::positive, x3::positive, y3::negative, r4::positive;

x = r0*cos(PI/6)
y = r0*sin(PI/6)
r1 = y
r2 = (r0 - y)/2
r4 = r0 - r1
eq1 = distance(0, -r0, x, y, x3, y3) - r3^2
eq2 = x3^2 + (y3 + r0)^2 - (r3 + r4)^2
eq3 = x3^2 + y3^2 - (r0 - r3)^2
res = solve([eq1, eq2, eq3], (r3, x3, y3))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (r0*(1 + 3*sqrt(5))/32, r0*(5*sqrt(3) + 7*sqrt(15))/64, r0*(-53/64 + 9*sqrt(5)/64))

外円と円弧の交点座標 (x0, y0) を求める。

@syms x0::positive, y0::negative
eq1 = x0^2 + y0^2 - r0^2
eq2 = x0^2 + (r0 + y0)^2 - (r0 - y)^2
res2 = solve([eq1, eq2], (x0, y0))

   1-element Vector{Tuple{Sym, Sym}}:
    (sqrt(15)*r0/8, -7*r0/8)

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 1
   x = r0*cos(pi/6)
   y = r0*sin(pi/6)  
   r1 = y
   r2 = (r0 - y)/2
   r4 = r0 - r1
   (r3, x3, y3) = (r0*(1 + 3*sqrt(5))/32, r0*(5*sqrt(3) + 7*sqrt(15))/64, r0*(-53/64 + 9*sqrt(5)/64))
   (x0, y0) = (sqrt(15)*r0/8, -7*r0/8)
   θ = atand((r0 + y0)/x0)
   @printf("r0 = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  y3 = %g;  r4 = %g\n",
       r0, r1, r2, r3, x3, y3, r4)
   plot([0, x, -x, 0], [-r0, y, y, -r0], color=:black, lw=0.5)
   circle(0, 0, r0, :blue)
   circle(0, 0, r1, :green)
   circle(0, r0 - r2, r2, :orange)
   circle(x3, y3, r3, :magenta)
   circle(-x3, y3, r3, :magenta)
   circle(0, -r0, r4, beginangle=θ, endangle=181-θ)
   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, y, "(x,y)  ", :blue, :right, :bottom, delta=delta/2)
       point(x0, y0, "(x0,y0)")
       point(r0, 0, "r0 ", :blue, :right, :bottom, delta=delta/2)
       point(0, r1, " r1", :green, :left, delta=-delta/2)
       point(0, -r1, " -r1", :green, :left, :bottom, delta=delta)
       point(0, r0 - r2, " 中円:r2,(0,r0-r2)", :black, :left, :vcenter)
       point(r1/2, r1/2, "大円:r1,(0,0)", :green, :center, delta=-2delta, mark=false)
       point(x3, y3, "小円:r3,(x3,y3)", :magenta, :center, delta=-delta/2)
   else
       plot!(showaxis=false)
   end
end;