算額（その787）
寛政十年戊午九月(1798) 藤田貞資門人 肥後人吉 西九郎左衛門孟周
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf
正方形内に円弧（弧背）と等円 3 個を入れる。正方形の一辺の長さが 971 寸のとき，等円の直径はいかほどか。
最初に見たとき，正方形の上辺で 2 つの円弧は互いに外接しているかと思ったが，そうではなかった。そこで交わっている。
正方形の一辺の長さを 2a
等円の半径と中心座標を r, (0, r), (a - r, 2a - r)
円弧のもとの円の半径と中心座標を r0, (x0, y0)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r0::positive, x0::positive, y0::positive, a::positive, r::positive
#a = 971//2
eq1 = (x0 - (a - r))^2 + (y0 - (2a - r))^2 - (r0 - r)^2
eq2 = x0^2 + (y0 - r)^2 - (r0 + r)^2
eq3 = (x0 - a)^2 + y0^2 - r0^2
eq4 = x0^2 + (y0 - 2a)^2 - r0^2;
res = solve([eq1, eq2, eq3, eq4], (r, r0, x0, y0))
1-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
(a*(-5*sqrt(5) + sqrt(14 + 10*sqrt(5)) + 7)/4, 5*a*(-34 - 14*sqrt(5) + 5*sqrt(70 + 50*sqrt(5)) + 13*sqrt(14 + 10*sqrt(5)))/(8*(-sqrt(2)*sqrt(7 + 5*sqrt(5)) + 5*sqrt(5) + 13)), a*(-sqrt(14 + 10*sqrt(5)) + 5*sqrt(5) + 17)/8, a*(-sqrt(7/128 + 5*sqrt(5)/128) + 5*sqrt(5)/16 + 29/16))
等円の直径は，正方形の一辺の長さの (sqrt(14 + 10*√5) + 7 - 5√5)/4 倍である。
正方形の一辺の長さが 971 寸のとき，等円の直径は 449.00055949866737 寸である。
数式はそんなに長くないと思うが，術は結構長い。
971*(sqrt(14 + 10*√5) + 7 - 5√5)/4
449.00055949866737
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
a = 971//2
(r, r0, x0, y0) = (a*(-5*sqrt(5) + sqrt(14 + 10*sqrt(5)) + 7)/4, 5*a*(-34 - 14*sqrt(5) + 5*sqrt(70 + 50*sqrt(5)) + 13*sqrt(14 + 10*sqrt(5)))/(8*(-sqrt(2)*sqrt(7 + 5*sqrt(5)) + 5*sqrt(5) + 13)), a*(-sqrt(14 + 10*sqrt(5)) + 5*sqrt(5) + 17)/8, a*(-sqrt(7/128 + 5*sqrt(5)/128) + 5*sqrt(5)/16 + 29/16))
@printf("等円の直径 = %g; r = %g; r0 = %g; x0 = %g\n", 2r, r, r0, x0)
plot([a, a, -a, -a, a], [0, 2a, 2a, 0, 0], color=:blue, lw=0.5)
circle2(a - r, 2a - r, r)
circle(0, r, r)
circle(x0, y0, r0, :green, beginangle=180, endangle= 235)
circle(-x0, y0, r0, :green, beginangle=305, endangle= 360)
if more
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, " 等円:r\n (0,r)", :red, :left, :vcenter)
point(a - r, 2a - r, " 等円:r\n (a-r,2a-r)", :red, :left, :vcenter)
point(0, 2a, " 2a", :blue, :left, :bottom, delta=delta/2)
point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
point(0.3a, a, "弧背:r0,(x0,y0)", :green, :left, mark=false)
end
end;
