算額（その1145）

算額（その1145）

一〇一 大宮市高鼻町 氷川神社 明治31年(1898)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円2個，楕円2個，長方形，斜線2本

長方形の中に2本の斜線を引き，大円 2 個と楕円 2 個を容れる。楕円の長径が 4 寸，大円の直径が 8 寸のとき，楕円の短径はいかほどか。

長方形の長辺の長さを 2c, （短辺の長さは 2r）
斜線と長方形の長辺との交点の座標を (d, 0)
大円の半径と中心座標を r, (c - r, r)
楕円の長半径，短半径，中心座標を a, b, (0, b), (0, 2b + a)
楕円と斜線の接点座標を (x1, y1), (x2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms a::positive, b::positive, r::positivr,
     c::positive, d::positive,
     x1::positive, y1::positive,
     x2::positive, y2::positive
@syms a, b, r, c, d, x1, y1, x2, y2
slope = -2r/d
eq1 = x1^2/a^2 + (y1 - b)^2/b^2 - 1
eq2 = -b^2*x1/(a^2*(y1 - b)) - slope
eq3 = y1/(d - x1) + slope

eq4 = x2^2/b^2 + (y2 - a - 2b)^2/a^2 - 1
eq5 = -a^2*x2/(b^2*(y2 - a - 2b)) - slope
eq6 = y2/(d - x2) + slope

eq7 = dist2(0, 2r, d, 0, c - r, r, r);

まず，eq1, eq2, eq3 を解き，b, x1, x2 を求める。

res1 = solve([eq1, eq2, eq3], (b, x1, y1))[1]

   (-a^2*r/d^2 + r, 2*a^2*d/(a^2 + d^2), 2*r*(-a^2 + d^2)/(a^2 + d^2))

eq4, eq5 に b を代入し，eq6, eq7 とともに連立方程式を解き，x2, y2, c, d  を求める。

eq4 = eq4(b => res1[1]);
eq5 = eq5(b => res1[1]);
res2 = solve([eq4, eq5, eq6, eq7], (x2, y2, c, d))[8]  # 8 of 8

   (2*sqrt(a^2*r^2/(a^2 - a*r + r^2))*(a^4 - 3*a^3*r + 4*a^2*r^2 - 3*a*r^3 + r^4)/(r^2*(2*a^2 - 3*a*r + 2*r^2)), 2*a*(-2*a^3 + 6*a^2*r - 6*a*r^2 + 3*r^3)/(r*(2*a^2 - 3*a*r + 2*r^2)), r*sqrt((5*a^2 - 4*a*r + 4*r^2)/(a^2 - a*r + r^2))/2 + r + sqrt(a^2*r^2/(a^2 - a*r + r^2))/2, sqrt(a^2*r^2/(a^2 - a*r + r^2)))

短半径は最初の連立方程式の解で -a^2*r/d^2 + r となっており，まだ未知数 d を含んでいる。
d は次の連立方程式の解で sqrt(a^2*r^2/(a^2 - a*r + r^2)) であることがわかる。
そこで，d を短半径の式に代入すれば a, r のみを含む式 a*(r - a)/r = (a*r - a^2)/r になる。これは「術」と同じである。
長半径が a = 4/2，大円の半径が r = 8/2 のとき，短半径は　a*(-a + r)/r = 1 である（短径は 2）。

res1[1](d => res2[4]) |> simplify |> println

   a*(-a + r)/r

術は以下のようになっている。

長径 = 4
大径 = 8
(長径*大径 - 長径^2)/大径

   2.0

その他のパラメータは以下の通りである。

  a = 2;  r = 4;  b = 1;  x1 = 1.97949;  y1 = 1.14286;  x2 = 0.866025;  y2 = 5;  c = 9.31803;  d = 2.3094

function draw(a, r, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x2, y2, c, d) = (2*sqrt(a^2*r^2/(a^2 - a*r + r^2))*(a^4 - 3*a^3*r + 4*a^2*r^2 - 3*a*r^3 + r^4)/(r^2*(2*a^2 - 3*a*r + 2*r^2)), 2*a*(-2*a^3 + 6*a^2*r - 6*a*r^2 + 3*r^3)/(r*(2*a^2 - 3*a*r + 2*r^2)), r*sqrt((5*a^2 - 4*a*r + 4*r^2)/(a^2 - a*r + r^2))/2 + r + sqrt(a^2*r^2/(a^2 - a*r + r^2))/2, sqrt(a^2*r^2/(a^2 - a*r + r^2)))
   (b, x1, y1) = (-a^2*r/d^2 + r, 2*a^2*d/(a^2 + d^2), 2*r*(-a^2 + d^2)/(a^2 + d^2))
   @printf("楕円の長径が %g，大円の直径が %g のとき，楕円の短径は %g である。\n", 2a, 2r, 2b)
   @printf("a = %g;  r = %g;  b = %g;  x1 = %g;  y1 = %g;  x2 = %g;  y2 = %g;  c = %g;  d = %g\n", a, r, b, x1, y1, x2, y2, c, d)
   plot([c, c, -c, -c, c], [0, 2r, 2r, 0, 0], color=:blue, lw=0.5)
   plot!([-d, 0, d], [0, 2r, 0], color=:magenta, lw=0.5)
   circle2(c - r, r, r)
   ellipse(0, b, a, b, color=:green)
   ellipse(0, a + 2b, b, a, color=:green)
   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(c, 0, " c", :blue, :left, :bottom, delta=delta/2)
       point(d, 0, " d", :magenta, :left, :bottom, delta=delta/2)
       point(c - r, r, "大円:r,(c-r,r)", :red, :center, delta=-delta/2)
       point(0, 2b, "2b", :green, :center, delta=-2delta)
       point(0, 2b + 2a, "2b+2a", :black, :center, :bottom, delta=delta)
   end
end;