算額（その1390）
十七 大里郡岡部村岡 稲荷社
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円8個，直線上，斜線2本
#Julia, #SymPy, #算額, #和算
直線と斜線 2 本を引き，円 8 個を載せる。丁円の直径が 33 寸のとき，丙円の直径はいかほどか。
甲円の半径と中心座標を r1, (-r1, r1), (r1, r1)
乙円の半径と中心座標を r2, (0, 2r4 + r2)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (0, r4), (2r1, r4)
とおく。
甲円と丁円が接することから，関係式 eq1 が成り立つ。
それを解けば，r1 は r4 の 4 倍であることが簡単にわかる。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r1::positive, r2::positive,
r3::positive, x3::positive,
r4::positive, a::positive
eq1 = r1 - 2sqrt(r1*r4)
solve(eq1, r1)[1] |> println
4*r4
甲円と同じ大きさの補助円（灰色）を描くと，左側の甲円の中心 (-r1, r1) と右側の丁円の中心 (2r1, r4) と補助円と x 軸の接点 (3r1, 0) は一直線上にある。
直角三角形 ⊿ABC と ⊿DBE の相似関係から sqrt((r1 + a)^2 + 16r1^2) = 4a が成り立つので a = 68r4/15 がわかる。
eq2 = sqrt((r1 + a)^2 + 16r1^2) - 4a
solve([eq1, eq2], (a, r1))[1] |> println
(68*r4/15, 4*r4)
乙円と丙円，丙円と丁円が外接するので，以下の関係式が成り立つ。
eq3 = x3^2 + (2r4 + r2 - r3)^2 - (r2 + r3)^2
eq4 = x3^2 + (r3 - r4)^2 - (r3 + r4)^2;
左側の甲円の中心，乙円の中心，右側の丁円の中心 から (-r1, r1 + a) と (3r1, 0) を結ぶ直線までの距離はそれぞれ r1, r2, r4 である。
この中から，これまでの eq1, eq2, eq3 ,eq4 と独立な関係式を用いる。
eq5 = dist2(-r1, r1 + a, 3r1, 0, 0, 2r4 + r2, r2);
# eq6 = dist2(-r1, r1 + a, 3r1, 0, -r1, r1, r1)
# eq7 = dist2(-r1, r1 + a, 3r1, 0, 2r1, 0, r4);
以上の 5 本の方程式を連立させて解き，a, r1, r2, r3, x3 を求める。
solve([eq1, eq2, eq3, eq4, eq5], (a, r1, r2, r3, x3))[1]
(68*r4/15, 4*r4, 33*r4/16, 49*r4/33, 14*sqrt(33)*r4/33)
丙円の半径 r3 は，丁円 の半径の 49/33 倍である。
丁円の直径が 33 寸のとき，丙円の直径は 49 寸である。
全てのパラメータは以下の通りである。
r4 = 16.5; a = 74.8; r1 = 66; r2 = 34.0312; r3 = 24.5; x3 = 40.2119
function draw(r4, more)
pyplot(size=(800, 400), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, r1, r2, r3, x3) = r4 .* (68/15, 4, 33/16, 49/33, 14/√33)
@printf("r4 = %g; a = %g; r1 = %g; r2 = %g; r3 = %g; x3 = %g\n", r4, a, r1, r2, r3, x3)
plot()
circle2(r1, r1, r1)
circle(3r1, r1, r1, :gray60)
circle(0, r4, r4, :blue)
circle2(2r1, r4, r4, :blue)
circle(0, 2r4 + r2, r2, :orange)
circle2(x3, r3, r3, :green)
segment(-r1, r1 + a, 3r1, 0, :magenta)
segment(r1, r1 + a, -3r1, 0, :magenta)
segment(-r1, r1, 3r1, 0, :gray60)
segment(-r1, 0, -r1, a + r1)
if more == true
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(-r1, r1 + a, "B:(-r1,r1+a) ", :black, :right, :vcenter)
point(-r1, r1, "C:(-r1,r1) ", :red, :right, :vcenter)
point(3r1, 0, "E:3r1", :red, :center, delta=-3delta)
point(2r1, 0, "2r1", :red, :center, delta=-3delta)
point(r1, r1, "甲円:r1,(r1,r1)", :red, :left, :bottom, delta=2delta)
point(0, 2r4 + r2, "乙円:r2\n(0,2r4+r2)", :orange, :left, :vcenter, deltax=delta)
point(0, r4, "丁円:r4\n(0,r4)", :blue, :center, :vcenter)
point(2r1, r4, "丁円:r4\n(2r1,r4)", :blue, :center, :vcenter)
point(x3, r3, " 丙円:r3,(x3,r3)", :green, :left, :vcenter)
θ = atand((r1 + a)/4r1)
segment(-r1, r1, -r1 + r1*sind(θ), r1 + r1*cosd(θ), :gray60)
point(-r1 + r1*sind(θ), r1 + r1*cosd(θ), " A", :green, :left, :bottom)
point(-r1, 0, "D:-r1", :green, :center, delta=-3delta)
plot!(xlims=(-5delta - 2r1 - r4, 4r1 + 5delta), ylims=(-10delta, r1 + a + 3delta))
end
end;
draw(33/2, true)
