算額（その445）

算額（その445）

十二 群馬県高崎市石原町 清水寺 享和3年(1803)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

享和3年癸亥4月(1803) 上州中里邑 黒崎九兵祀則
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

キーワード：円4個，外円，菱形
#Julia, #SymPy, #算額, #和算

外円内に甲円 1 個，乙円 2 個，菱形 1 個が入っている。外円と甲円の直径がそれぞれ 7 寸，3 寸のとき，乙円の直径を求めよ。

外円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (r0 - r1, 0)
乙円の半径と中心座標を r2, (x2, y2)
とおき，連立方程式の解を求める。

include("julia-source.txt")

using SymPy

@syms r0::positive, r1::positive,
     r2::positive, x2::positive, y2::positive;

(r0, r1) = (7, 3) .// 2
eq1 = x2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + y2^2 - (r0 - r2)^2
eq3 = distance(0, r0 - 2r1, sqrt(r0^2 - r1^2), -r1, x2, y2) - r2^2 
res = solve([eq1, eq2, eq3], (r2, x2, y2));

解は 1 組得られるが，係数が分数の長い式になるので，evalf() した値を示しておく。

乙円の直径は 2 * 1.2 = 2寸4分である。

(res[1][1].evalf(), res[1][2].evalf(), res[1][3].evalf())

   (1.20000000000000, 2.24499443206436, 0.500000000000000)

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1) = (7, 3) .// 2
   (r2, x2, y2) = (res[1][1].evalf(), res[1][2].evalf(), res[1][3].evalf())
   @printf("r2 = %g;  x2 = %g;  y2 = %g\n", r2, x2, y2)
   @printf("乙円の直径 = %g\n", 2r2)
   plot()
   circle(0, 0, r0, :red)
   circle(0, r0 - r1, r1, :blue)
   circle2(x2, y2, r2, :brown)
   b = sqrt(r0^2 - r1^2)
   plot!([0, b, 0, -b, 0], [-r0, -r1, r0 - 2r1, -r1, -r0], color=:orange, lw=0.5)
   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(0, r0 - r1, " 甲円:r1, (0,r0-r1)", :black, :center, delta=-delta)
       point(x2, y2, "乙円:r2,(x2,y2)", :brown, :center, delta=-delta)
       point(0, r0 - 2r1, " r0-2r1")
       point(0, -r1, " -r1")
       point(0, r0, " r0", :red, :left, :bottom)
       point(sqrt(r0^2 - r1^2), -r1, "(√(r0^2-r1^2),-r1) ", :orange, :right, :vcenter)
   end
end;