算額(その208)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2(113)
長野県長野市信更町田野口 清水神社 文政11年(1828)
鉤が 60 寸,股が 80 寸の鉤股弦内に径が 40 寸の全円を入れる。全円と中斜と股に接する甲円,全円と中斜と鉤に接する乙円の径を求めよ。
甲円,乙円の半径を r1, r2 とする。中心座標は (x1, r1), (r2, y2) とする。
中斜と弦は点 A で垂直に交わる。A の座標は (48*3/5, 48*4/5) である。
甲円,乙円が全円と外接する条件式と甲円,乙円の中心から中斜までの距離に関する条件式の連立方程式を解く。
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 r1::positive, x1::positive, r2::positive, y2::positive;
eq1 = (20 - x1)^2 + (20 - r1)^2 - (20 + r1)^2
eq2 = (20 - r2)^2 + (20 - y2)^2 - (20 + r2)^2
eq3 = distance(0, 0, 48*3/5, 48*4/5, x1, r1) - r1^2
eq4 = distance(0, 0, 48*3/5, 48*4/5, r2, y2) - r2^2
res = solve([eq1, eq2, eq3, eq4], (r1, x1, r2, y2))
4-element Vector{NTuple{4, Sym}}:
((20 - 20*sqrt(3))^2/80, 40 - 20*sqrt(3), 20/9, 20/3)
((20 - 20*sqrt(3))^2/80, 40 - 20*sqrt(3), 20, 60)
((20 + 20*sqrt(3))^2/80, 20*sqrt(3) + 40, 20/9, 20/3)
((20 + 20*sqrt(3))^2/80, 20*sqrt(3) + 40, 20, 60)
using Printf
for i in 1:4
@printf("r1 = %.5f, x1 = %.5f, r2 = %.5f, y2 = %.5f\n", res[i][1], res[i][2], res[i][3], res[i][4])
end
r1 = 2.67949, x1 = 5.35898, r2 = 2.22222, y2 = 6.66667
r1 = 2.67949, x1 = 5.35898, r2 = 20.00000, y2 = 60.00000
r1 = 37.32051, x1 = 74.64102, r2 = 2.22222, y2 = 6.66667
r1 = 37.32051, x1 = 74.64102, r2 = 20.00000, y2 = 60.00000
最初の組が適切解である。「甲円,乙円の直径は 5寸3分5厘あまりあり,4寸4分4厘あまりあり」
(20 - 20*sqrt(3))^2/80*2, 20/9*2
(5.358983848622452, 4.444444444444445)
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(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, x1, r2, y2) = res[1]
@printf("r1 = %.5f, x1 = %.5f, r2 = %.5f, y2 = %.5f\n", r1, x1, r2, y2)
plot([0, 80, 0, 0, 48*3/5], [0, 0, 60, 0, 48*4/5], color=1, lw=0.5)
circle(20, 20, 20)
circle(x1, r1, r1, :blue)
circle(r2, y2, r2, :green)
if more == true
point(20, 20, "(20,20) 全円", :red)
point(x1, r1, " (x1,r1) 甲円", :blue)
point(r2, y2, " (r2,y2) 乙円", :green, :left, :bottom)
point(48*3/5, 48*4/5, " A (48*3/5, 48*4/5)", :orange, :left, :bottom)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;