裏 RjpWiki

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

算額(その222)

2023年05月07日 | Julia

算額(その222)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3(200)
長野県北佐久郡軽井沢町峠 熊野神社 安政4年(1857)

外円の中に大円,中円,小円とそれぞれに外接する正三角形が入っている。大円は中円と小円にも外接している。大円と正三角形の1辺の長さが等しく 625 寸であるとき,中円と小円の径はいかほどか。

正三角形の底辺と,大円の中心は外円の直径上にある。

外円の半径を r0 とする。
大円の中心座標と半径が (-r0, 0), r1
中円の中心座標と半径が (x2, -r2), r2
小円の中心座標と半径が (0, r0 - r3), r3
とする。
以下の連立方程式を解く。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
   p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
   l = sympy.Line(p1, p2)
   l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す!!
end;

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

r0 = 625
r1 = r0 // 2

eq1 = (x2 + r1)^2 + r2^2 - (r1 + r2)^2
eq2 = r1^2 + (r0 -r3)^2 - (r1 + r3)^2
eq3 = distance(0, 0, r0/2, sqrt(Sym(3))/2*r0, 0, r0 - r3) - r3^2
eq4 = x2^2 + r2^2 - (r0 - r2)^2;

res = solve([eq1, eq4, eq3], (x2, r2, r3))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (625/3, 2500/9, 625/3)

中円の直径 = 555.556,  小円の直径 = 416.667;  中円の中心座標 = (208.333, -277.778)

中円の径は 2\*2500/9 = 555 + 5/9 である。術では (505 + 2/9)寸となっているが,転記ミスであろう。

小円の径は 2\*625/3 = 416 + 2/3 である。術では正しく (416 + 2/3)寸となっている。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
end;

function point(x, y, string="", color=:green, 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(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 625
   r1 = r0 // 2
   (x2, r2, r3) = res[1]
   @printf("中円の直径 = %.3f,  小円の直径 = %.3f;  中円の中心座標 = (%.3f, %.3f)\n", 2r2, 2r3, x2, -r2)
   plot([r0, r0/2, 0, r0], [0, √3r0/2, 0, 0], color=:orange, lw=0.5)
   circle(0, 0, r0, :blue)
   circle(-r1, 0, r1, :magenta)
   circle(x2, -r2, r2)
   circle(0, r0 - r3, r3, :green)
   if more == true
       point(r0, 0, " r0", :blue, :left, :bottom)
       point(0, r0 - r3, " ro-r3")
       point(0, r0, " r0", :blue, :left, :bottom)
       point(-r1, 0, " -r1", :magenta, :left, :bottom)
       point(r0/2, √3*r0/2, "(r0/2,√3/2*r0)", :orange, :left, :bottom)
       point(x2, -r2, "(x2,-r2)", :red)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

 

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

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

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