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