算額（その1384）

算額（その1384）

一八 大里郡岡部村岡 稲荷社 文化14年(1817)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円9個，外円，弦2本
#Julia, #SymPy, #算額, #和算

外円の中に 2 本の弦を引き，甲円 1 個，乙円 1 個，丙円 6 個を容れる。甲円，乙円の直径が 22 寸，11 寸のとき，丙円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, r1 - R)
乙円の半径と中心座標を r2, (0, R - r2)
丙円の半径と中心座標を r3, (x31, y31), (x32, y32), (x33, y33)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms R::positive, x01::positive, x02::positive,
     r1::positive, r2::positive, r3::positive,
     x31::positive, y31::positive,
     x32::positive, y32::positive,
     x33::positive, y33::positive;

R = r1 + r2
eq1 = x31^2 + y31^2 - (R - r3)^2
eq2 = x32^2 + y32^2 - (R - r3)^2
eq3 = x33^2 + y33^2 - (R - r3)^2
eq4 = (x31 - x32)^2 + (y31 - y32)^2 - 4r3^2
eq5 = (x31 - x33)^2 + (y31 - y33)^2 - 4r3^2
eq6 = dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), 0, r1 - R) - r1^2
eq7 = dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), 0, R - r2) - r2^2
eq8 = dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), x32, y32) - r3^2
eq9 = dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), x33, y33) - r3^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (x01, x02, r3, x31, y31, x32, y32, x33, y33))

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)
   (x01, x02, r3, x31, y31, x32, y32, x33, y33) = u
   return [
       x31^2 + y31^2 - (R - r3)^2,
       x32^2 + y32^2 - (R - r3)^2,
       x33^2 + y33^2 - (R - r3)^2,
       (x31 - x32)^2 + (y31 - y32)^2 - 4r3^2,
       (x31 - x33)^2 + (y31 - y33)^2 - 4r3^2,
       dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), 0, r1 - R) - r1^2,
       dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), 0, R - r2) - r2^2,
       dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), x32, y32) - r3^2,
       dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), x33, y33) - r3^2
   ]
end;
(r1, r2) = (22, 11) ./ 2
R = r1 + r2

iniv = BigFloat[4, 13, 3.1, 13, 4, 13, -1.5, 10, 10]
res = nls(H, ini=iniv)

   ([4.069279408445208, 13.215553020559287, 3.0, 12.727922061357855, 4.5, 13.420835425335575, -1.459854953773714, 9.520851253161299, 9.570966064884825], true)

甲円，乙円の直径が 22 寸, 11 寸のとき，丙円の直径は 6 寸である。

全てのパラメータは以下のとおりである。
   R = 16.5;  x01 = 4.06928;  x02 = 13.2156;  r1 = 11;  r2 = 5.5;  r3 = 3;  x31 = 12.7279;  y31 = 4.5;  x32 = 13.4208;  y32 = -1.45985;  x33 = 9.52085;  y33 = 9.57097

function draw(r1, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
  (x01, x02, r3, x31, y31, x32, y32, x33, y33) = res[1]
   R = r1 + r2
   @printf("甲円，乙円の直径が %g, %g のとき，丙円の直径は %g である。\n", 2r1, 2r2, 2r3)
   @printf("R = %g;  x01 = %g;  x02 = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x31 = %g;  y31 = %g;  x32 = %g;  y32 = %g;  x33 = %g;  y33 = %g\n",
       R, x01, x02, r1, r2, r3, x31, y31, x32, y32, x33, y33)
   plot()
   circle(0, 0, R, :magenta)
   circle(0, r1 - R, r1, :blue)
   circle(0, R - r2, r2, :green)
   circle2(x31, y31, r3)
   circle2(x32, y32, r3)
   circle2(x33, y33, r3)
   segment(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2))
   segment(-x01, sqrt(R^2 - x01^2), -x02, -sqrt(R^2 - x02^2))
   if more == true
       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 - r2, "乙円:r2,(0,R-r2)", :black, :center, :bottom, delta=delta)
       point(0, r1 - R, "甲円:r1,(0,r1-R)", :blue, :center, delta=-delta)
       point(0, R, "R", :green, :center, :bottom, delta=delta)
       point(x01, sqrt(R^2 - x01^2), "(x01,y01)", :black, :left, :bottom, delta=delta)
       point(x02, -sqrt(R^2 - x02^2), "(x02,y02)", :black, :left, delta=-delta, deltax=-delta)
       point(x31, y31, "丙円:r3\n(x31,y31)", :red, :center, :bottom, delta=delta/2)
       point(x32, y32, "丙円:r3\n(x32,y32)", :red, :center, :bottom, delta=delta/2)
       point(x33, y33, "丙円:r3\n(x33,y33)", :red, :center, :bottom, delta=delta/2)
   end
end;

draw(22/2, 11/2, true)