算額（その6）
香川県高松市 石清尾八幡宮 文化三年
http://www.wasan.jp/kagawa/iwaseo.html
千葉県富津市 吾妻神社 明治十年
http://www.wasan.jp/chiba/azuma.html
- 外側の円の半径を 1 とする。
- 内部の三角形は正三角形である。
- 円 A の半径を 0.17 とする。
- それぞれの円と正三角形は互いに接している。
- 円 B の中心 (x1, y1)，半径 r1 を求める。
using SymPy
@syms x, y, r::positive, r1::positive, r2::positive, x1::positive, y1::positive, x2::positive, y2::positive;
eq1 = (1 - r2 - y1)^2 + x1^2 - (r1 + r2)^2; # 円 A，円 B が接する
eq2 = y2 - y1 - (x2 - x1)/√Sym(3); # 円 B の中心と接点 C の直線の方程式
eq3 = y2 - (1 - 2r2) + x2√Sym(3); # 接点 C を通る正三角形の辺の方程式
res = solve([eq2, eq3], (x2, y2)); # 接点 x2, y2 の座標
eq4 = (res[y2] - y1)^2 + (res[x2] - x1)^2 - r1^2; # BC が半径 r1
eq5 = sqrt(x1^2 + y1^2) + r1 - 1; # 外側の円と円 B が接する
a = 0.17;
EQ1 = eq1(r2 => a); # 方程式は繊細で，r2 = 0.17 でないと解けない？
EQ4 = eq4(r2 => a);
res2 = solve([EQ1, EQ4, eq5], (x1, y1, r1)) # x1, y1, r1 を求める
1-element Vector{Tuple{Sym, Sym, Sym}}:
(0.396484958572158, 0.575520654261460, 0.301126373472639)
using Plots
function circle(ox, oy, r; color=:red)
θ = 0:0.01:2pi
x = r.*cos.(θ)
y = r.*sin.(θ)
plot!(ox .+ x, oy .+ y, color=color)
end;
function triangle(ox, oy, r; color=:blue)
θ = [90, 210, 330, 90]
x = r .* cosd.(θ)
y = r .* sind.(θ)
plot!(ox .+ x, oy .+ y, color=color)
end;
function point(x, y, string="", position=:left)
scatter!([x], [y])
annotate!(x, y, text(string, 10, position))
end;
function draw(more=false)
pyplot(size=(500, 500), aspectratio=1, label="", fontfamily="IPAMincho")
r, r2 = (1, 0.17)
x1, y1, r1 = (0.396484958572158, 0.575520654261460, 0.301126373472639)
plot()
circle(0, 0, r)
triangle(0, 0, r - 2r2)
circle(0, r - r2, r2, color=:black)
circle((r - r2)*cosd(330), (r - r2)*sind(330), r2, color=:black)
circle((r - r2)*cosd(210), (r - r2)*sind(210), r2, color=:black)
circle(x1, y1, r1)
circle(-x1, y1, r1)
circle(x1+2*r1*sind(30), y1-2*r1*cosd(30), r1)
circle(-(x1+2*r1*sind(30)), y1-2*r1*cosd(30), r1)
circle(-r1, -r1 - (r - 2r2)*sind(30), r1)
circle(r1, -r1 - (r - 2r2)*sind(30), r1)
if more
point(0, 0, " O")
point(0, 1 - r2, " A")
point(x1, y1, "B ", :right)
point(x1, y1, " (x1, y1)")
x2 = x1/4 + sqrt(3)*(1/4 - r2/2 - y1/4)
y2 = (3*y1 + 1 - sqrt(3)*x1)/4 - r2/2
point(x2, y2, "C ", :right)
point(x2, y2, " (x2, y2)")
plot!([x1, x2], [y1, y2])
annotate!((x1 + x2)/2, (y1 + y2)/1.8, text("r1", 10, :right))
end
end;
