算額（その511）

算額（その511）

長野県長野市安茂里正覚院 久保寺観音堂 文政13年(1830)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html

半円内に2本の斜線を引き，甲円 2 個と乙円を入れる。乙円の径が 1 寸のとき，甲円の径はいかほどか。

各円の半径，中心座標を以下のようにする。
外円: R, (0, 0)
甲円: r1, (x1, y1) および (x2, r1)
乙円: r3, (x3, y3)
外円の円周上にある2本の斜線の交点座標を (x, y)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms x1, y1, r1, x2, r3, x3, y3, x, R

r3 = 1//2
y = sqrt(R^2 - x^2)  # (x, y) は外円の円周上にある
eq1 = x1^2 + y1^2 - (R - r1)^2  # 右の甲円が外円に内接する
eq2 = x3^2 + y3^2 - (R - r3)^2  # 乙円が外円に内接する
eq3 = sqrt((R + x)^2 + y^2) + sqrt((R - x)^2 + y^2) - 2R - 2r1  # 鉤股弦に内接する円の直径
eq4 = (y3/x3)*(y/(x + R)) + 1
eq5 = numerator(distance(-R, 0, x, y, x2, r1) - r1^2) |> simplify
eq6 = numerator(distance(R, 0, x, y, x1, y1) - r1^2) |> simplify
eq7 = ((x - R)/2 - x3)^2 + (y/2 - y3)^2 - r3^2
eq8 = (y1/x1) * (y/(x - R)) + 1;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (x1, y1, r1, x2, x3, y3, x, R))

using Printf
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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (x1, y1, r1, x2, x3, y3, x, R) = u
   t = sqrt(R^2 - x^2)
   return [
       x1^2 + y1^2 - (R - r1)^2,  # eq1
       x3^2 + y3^2 - (R - 1/2)^2,  # eq2
       -2*R - 2*r1 + sqrt(R^2 - x^2 + (R - x)^2) + sqrt(R^2 - x^2 + (R + x)^2),  # eq3
       1 + y3*sqrt(R^2 - x^2)/(x3*(R + x)),  # eq4
       (-4*R^2*r1^2 + (R*r1 - R*t + r1*x - x2*t)^2 + (R^2 - R*x + R*x2 - r1*t - x*x2)^2)/(4*R^2),  # eq5
       (-4*R^2*r1^2 + (R*y1 - R*t - x*y1 + x1*t)^2 + (-R^2 - R*x + R*x1 + x*x1 + y1*t)^2)/(4*R^2),  # eq6
       (-y3 + t/2)^2 + (-R/2 + x/2 - x3)^2 - 1/4,  # eq7
       1 + y1*t/(x1*(-R + x)),  # eq8
   ]
end;

iniv = BigFloat[1.7, 4.2, 2.1, -3.5, -5.8, 2.4, -4.6, 6.4]
res = nls(H, ini=iniv)
println(res)

   (BigFloat[3.461538461538461538461538509169459561367959236632964526959234465030510845593508, 8.307692307692307692307692429999766197094115669634049393814154165178838427207838, 4.000000000000000000000000058695639811160390953981958412762403947086964336791015, -7.000000000000000000000000103675105612121327124114588925459978536567327493943071, -11.53846153846153846153846171778119379511106203957672690108302039508364569817559, 4.807692307692307692307692375277671136103845336879392455246078657573883943485054, -9.153846153846153846153846290689476553117442078831806765600936112963486841120747, 13.00000000000000000000000019241303418103091110455672021661459202029874280726291], true)

   r3 = 0.5;  x1 = 3.46154;  y1 = 8.30769; r1 = 4;  x2 = -7;  x3 = -11.5385;  y3 = 4.80769;  x = -9.15385;  y = 9.23077;  R = 13
   甲円の径 = 2r1 = 8

甲円の半径は 4 である。元の単位では甲円の直径は 8 寸である。

using Plots
function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1//2
   (x1, y1, r1, x2, x3, y3, x, R) = res[1]
   y = sqrt(R^2 - x^2)
   @printf("r3 = %g;  x1 = %g;  y1 = %g; r1 = %g;  x2 = %g;  x3 = %g;  y3 = %g;  x = %g;  y = %g;  R = %g\n",
       r3, x1, y1, r1, x2, x3, y3, x, y, R)
   @printf("甲円の径 = 2r1 = %g\n", 2r1)
   plot([R, x, -R, R], [0, y, 0, 0], color=:brown, lw=0.5)
   circle(0, 0, R, beginangle=0, endangle=180)
   circle(x1, y1, r1, :blue)
   circle(x2, r1, r1, :green)
   circle(x3, y3, r3, :magenta)
   if more == true
       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\n (x1,y1)", :blue, :left, :vcenter)
       point(x2, r1, " 甲円:r1\n (x2,r1)", :blue, :left, :top, delta=-delta)
       point(x3, y3, "  乙円:r3,(x3,y3)", :magenta, :left, :vcenter)
       point(x, y, "(x,y) ", :brown, :right, :vcenter)
   else
       plot!(showaxis=false)
   end
end;