算額（その1398）

算額（その1398）

十 胆沢町若柳市野々 個人宅 安政2年(1855)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円4個，楕円2個
#Julia, #SymPy, #算額, #和算

2 個の楕円が外接しており，それぞれの楕円の中に等円が 1 個ずつ内接し，2 個の等円が楕円に外接している。楕円の長径が 4 寸，短径が 1 寸のとき，等円の直径はいかほどか。

楕円の長径，短径を 2a, 2b
等円の半径と中心座標を r, (x, 0), (0, y)
等円と楕円の接点の座標を (x0, y0)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms a, b, x, y, r
(x0, y0) = (x, y) ./ 2
eq1 = (a - x)^2 - (a^2 - b^2)*(b^2 - r^2)/b^2  # 算法助術公式84
eq2 = x^2 + y^2 - 4r^2  # 等円が外接
eq3 = (x0 - a)^2/a^2 + y0^2/b^2 - 1;  # 接点が楕円の周上にある

res = solve([eq1, eq2, eq3], (r, x, y))[4]  # 4 of 4

   (b^2*sqrt(3*a^2 + b^2)/(a^2 + b^2), 2*a*b^2/(a^2 + b^2), sqrt(8*a^2*b^4/(a^4 + 2*a^2*b^2 + b^4) + 4*b^6/(a^4 + 2*a^2*b^2 + b^4)))

等円の半径は b^2*sqrt(3*a^2 + b^2)/(a^2 + b^2) である。
a = 4/2, b = 1/2 のとき，b^2*sqrt(3*a^2 + b^2)/(a^2 + b^2) = 7/34 である。
直径は 7/17 寸である。

@syms a, b
a = Sym(4)//2
b = Sym(1)//2
b^2*sqrt(3*a^2 + b^2)/(a^2 + b^2) |> println

   7/34

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

  a = 2;  b = 0.5;  r = 0.205882;  x = 0.235294;  y = 0.337915;  x0 = 0.117647;  y0 = 0.168958

function draw(a, b, more, magnify)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho") 
   (r, x, y) = (b^2*sqrt(3*a^2 + b^2)/(a^2 + b^2), 2*a*b^2/(a^2 + b^2), sqrt(8*a^2*b^4/(a^4 + 2*a^2*b^2 + b^4) + 4*b^6/(a^4 + 2*a^2*b^2 + b^4)))
   x0 = x/2
   y0 = y/2
   @printf("楕円の長半径，短半径が %g，%g のとき，等円の直径は %g である。\n", 2a, 2b, 2r)
   @printf("a = %g;  b = %g;  r = %g;  x = %g;  y = %g;  x0 = %g;  y0 = %g\n", a, b, r, x, y, x0, y0)
   plot()
   ellipse(a, 0, a, b, color=:blue)
   ellipse(-a, 0, a, b, color=:blue)
   circle2(x, 0, r)
   circle22(0, y, r)
   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)
       if magnify
           plot!(xlims=(-0.6, 0.6), ylims=(-0.6, 0.6))
           point(x, 0, "等円:r,(x,0)", :red, :center, delta=-delta)
           point(0, y, "等円:r,(0,y)", :red, :center, delta=-delta)
           point(0.35, 0.28, "楕円:a,b,(a,0)", :blue, :left, delta=-delta, mark=false)
           point(x0, y0, "(x0,y0)", :red, :left, delta=-delta)
           point(0.2, 0.45, "中心部の拡大図", :black, mark=false)
       else
           point(a, 0, "楕円:a,b,(a,0)", :blue, :center, delta=-6delta)
           point(x, 0, "", :red)
           point(0, y, "", :red)
       end
   end  
end;

draw(4/2, 1/2, true, true)