算額（その1647）
三四 武蔵国埼玉郡下忍村遍照院境内 金毘羅社（神楽堂） 天保11年(1840)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円7個，楕円3個，正方形
#Julia, #Julia, #SymPy, #算額, #和算, #数学
正方形の中に，交差する大楕円 2 個を設け，四隅に丙円 4 個，中央部に小楕円，更にその中に甲円 1 個，乙円 2 個を容れる。乙円とは丙円はそれぞれの楕円の長径端において楕円と 1 点で接する最大の円である（曲率円）。小楕円の短径が 1 寸のとき，大楕円の短径はいかほどか。
大楕円の長半径，短半径，中心座標を a, b, (0, 0)
小楕円の長半径，短半径，中心座標を c, d, (0, 0)
甲円の半径と中心座標を r1, (0, 0); r1 = b
乙円の半径と中心座標を r2, (b - r2, 0)
丙円の半径と中心座標を r4, ((c - r4)/√2, (c - r4)/√2)
大楕円と小楕円の接点座標を A:(x0, y0)
接点から長軸へおろした垂線の脚を B，原点を O としたとき AB = y1, 0B = x1
とおき，以下の連立方程式を解く。
小楕円の中の乙円は曲率円なので，r2 = b^2/a
また，a = b + 2r2 である。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, b::positive, r2::positive,
c::positive, d::positive, r4::positive,
x0::positive, y0::positive, x1::positive, y1::positive
eq01 = r2 - b^2/a
eq02 = a - (b + 2r2)
solve([eq01, eq02], (a, r2))[1]
(2*b, b/2)
a = 2b, r2 = b/2 である。
eq03 = r4 - d^2/c
eq04 = c - (d + 2r4)
solve([eq03, eq04], (c, r4))[1]
(2*d, d/2)
c = 2d, r4 = d/2 である。
大楕円，小楕円は相似である。
y1 = (x0 + y0)/√Sym(2)
eq0 = y0^2 - (4b^2 - x0^2)/4
eq4 = (x0^2 + y0^2) - ((x0 + y0)^2/2 + x1^2)
eq5 = 3√Sym(2)x0*y0/(4y0 + x0) - 3x1/4
eq6 = x1^2 + 4y1^2 - 4d^2
res = solve([eq0, eq4, eq5, eq6], (d, x0, y0, x1))[2] # 2 of 2
(b*sqrt(3*sqrt(41) + 25)/4, b*sqrt(205 - sqrt(41))*(41*sqrt(2) + 13*sqrt(82))/1312, b*sqrt(205 - sqrt(41))*(-3*sqrt(82) + 41*sqrt(2))/1312, b*sqrt(8405 - 41*sqrt(41))/82)
# d
res[1] |> println
res[1].evalf() |> println
b*sqrt(3*sqrt(41) + 25)/4
1.66225322815709*b
大楕円の短半径 d は，小楕円の短半径 b の sqrt(3√41+ 25)/4 = 1.66225322815709 倍」である。
小楕円の短径が 1 寸のとき，大楕円の短径は 1.66225322815709 寸である。
術は「大楕円の短径 d は，小楕円の短径 b の sqrt(√369 + 25)/4 倍」である。
sqrt(√369 + 25)/4 = sqrt(3√41 + 25)/4 なので，同じである。
function draw(b, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(d, x0, y0, x1) = (b*sqrt(3*sqrt(41) + 25)/4,
b*sqrt(205 - sqrt(41))*(41*sqrt(2) + 13*sqrt(82))/1312,
b*sqrt(205 - sqrt(41))*(-3*sqrt(82) + 41*sqrt(2))/1312,
b*sqrt(8405 - 41*sqrt(41))/82)
f = d/b
a = 2b
c = 2d
r1 = b
@printf("factor = %g\n", f)
@printf("a = %g; b = %g; c = %g; d = %g\n", a, b, c, d)
r2 = b^2/a # 曲率円
r3 = r2*f
r4 = d^2/c # 曲率円
s = sqrt(2(c^2 + d^2)) # 算法助術−94
plot(s/2 .*[1, 1, -1, -1, 1], s/2 .*[-1, 1, 1, -1, -1], color=:orange, lw=0.5)
ellipse(0, 0, a, b, color=:red)
ellipse(0, 0, a, b, color=:pink, φ=90)
circle(0, 0, b, :green)
circle2(a - r2, 0, r2, :brown)
ellipse(0, 0, c, d, color=:blue, lw=0.5, φ=45)
ellipse(0, 0, c, d, color=:blue, lw=0.5, φ=135)
ellipse(0, 0, c, d, color=:gray80, lw=0.5, φ=0)
circle(0, 0, d, :lightgreen)
circle4((c - r4)/√2, (c - r4)/√2, r4, :deeppink)
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, 0, "甲円:r1,(0,0)", :green, :center, :bottom, delta=2delta)
point(a - r2, 0, "乙円:r2\n(a-r2,0)", :brown, :center, delta=-delta/2)
point((c - r4)/√2, (c - r4)/√2, "丙円:r4\n((c-r4)/√2,\n(c-r4)/√2)", :magenta, :center, delta=-delta/2)
point(x0, y0, "(x0,y0)", :red, :left, :bottom, delta=delta, deltax=-delta/2)
point(0, b, "b", :red, :center, :bottom, delta=delta/2)
point(0, d, "d", :red, :center, :bottom, delta=delta/2)
point(a, 0, "a", :red, :right, :vcenter, deltax=-delta/2)
point(c, 0, "c", :red, :right, :vcenter, deltax=-delta/2)
segment(0, 0, -s/2, -s/2, :gray80)
point(0, 0, "O", :black, :left, delta=-delta)
segment(0, 0, -x0, y0)
point(-x0, y0, "A", :black, :right, :vcenter, deltax=-delta)
segment(0, 0, -x1/√2, -x1/√2)
point(-x1/√2, -x1/√2, "B", :black, :right, :vcenter, deltax=-delta)
segment(-x1/√2, -x1/√2, -x0, y0)
end
end;
draw(1, true)
