算額（その1370）

算額（その1370）

六十 群馬県佐波郡玉村町飯塚 光琳寺 天保10年(1837)
六十一 群馬県高崎市石原町 清水寺 天保10年(1837)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：累円
#Julia, #SymPy, #算額, #和算

等脚台形の中に半円，全円，次円，三円，四円，五円を容れる。上底（大頭）が 90 寸のとき，各円の直径はいかほどか。

上底，下底を 2b, 2a
半円の半径と中心座標を r0, (0, r0); r0 = 2r1 = b/2
全円の半径と中心座標を r1, (0, r1)
次円の半径と中心座標を r2, (x21, r2), (x22, r0 - r2)
三円の半径と中心座標を r3, (x3, r3)
四円の半径と中心座標を r4, (x4, r4)
五円の半径と中心座標を r5, (x5, r5)
とおき，以下の連立方程式を解く。
四円のパラメータまでは一度に求めることができるが，SymPy の性能上，五円のパラメータも一緒に求めることができないので，分けて解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, r0::positive, r1::positive,
     r2::positive, x21::positive, x22::positive,
     r3::positive, x3::positive,
     r4::positive, x4::positive,
     r5::positive, x5::positive;
r0 = 2r1
eq1 = x21^2 + (r0 - r2)^2 - (r0 + r2)^2
eq2 = x22^2 + r2^2 - (r0 - r2)^2
eq3 = x22^2 + (r0 - r2 - r1)^2 - (r1 + r2)^2
eq4 = dist2(r0, r0, a, 0, x21, r2, r2)
eq5 = x3^2 + (r0 - r3)^2 - (r0 + r3)^2
eq6 = (x21 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq7 = x4^2 + (r0 - r4)^2 - (r0 + r4)^2
eq8 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (r2, x21, x22, r3, x3, r4, x4, a))[4]  # 4 of 6

   (r1/2, 2*r1, sqrt(2)*r1, 2*r1/9, 4*r1/3, r1/8, r1, r1*(sqrt(2)/2 + 2))

@syms r5::positive, x5::positive
(r2, x21, x22, r3, x3, r4, x4, a) = (r1/2, 2*r1, sqrt(Sym(2))*r1, 2*r1/9, 4*r1/3, r1/8, r1, r1*(sqrt(Sym(2))/2 + 2))
eq9 = x5^2 + (r0 - r5)^2 - (r0 + r5)^2
eq10 = (x4 - x5)^2 + (r4 - r5)^2 - (r4 + r5)^2
#res2 = solve([eq7, eq8], (r4, x4,))
res2 = solve([eq9, eq10], (r5, x5))[1]  # 1 of 2

   (2*r1/25, 4*r1/5)

次円（二円），三円，四円，五円の半径は 全円の半径の 1/2, 2/9, 1/8, 2/25 倍である。

n 番目の円の半径は，全円の半径の 2/n^2 である。

上底が 90 のとき，次円の直径 = 22.5，三円の直径 = 10，四円の直径 = 5.625，五円の直径 = 3.6 である。

全てのパラメータは以下のとおりである。

   r2 = 11.25;  x21 = 45;  x22 = 31.8198;  r3 = 5;  x3 = 30;  r4 = 2.8125;  x4 = 22.5;  r5 = 1.8;  x5 = 18;  a = 60.9099

function draw(b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = b/2
   r1 = r0/2
   (r2, x21, x22, r3, x3, r4, x4, a) = (r1/2, 2*r1, sqrt(2)*r1, 2*r1/9, 4*r1/3, r1/8, r1, r1*(sqrt(2)/2 + 2))
   (r5, x5) = (2*r1/25, 4*r1/5)
   @printf("上底が %g のとき，次円の直径 = %g，三円の直径 = %g，四円の直径 = %g，五円の直径 = %g である。\n", b, 2r2, 2r3, 2r4, 2r5)
   @printf("r2 = %g;  x21 = %g;  x22 = %g;  r3 = %g;  x3 = %g;  r4 = %g;  x4 = %g;  r5 = %g;  x5 = %g;  a = %g\n", r2, x21, x22, r3, x3, r4, x4, r5, x5, a)
   plot([a, r0, -r0, -a, a], [0, r0, r0, 0, 0], color=:magenta, lw=0.5)
   circle(0, r0, r0, beginangle=180, endangle=360)
   circle(0, r1, r1, :blue)
   circle2(x21, r2, r2, :green)
   circle2(x22, r0 - r2, r2, :green)
   circle2(x3, r3, r3, :orange)
   circle2(x4, r4, r4, :brown)
   circle2(x5, r5, r5, :tomato)
   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(b/2, b/2, "(b/2,b/2)", :magenta, :right, :bottom, delta=delta)
       point(a, 0, "a", :magenta, :center, delta=-4delta)
       point(0, r1, "全円:r1,(0,r1)", :blue, :center, delta=-delta)
       point(x21, r2, "次円:r2\n(x21,r2)", :green, :center, :bottom, delta=delta)
       point(x22, r0 - r2, "次円:r2\n(x22,r0-r2)", :green, :center, :bottom, delta=delta)
       point(x3, r3, "三円:r3,(x3,r3)", :orange, :left, delta=-18delta)
       point(x4, r4, "四円:r4,(x4,r4)", :brown, :left, delta=-22delta)
       point(x5, r5, "五円:r5,(x5,r5)", :tomato, :left, delta=-28delta)
   end
end;

draw(90, true)