裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

算額(その198)

2023年04月14日 | Julia

算額(その198)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1(92)
長野県下高井郡木島平村穂高 一川谷大元神社 文化8年(1811)

問4 全円の中に大円,中円,小円が入っている。大円,中円,小円の径が与えられたとき,全円の径を求めよ。

図のように,全円,大円,中円,小円の半径を R, r1, r2, r3 とし,小円の中心座標を (x3, y3) とする。連立方程式を立て,R, x3, y3 について解く。

using SymPy

@syms r1::positive, r2::positive, r3::positive, R::positive,
     x3::positive, y3::negative;

eq1 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2
eq2 = x3^2 + y3^2 - (R - r3)^2
eq3 = x3^2 + (R - 2r1 - r2 - y3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (R, x3, y3));

println("R = ", res[1][1])
println("x3 = ", res[1][2])
println("y3 = ", res[1][3])

   R = r1^2*(r1 + r2 + r3)/(r1^2 + r1*r2 - r2*r3)
   x3 = 2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 + r2 + r3)/(r1 + r2)
   y3 = -(r1^4 + 2*r1^3*r2 + r1^2*r2^2 - 3*r1^2*r2*r3 - 3*r1*r2^2*r3 - r1*r2*r3^2 + r2^2*r3^2)/((r1 + r2)*(r1^2 + r1*r2 - r2*r3))

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:blue, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(r1, r2, r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = r1^2*(r1 + r2 + r3)/(r1^2 + r1*r2 - r2*r3)
   x3 = 2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 + r2 + r3)/(r1 + r2)
   y3 = -(r1^4 + 2*r1^3*r2 + r1^2*r2^2 - 3*r1^2*r2*r3 - 3*r1*r2^2*r3 - r1*r2*r3^2 + r2^2*r3^2)/((r1 + r2)*(r1^2 + r1*r2 - r2*r3))
   @printf("R = %.5f,  x3 = %.5f,  y3 = %.5f\n", R, x3, y3)
   plot()
   circle(0, 0, R)
   circle(0, R - r1, r1, :blue)
   circle(x3, y3, r3, :green)
   circle(-x3, y3, r3, :green)
   circle(0, R - 2r1 - r2, r2, :orange)
   if more == true
       point(0, R - r1, " R-r1", :blue)
       point(0, R - 2r1 - r2, " R-2r1-r2", :orange)
       point(x3, y3, "(x3,y3)", :green)
       point(0.4r1, R - 0.6r1, "大円", :blue, mark=false)
       point(0.4r2, R - 2r1 - 0.6r2, "中円", :orange, mark=false)
       point(x3 + 0.4r3, y3 + 0.4r3, "小円", :green, mark=false)
       point(0.75R, 0.75R, "全円", :red, mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

draw(6, 2, 4, false)

   R = 10.80000,  x3 = 6.00000,  y3 = -3.20000

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村