算額（その221）
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（199）
長野県北佐久郡軽井沢町峠 熊野神社 安政4年(1857)
半外円内に4本の斜線を引き，大円，中円，小円を入れる。
図のように記号を定め，方程式を解き r3 を求める。
外円の半径を r0 とする。
大円の中心座標，半径を (0, r1), r1
中円の中心座標，半径を (0, r2), r2
小円の中心座標，半径を (0, r0 - r3), r3
斜線と円の交点を (x, y) とおく。
熊野神社の算額では r0 = 100 寸のとき，r2, r3 を求めよというものであるが答えは完全ではない。
同じ問題が
県内の算額3（232）
長野県上田市塩田町五加公民館 絵堂地蔵堂 明治4年(1871)
にあり， r1 = 9 寸のとき r2, 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, x::positive, y::positive;
r0 = 100
r1 = r0 // 2
r1 = 9
r0 = 2r1
eq1 = x^2 + y^2 - r0^2
eq2 = distance(-r0, 0, x, y, 0, r0 - r3) - r3^2
eq3 = distance(-r0, 0, x, y, 0, r2) - r2^2
eq4 = distance(r0, 0, x, y, 0, r1) - r1^2
res = solve([eq1, eq2, eq3, eq4], (x, y, r2, r3))
1-element Vector{NTuple{4, Sym}}:
(126/25, 432/25, 6, 2)
ちなみに上記のようにすると処理時間が長くなる。
res = solve([eq1, eq2, eq3, eq4]) とすると，速い。
solve([eq1, eq2, eq3, eq4])
1-element Vector{Dict{Any, Any}}:
Dict(r3 => 2, y => 432/25, x => 126/25, r2 => 6)
外円の半径 = 18.000, 大円の半径 = 9.000; 中円の半径 = 6.000; 小円の半径 = 2.000
大円，中円，小円の径はそれぞれ外円の径の 1/2，1/3, 1/9 である。
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")
(x, y, r2, r3) = res[1]
r1 = 9
r0 = 2r1
@printf("外円の半径 = %.3f, 大円の半径 = %.3f; 中円の半径 = %.3f; 小円の半径 = %.3f\n", r0, r1, r2, r3)
plot([r0, -x, -r0, x, r0, -r0], [0, y, 0, y, 0, 0], color=:black, lw=0.5)
circle(0, 0, r0, :blue, beginangle=0, endangle=180)
circle(0, r1, r1, :magenta)
circle(0, r2, r2)
circle(0, r0 - r3, r3, :green)
if more == true
point(0, r1, " r1")
point(0, r2, " r2")
point(0, r0 - r3, " r3")
point(0, r0, " r0", :green, :left, :bottom)
point(-r0, 0, " -r0", :green, :left, :bottom)
point(r0, 0, " r0", :green, :left, :bottom)
point(x, y, " (x,y)", :green, :left, :bottom)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
