算額（その1434）
福島県伊達市梁川町二野袋庭渡 庭渡神社 明治27年(1894)
街角の数学 Street Wasan ～落書き帳「○△□」～ 40．正方形と四分円（5）
http://streetwasan.web.fc2.com/math15.6.22.html
キーワード：円2個，四分円2個，長方形，斜線2本
#Julia, #SymPy, #算額, #和算
長方形の中に 2 個の四分円を容れ，長方形の頂点から接戦を引く。区画された領域に乾円，坤円を 1 個ずつ容れる。坤円の半径を与えて，乾円の半径を求めよ。
長方形の長辺，短辺を 2a, b
接線と長方形の長辺の交点座標を (c, 0)
四分円の半径と中心座標を b, (a, b), (-a, b)
乾円の半径と中心座標を r1, (0, b - r1)
坤円の半径と中心座標を r2, (0, r2)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, b::positive, c::positive, r1::positive, r2::positive
eq1 = a^2 + r1^2 - (r1 + b)^2
eq2 = dist2(-a, b, c, 0, 0, b - r1, r1)
eq3 = dist2(-a, b, c, 0, 0, r2, r2)
eq4 = dist2(-a, b, c, 0, a, b, b);
res = solve([eq1, eq2, eq3, eq4], (a, b, c, r1))[1] # 1 of 2
(-r2*sqrt(11 + 8*sqrt(2))*sqrt(56*sqrt(2) + 81)/6 + r2*sqrt(11 + 8*sqrt(2))/3 + sqrt(2)*r2*sqrt(11 + 8*sqrt(2))*sqrt(56*sqrt(2) + 81)/6, 3*sqrt(2)*r2/2 + 5*r2/2, r2*sqrt(11 + 8*sqrt(2)), -r2*sqrt(112*sqrt(2) + 162)/9 - r2*sqrt(56*sqrt(2) + 81)/9 - 5*r2/18 + sqrt(2)*r2/18 + sqrt(2)*r2*sqrt(56*sqrt(2) + 81)/9 + sqrt(2)*r2*sqrt(112*sqrt(2) + 162)/9)
# a
res[1] |> simplify |> println
res[1](r2 => 1).evalf() |> println
r2*(-sqrt(1787 + 1264*sqrt(2)) + 2*sqrt(11 + 8*sqrt(2)) + sqrt(3574 + 2528*sqrt(2)))/6
5.70205716976694
# b
res[2] |> simplify |> println
res[2](r2 => 1).evalf() |> println
r2*(3*sqrt(2) + 5)/2
4.62132034355964
# c
res[3] |> println
res[3](r2 => 1).evalf() |> println
r2*sqrt(11 + 8*sqrt(2))
4.72373882628843
# r1
res[4] |> simplify |> sympy.sqrtdenest |> simplify |> println
res[4](r2 => 1).evalf() |> println
r2*(1 + sqrt(2))/2
1.20710678118655
乾円の半径 r1 は，坤円の半径 r2 の (1 + √2)/2 倍である。かなりきれいな倍数だと思う。
坤円の半径 r2 が 1 のとき，乾円の半径 r1 は (1 + √2)/2 = 1.2071067811865475 である。
全てのパラメータは以下のとおりである。
r2 = 1; r1 = 1.20711; a = 5.70206; b = 4.62132; c = 4.72374
function draw(r2, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
a = r2*(-sqrt(1787 + 1264√2) + 2*sqrt(11 + 8√2) + sqrt(3574 + 2528√2))/6
b = r2*(3√2 + 5)/2
c = r2*sqrt(11 + 8√2)
r1 = r2*(1 + √2)/2
@printf("r2 = %g; r1 = %g; a = %g; b = %g; c = %g\n", r2, r1, a, b, c)
plot([a, a, -a, -a, a], [0, b, b, 0, 0], color=:green, lw=0.5)
circle(a, b, b, :magenta, beginangle=180, endangle=270)
circle(-a, b, b, :magenta, beginangle=270, endangle=360)
circle(0, b - r1, r1)
circle(0, r2, r2, :blue)
segment(-a, b, c, 0, :orange)
segment(a, b, -c, 0, :orange)
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(a, 0, "a", :green, :center, delta=-delta)
point(c, 0, "c", :black, :center, delta=-delta)
point(a - b, b, "(a-b,b)", :magenta, :center, :bottom, delta=delta)
point(0, b - r1, "乾円:r1\n(0,b-r1)", :red, :center, delta=-delta)
point(0, r2, "坤円:r2\n(0,r2)", :blue, :center, delta=-delta)
ylims!(-8delta, b + 2delta)
end
end;
draw(1, true)
