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