算額（その674）
『不朽算法』上巻第六問 [日下貞八郎誠嗣辺・安島万蔵直円遺稿]
米山忠興: 側円術(IV) ---『安子遺稿側円解二条第一』の現代語訳 ---，東洋大学紀要. 自然科学篇，Vol. 52, pp. 97-116, 2008-03
http://id.nii.ac.jp/1060/00002535/
直角三角形（鈎股弦）内に円と楕円が入っている。円は直角三角形の三辺に接しており，楕円と外接している。楕円の長径は底辺（股）に平行で，直角三角形の斜辺（弦）に接している。
直角三角形の辺の長さと楕円の短径を与えて，長径を求めよ。
直角三角形の直角を挟む二辺のうち，短い方を「鈎」，長い方を「股」とする（斜辺「弦」は sqrt(鈎^2 + 股^2) である）。
円の半径 r は「(鈎 + 股 - 弦)/2」，中心座標は (r, r)
楕円の長径と短径を 2a, 2b，中心座標を (x0, y0) とする。
楕円と斜辺，楕円と円の接点座標を (x1, y1), (x2, y2)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms 鈎, 股, 弦, a, b, r, x, x1, y1, x2, y2, x0, y0, slope0
(鈎, 股, b) = (2240, 2808, 455//2)
弦 = 3592 # sqrt(鈎^2 + 股^2)
r = (鈎 + 股 - 弦)//2
y0 = b
長径，短径が a b，中心座標が (x0, y0) の楕円
(x1, y1) が楕円上にある
eq1 = (x1 - x0)^2/a^2 + (y1 - y0)^2/b^2 - 1
(x1, y1) における楕円の接線の傾き
slope0 = -鈎//股 # 目的とする接線。符号に注意
eq2 = b^2*(x1 - x0)/(a^2*(y0 - y1)) - slope0 # slope
もう一個の条件
eq3 = y1/(股 - x1) + slope0
楕円と円が外接する
eq4 = (x2 - x0)^2/a^2 + (y2 - y0)^2/b^2 - 1
slope2 = (x2 - r)/(r - y2)
eq5 = b^2*(x2 - x0)/(a^2*(y0 - y2)) - slope2
eq6 = (x2 - r)^2 + (y2 - r)^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 v, r.f_converged
end;
function H(u)
(x1, y1, x0, x2, y2, a) = u
return [
4*(y1 - 455/2)^2/207025 - 1 + (-x0 + x1)^2/a^2, # eq1
280/351 + (-207025*x0/4 + 207025*x1/4)/(a^2*(455/2 - y1)), # eq2
y1/(2808 - x1) - 280/351, # eq3
4*(y2 - 455/2)^2/207025 - 1 + (-x0 + x2)^2/a^2, # eq4
-(x2 - 728)/(728 - y2) + (-207025*x0/4 + 207025*x2/4)/(a^2*(455/2 - y2)), # eq5
(x2 - 728)^2 + (y2 - 728)^2 - 529984, # eq6
]
end;
iniv = BigFloat[2398, 327, 1870, 1300, 230, 580]
res = nls(H, ini=iniv)
(BigFloat[2397.842696629213493816709019056387960138168817704006327516310149222101255756478, 327.1910112359550540245449460804349431777480338181361396304024428466875180291161, 1872.000000000000004115470053719782306992126301163385447666330986804201857308351, 1310.399999999999999733344750238321544467653247568361438815996098107649959134677, 291.1999999999999996444596669844287261829629732073189591407411915504426299859508, 585.0000000000000042869479726247732354323812232955671058370412692641114894978273], true)
楕円の長径(a) = 585（差渡し径(2a) = 1170）である。
その他のパラメータは以下の通り。
鈎 = 2240; 股 = 2808; b = 227.5
x1 = 2397.84; y1 = 327.191; x0 = 1872; y0 = 227.5
x2 = 1310.4; y2 = 291.2; r = 728
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(鈎, 股, b) = (2240, 2808, 455//2)
弦 = sqrt(鈎^2 + 股^2)
r = (鈎 + 股 - 弦)/2
y0 = b
(x1, y1, x0, x2, y2, a) = res[1]
@printf("鈎 = %g; 股 = %g; b = %g\n", 鈎, 股, b)
@printf("楕円の長径(a) = %g （差渡し径(2a) = %g\n", a, 2a)
@printf("x1 = %g; y1 = %g; x0 = %g; y0 = %g\n", x1, y1, x0, y0)
@printf("x2 = %g; y2 = %g; r = %g\n", x2, y2, r)
plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
ellipse(x0, b, a, b, color=:red)
circle(r, r, r, :blue)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
point(x0, y0, "(x0,y0)")
point(r, r, "(r,r)")
point(0, 鈎, " 鈎", :black, :left, :bottom)
point(股, 0, " 股", :black, :left, :bottom, delta=delta/2)
point(x1, y1, "(x1,y1) ", :red, :right)
point(x2, y2, "(x2,y2) ", :red, :right)
end
end;
