算額（その232）
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（222)
長野県東筑摩郡筑北村青柳 碩水寺豊川稲荷社 慶応4年(1868)
団扇の中に大円 2 個、小円 １ 個がある。扇の径が 9 寸、大円の径が 4寸5分のとき、小円の径はいかほどか。
大円の半径と中心座標を r1, (x1, y1), r1 = r0//2, y1 = 0
小円の半径と中心座標を r2, (0, r0 - r2)
団扇面を切り欠く円の中心座標を (0, y), y < 0
団扇面の円と扇面を切り欠く円の交点座標を (x0, y0)
とし，連立方程式を解く。
なお，団扇面を切り欠く円の半径，中心については条件式が与えられていないので解は y を含む式になる。
include("julia-source.txt")
using SymPy
@syms r0::positive, r1::positive, y1,
r2::pisitive, r::positive, y::negative,
x0::positive, y0::negative;
r0 = 9
r1 = r0 // 2
y1 = 0
y0 = -sqrt(r0^2 - x0^2)
eq1 = r1^2 + (r0 - r2 -y1)^2 - (r1 + r2)^2
eq3 = r1^2 + (y1 - y)^2 - (r + r1)^2
eq5 = sqrt(r^2 - x0^2) + y + sqrt(r0^2 - x0^2)
res = solve([eq1, eq3, eq5], (r2, r, x0))
1-element Vector{Tuple{Sym, Sym, Sym}}:
(3, sqrt(4*y^2 + 81)/2 - 9/2, -9*sqrt(12*y^2 - 18*sqrt(4*y^2 + 81) - 162)/(4*y))
たとえば y = -9 のとき，
r2 = 3
r = 5.562305898749054
x0 = 5.2900672706322585
になる。
y = -9
(3, sqrt(4*y^2 + 81)/2 - 9/2, -9*sqrt(12*y^2 - 18*sqrt(4*y^2 + 81) - 162)/(4*y))
(3, 5.562305898749054, 5.2900672706322585)
小円半径 = 3.00000
団扇面の円を切り欠く円の中心 = (0, -9.00000); 半径 = 5.56231
2円の交点 = (5.29007, -7.28115)
using Printf
using Plots
function draw(y, more)
pyplot(size=(500, 500), grid=false, aspect_ratio=1, label="", fontfamily="IPAMincho")
r0 = 9
r1 = 9//2
y1 = 0
(r2, r, x0) = (3, sqrt(4*y^2 + 81)/2 - 9/2, -9*sqrt(12*y^2 - 18*sqrt(4*y^2 + 81) - 162)/(4*y))
y0 = -sqrt(r0^2 - x0^2)
@printf("小円半径 = %.5f\n", r2)
@printf("団扇面の円を切り欠く円の中心 = (0, %.5f); 半径 = %.5f\n", y, r)
@printf("2円の交点 = (%.5f, %.5f)\n", x0, y0)
θ = Int(round(atand((y0 - y)/x0), digits=0))
plot()
circle(0, 0, r0, :black)
circle(0, r0 - r2, r2)
circle(r1, y1, r1, :blue)
circle(-r1, y1, r1, :blue)
circle(0, y, r, :black, beginangle=θ, endangle=180 - θ)
for deg = 90:10:180-θ
plot!([-r*cosd(deg), 0, r*cosd(deg)], [r*sind(deg) + y, -r0, r*sind(deg) + y], color=:black, lw=2)
end
if more
point(0, r0, " r0", :black, :left, :bottom)
point(0, r0 - r2, " r0-r2", :red, :left, :bottom)
point(x0, y0, "(x0,y0)")
point(r1, y1, " r1", :blue)
point(r0, y1, " r0 ", :black, :right)
point(0.2r1, (2r0 - r2)/2, "小円", :red, mark=false)
point(r1, r1/2, "大円", :blue, mark=false)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
