算額（その895）

算額（その895）

七二 加須市大字外野 棘脱地蔵堂 明治7年(1874)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

仁円 1 個，義円 2 個，禮円 4 個と 2 本の斜線がある。義円の直径が 3 寸のとき，禮円の直径はいかほどか。

仁円の半径と中心座標を R, (0, 0)
義円の半径と中心座標を r1, (0, R + r1), (0, R - r1)
禮円の半径と中心座標を r2, (x21, y21), (x22, y21)
斜線と人円の交点座標を (x, y)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms x::positive, y::negative, R::positive, r1::positive, r2::positive,
     x21::positive, y21::positive, x22::positive, y22::positive
eq1 = dist2(0, R + 2r1, x, y,0, R - r1, r1)
eq2 = dist2(0, R + 2r1, x, y, x21, y21, r2)
eq3 = dist2(0, R + 2r1, x, y, x22, y22, r2)
eq4 = x21^2 + y21^2 - (R - r2)^2
eq5 = x22^2 + (y22 - R - r1)^2 - (r1 - r2)^2
eq6 = y21/x21 - x/(R + 2r1 - y)
eq7 = (y22 - R - r1)/x22 - x/(R + 2r1 - y)
eq8 = x^2 + y^2 - R^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 Float64.(v), r.f_converged
end;

function H(u)
   (R, r2, x21, y21, x22, y22, x, y) = u
   return [
       -R^2*r1^2 - 4*R*r1^3 + 2*R*r1^2*y - 4*r1^4 + 4*r1^3*y + 8*r1^2*x^2 - r1^2*y^2,  # eq1
       -R^2*r2^2 + R^2*x^2 - 2*R^2*x*x21 + R^2*x21^2 - 4*R*r1*r2^2 + 4*R*r1*x^2 - 8*R*r1*x*x21 + 4*R*r1*x21^2 + 2*R*r2^2*y - 2*R*x^2*y21 + 2*R*x*x21*y + 2*R*x*x21*y21 - 2*R*x21^2*y - 4*r1^2*r2^2 + 4*r1^2*x^2 - 8*r1^2*x*x21 + 4*r1^2*x21^2 + 4*r1*r2^2*y - 4*r1*x^2*y21 + 4*r1*x*x21*y + 4*r1*x*x21*y21 - 4*r1*x21^2*y - r2^2*x^2 - r2^2*y^2 + x^2*y21^2 - 2*x*x21*y*y21 + x21^2*y^2,  # eq2
       -R^2*r2^2 + R^2*x^2 - 2*R^2*x*x22 + R^2*x22^2 - 4*R*r1*r2^2 + 4*R*r1*x^2 - 8*R*r1*x*x22 + 4*R*r1*x22^2 + 2*R*r2^2*y - 2*R*x^2*y22 + 2*R*x*x22*y + 2*R*x*x22*y22 - 2*R*x22^2*y - 4*r1^2*r2^2 + 4*r1^2*x^2 - 8*r1^2*x*x22 + 4*r1^2*x22^2 + 4*r1*r2^2*y - 4*r1*x^2*y22 + 4*r1*x*x22*y + 4*r1*x*x22*y22 - 4*r1*x22^2*y - r2^2*x^2 - r2^2*y^2 + x^2*y22^2 - 2*x*x22*y*y22 + x22^2*y^2,  # eq3
       x21^2 + y21^2 - (R - r2)^2,  # eq4
       x22^2 - (r1 - r2)^2 + (-R - r1 + y22)^2,  # eq5
       -x/(R + 2*r1 - y) + y21/x21,  # eq6
       -x/(R + 2*r1 - y) + (-R - r1 + y22)/x22,  # eq7
       -R^2 + x^2 + y^2,  # eq8
   ]
end;

r1 = 6/2
iniv = BigFloat[6.1, 1.1, 4.8, 1.6, 1.8, 9.6, 5.2, -2.8]
res = nls(H, ini=iniv)

   ([6.0, 1.0, 4.714045207910317, 1.6666666666666667, 1.8856180831641267, 9.666666666666666, 5.261948151328113, -2.883036880224506], true)

禮円の半径は義円の半径の 1/3 である。
義円の直径が 6 寸のとき，禮円の直径は 2 寸である。

ちなみに，仁円の直径は義円の直径の 2 倍の 12 寸である。その他のパラメータは以下のとおりである。
R = 6;  r2 = 1;  x21 = 4.71405;  y21 = 1.66667;  x22 = 1.88562;  y22 = 9.66667;  x = 5.26195;  y = -2.88304

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 6/2
   (R, r2, x21, y21, x22, y22, x, y) = res[1]
   @printf("義円の直径が %g 寸のとき，禮円の直径は %g 寸\n", 2r1, 2r2)
   @printf("R = %g;  r2 = %g;  x21 = %g;  y21 = %g;  x22 = %g;  y22 = %g;  x = %g;  y = %g\n", R, r2, x21, y21, x22, y22, x, y)
   plot()
   circle(0, 0, R, :blue)
   circle(0, R - r1, r1)
   circle(0, R + r1, r1)
   circle2(x22, y22, r2, :green)
   circle2(x21, y21, r2, :green)
   segment(0, R + 2r1, x, y, :magenta)
   segment(0, R + 2r1, -x, y, :magenta)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R + 2r1, "R+2r1", :red, :center, :bottom, delta=delta/2)
       point(x, y, "(x,y) ", :blue, :right, :vcenter)
       point(0, 0, "仁円:R,(0,0)", :blue, :center, delta=-delta)
       point(0, R + r1, "義円:r1,(0,R+r1)", :red, :center, delta=-1.5delta)
       point(0, R - r1, "義円:r1,(0,R-r1)", :red, :center, delta=-delta)
       point(x21, y21, "禮円:r2,(x21,y21) ", :green, :right, :vcenter)
       point(x22, y22, " 禮円:r2,(x22,y22) ", :green, :left, :bottom)
   end
end;