算額（その655）

算額（その655）

『蠡管算法　巻上』の第１３問 をまなぶ（１／１）
http://takasakiwasan.web.fc2.com/kaihouhenomiti_1/reikannsannpou_jyou_13.html

「算法直術正解」第35問
http://takasakiwasan.web.fc2.com/kaihouhenomiti_1/pdf/sanpochokujutuseikai_35_2.pdf

団扇の中に菱形と 2 つの等円が入っている。団扇の直径が 10 寸で，菱形の一辺の長さが 6 寸の場合，等円の直径はいかほどか。

後者では 14 ページにわたる解法を示している。

菱形の長径と短径を 2a, 2b
団扇の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x, y)
とする。

菱形の一辺の長さが 6 寸で，右端の頂点が半径 5 の円周上にあることから以下の方程式の eq4 を立てるが，R が与えられているので，b = 18/5 であることが容易に計算できる。

eq1 は等円が団扇に内接すること，eq2 は扇の下部の弧と等円が外接すること，eq3 は等円の中心と菱形の右下の辺までの距離が等円の直径であることを表している。この 3 方程式を解くと r, x, y が正確に求まる。

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

using SymPy
@syms R::positive, a::positive, b::positive,
     r::positive, x::positive, y::negative
R = 10//2
b = 18//5
eq1 = x^2 + y^2 - (R - r)^2
eq2 = x^2 + (y + R)^2 - (2R - 2b + r)^2
eq3 = distance(0, R - 2b, sqrt(R^2 - (R - b)^2), R - b, x, y) - r^2
eq4 = sqrt(6^2 - b^2) - sqrt(R^2 - (R - b)^2);
res = solve([eq1, eq2, eq3], (r, x, y))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (26208/17405, 51408/17405, -6499/3481)

等円の直径は 2*26208/17405 = 3.011548405630566
中心座標は (2.9536340132145935, -1.8669922436081585) である。

2*26208/17405, 51408/17405, -6499/3481

   (3.011548405630566, 2.9536340132145935, -1.8669922436081585)

団扇の円と下部の円弧の交点座標 (ax, ay) は以下の連立方程式を解いて求めることができる。

@syms ax::positive, ay::negative
eq5 = ax^2 + ay^2 - R^2
eq6 = ax^2 + (ay + R)^2 - (2R - 2b)^2
solve([eq5, eq6], (ax, ay))

   1-element Vector{Tuple{Sym, Sym}}:
    (336/125, -527/125)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 10//2
   b = 18//5
   (r, x, y) = res[1]
   a = sqrt(6^2 - b^2)
   @printf("等円の直径 = %g;  r = %g;  x = %g;  y = %g;  a = %g;   b = %g\n",
       2r, r, x, y, a, b)
   plot([a, 0, -a, 0, a], [R - b, R, R - b, R - 2b, R-b], color=:blue, lw=0.5)
   circle(0, 0, R)
   circle(x, y, r, :green)
   (ax, ay) = (336/125, -527/125)
   @printf("ax = %g; ay = %g\n", ax, ay)
   θ = atand((R + ay)/ax)
   println(θ)
   circle(0, -R, 2R - 2b, :orange, beginangle=θ, endangle=180-θ)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(0, R - b, " R-b")
       point(0, R - 2b, " R-2b")
       point(a, R - b, "(a,R-b) ", :blue, :right, :vcenter)
       point(x, y, "(x,y)")
       point(ax, ay, "(ax,ay)")
       point(0, R, " R", :red, :left, :bottom, delta=delta/2)
   end
end;