算額（その400）

算額（その400）

埼玉県川島町下小見野 光西寺観音堂 明治25年(1892)
https://yamabukiwasan.sakura.ne.jp/kousaiji.pdf

鈎股弦の中に菱形2個，甲円，乙円，丙円，全円を入れる。甲円，乙円，丙円の直径がそれぞれ 12.5, 7.5, 12 寸のとき，全円の直径はいかほどか。

鈎の上に ya，股の上に xa, xb をとる。その他 yb，yc を定める。
甲円の半径と中心座標を r1, (xb + r1, r1)
乙円の半径と中心座標を r2, (xa + r2, r2)
丙円の半径と中心座標を r3, (r3, r3)

以下の連立方程式を立て解を求める。

include("julia-source.txt");

using SymPy

@syms xa::positive, xb::positive, ya::positive, yb::positive, yc::positive,
     鈎::positive, 股::positive, 弦::positive,
     r0::positive, r1::positive, r2::positive, r3::positive;

(r1, r2, r3) = (125, 75, 120) .// 20
eq1 = ya/xa - 鈎/股
eq2 = yb/(xb - xa) - 鈎/股
eq3 = yc/(股 -xb) - 鈎/股
eq4 = ya + xa - sqrt(ya^2 + xa^2) - 2r3
eq5 = yb + (xb - xa) - sqrt(yb^2 + (xb - xa)^2) - 2r2
eq6 = yc + (股 - xb) - sqrt(yc^2 + (股 - xb)^2) - 2r1
eq7 = ya + yb + yc - 鈎
eq8 = 鈎 + 股 - 弦 - 2r0
eq9 = 鈎^2 + 股^2 - 弦^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9])

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)
   (xa, xb, ya, yb, yc, 鈎, 股, 弦, r0) = u
   return [
       -鈎/股 + ya/xa,  # eq1
       yb/(-xa + xb) - 鈎/股,  # eq2
       yc/(-xb + 股) - 鈎/股,  # eq3
       xa + ya - sqrt(xa^2 + ya^2) - 12,  # eq4
       -xa + xb + yb - sqrt(yb^2 + (-xa + xb)^2) - 15/2,  # eq5
       -xb + yc + 股 - sqrt(yc^2 + (-xb + 股)^2) - 25/2,  # eq6
       ya + yb + yc - 鈎,  # eq7
       -2*r0 - 弦 + 股 + 鈎,  # eq8
       -弦^2 + 股^2 + 鈎^2,  # eq9
   ]
end;

iniv = [big"48.7", 79.1, 14.0, 8.7, 14.5, 37.2, 129.9, 135.1, 16.0]
res = nls(H, ini=iniv);
println(res);

  (BigFloat[50.14137786721619150424442148054351701570094903819146137553978665178239848983374, 81.47973903422631119439718490588321515051404218706112472703481843278690018275391, 13.88771365970725570233435600656382262425686180304594397465088775657750893885396, 8.679821037317034813958972504102389140160538626903714984058309678348674823548669, 14.46636839552839135659828750683731523360089771150619163995304430858705258296429, 37.03390309255268187289161601750352699801829814145585059866224174351323634536706, 133.7103409792431773446517906147827120418691974351772302874933122640360306862704, 138.7442440717958592175434066322862390398874955766330808838267099554706218187839, 16.00000000000000000000000000000000000000000000000000000116442202603932260642599], true)

   xa = 50.1414;  xb = 81.4797;  ya = 13.8877;  yb = 8.67982;  yc = 14.4664
   鈎 = 37.0339;  股 = 133.71;  弦 = 138.744;  r0 = 16
   全円の直径 = 32;  甲円の直径 = 12.5;  乙円の直径 = 7.5;  丙円の直径 = 12

全円の半径は 16 寸（直径 32 寸）である。

using Plots

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = (125, 75, 120) .// 20
   (xa, xb, ya, yb, yc, 鈎, 股, 弦, r0) = res[1]
   @printf("xa = %g;  xb = %g;  ya = %g;  yb = %g;  yc = %g\n鈎 = %g;  股 = %g;  弦 = %g;  r0 = %g\n", xa, xb, ya, yb, yc, 鈎, 股, 弦, r0)
   @printf("全円の直径 = %g;  甲円の直径 = %g;  乙円の直径 = %g;  丙円の直径 = %g\n", 2r0, 2r1, 2r2, 2r3)
   x1 = xb + r1
   x2 = xa + r2
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(x1, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(r3, r3, r3, :green)
   circle(r0, r0, r0, :orange)
   segment(0, ya, xa, 0)
   segment(xa, yb, xb, 0)
   segment(xa, 0, xa, yb + yc)
   segment(xb, 0, xb, yc)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, 鈎, " 鈎=ya+yb+yc", :black, :left, :bottom)
       point(股, 0, " 股", :black, :left, :bottom, delta=2delta)
       point(xa, 0, "xa", :black, :center, :top, delta=-2delta)
       point(xb, 0, "xb", :black, :center, :top, delta=-2delta)
       point(0, ya, "ya ", :black, :right, :vcenter)
       point(xa, yb, "(xa,yb)", :black, :right, :vcenter)
       point(xa, yb + yc, "(xa,yb+yc)", :black, :left, :bottom, delta=2delta)
       point(xb, yc, "(xb,yc)", :black, :left, :bottom, delta=2delta)
       point(xb + r1, r1, " 甲円:r1,(xb+r1,r1)", :red, :left, :bottom)
       point(xa + r2, r2, " 乙円:r2,(xa+r2,r2)", :blue, :left, :top, delta=2delta)
       point(r3, r3, " 丙円:r3,(r3,r3)", :green, :left, :top, delta=2delta)
       point(r0, r0, " 全円:r0,(r0,r0)", :orange, :left, :vcenter)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       plot!(xlims=(-8, 140), ylims=(-6, 42))
   else
       plot!(showaxis=false)
   end
end;