算額（その1622）
～落書き帳「○△□」～ 392．○△□の新算額（その５）
http://streetwasan.web.fc2.com/math18.2.1.html
キーワード：円4個，外円，円弧
#Julia, #Julia, #SymPy, #算額, #和算, #数学
外円の中に，甲円 1 個，乙円 1 個，丙円 2 個を容れる。外円の直径が与えられ，甲円，乙円，丙円の中心を結ぶと直角三角形になるとき，丙円の直径はいかほどか。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (0, r2 - R)
丙円の半径と中心座標を r3, (x3, y3)
外円と乙円の交点座標を (x0, y0)
とおき，以下の連立方程式を解く。
注：「甲円，乙円，丙円の中心を結ぶと直角三角形になる」という条件を eq4 で記述する。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms R::integer, r1::positive, r2::positive,
r3::positive, x3::positive, y3::negative,
x0::positive, y0::positive
eq1 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2
eq2 = x3^2 + (y3 + R)^2 - (r2 + r3)^2
eq3 = x3^2 + y3^2 - (R - r3)^2
eq4 = (r1 + r2)^2 - ((r1 + r3)^2 + (r2 + r3)^2)
eq5 = 2r1 + r2 - 2R
eq6 = x0^2 + (y0 + R)^2 - r2^2
eq7 = x0^2 + y0^2 - R^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, r2, r3, x3, y3, x0, y0))[1]; # 1 of 3
@syms d, t
# r1
ans_r1 = res[1](629 + 48√Sym(177) => t) |> simplify |> println
R*(-10*sqrt(177)*t^(2/3) + 109*t^(2/3) - 46*sqrt(177)*t^(1/3) + 713*t^(1/3) + 3703)/3174
# r2
res[2](629 + 48√Sym(177) => t) |> simplify |> println
R*(-109*t^(2/3) + 10*sqrt(177)*t^(2/3) - 713*t^(1/3) + 46*sqrt(177)*t^(1/3) - 529)/1587
# r3
res[3](629 + 48√Sym(177) => t) |> simplify |> println
R*(-48*sqrt(177)*t^(2/3) + 629*t^(2/3) + 529*t^(1/3) - 3703)/3174
t = 629 + 48√177 とおいたとき，
外円の半径が R のとき，丙円の半径 r3 は R*(-48√177t^(2/3) + 629t^(2/3) + 529t^(1/3) - 3703)/3174 である。
たとえば，外円の直径が 1 のとき，丙円の直径は 0.282881192474901 である。
res[3](R => 1).evalf() |> println
0.282881192474901
# x3
res[4](629 + 48√Sym(177) => t) |> simplify |> println
2*sqrt(-3447*t^(2/3) + 258*sqrt(177)*t^(2/3) - 138*sqrt(177)*t^(1/3) - 1035*t^(1/3) + 33327)*Abs(R)/69
# y3
res[5](629 + 48√Sym(177) => t) |> simplify |> println
R*(-88*sqrt(177)*t^(2/3) + 1065*t^(2/3) - 184*sqrt(177)*t^(1/3) + 3381*t^(1/3) + 1587)/3174
# x0
res[6](629 + 48√Sym(177) => t) |> simplify |> println
sqrt(-1248*sqrt(177)*t^(2/3) + 15825*t^(2/3) - 1104*sqrt(177)*t^(1/3) + 28221*t^(1/3) - 46023)*Abs(R)/138
# y0
apart(res[7], d)(629 +48√Sym(177) => t) |> simplify |> println
R*(t^(1/3)*(t^(1/3) - 13) - 23)/(6*t^(1/3))
function draw(R, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
t = 629 +48√177
t13 = t^(1/3)
t23 = t13^2
r1 = R*(-10√177t23 + 109t23 - 46√177*t13 + 713t13 + 3703)/3174
r2 = R*(-109t23 + 10√177t23 - 713t13 + 46√177*t13 - 529)/1587
r3 = R*(-48√177t23 + 629t23 + 529t13 - 3703)/3174
x3 = 2R*sqrt(-3447t23 + 258√177*t23 - 138√177t13 - 1035t13 + 33327)/69
y3 = R*(-88√177t23 + 1065t23 - 184√177t13 + 3381t13 + 1587)/3174
x0 = R*sqrt(-1248√177t23 + 15825t23 - 1104√177t13 + 28221t13 - 46023)/138
y0 = R*(t13*(t13 - 13) - 23)/(6t13)
@printf("外円の直径が %g のとき，丙円の直径は %g である。\n", 2R, 2r3)
@printf("R = %g; r1 = %g; r2 = %g; r3 = %g; x3 = %g; y3 = %g\n", R, r1, r2, r3, x3, y3)
θ = atand(y0 + R, x0)
println(θ)
plot()
circle(0, 0, R, :green)
circle(0, R - r1, r1)
circle(0, -R, r2, :blue, beginangle=θ, endangle=180 - θ)
circle2(x3, y3, r3, :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(0, R, "R", :green, :center, :bottom, delta=delta/2)
point(0, R - r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta)
point(0, -R, "乙円:r2,(0,-R)", :blue, :center, :bottom, delta=2delta)
point(x3, y3, "丙円:r3\n(x3,y3)", :orange, :center, delta=-delta)
end
end;
draw(1/2, true)
