算額（その278）
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(239)
長野県松本市筑摩 筑摩神社 明治6年(1873)
等脚台形に2本の斜線を引き，分割された区画に2個の等円が入っている。上底，下底がそれぞれ2寸4分，1尺2寸のとき，等円の径をもとめよ。
小数点を避けるため，長さを「分」単位で表す。すなわち，上底，下底をそれぞれ 24（分），120（分）とし，y軸を対称軸とする。
等円の半径を r，台形の高さを y1, 上底と下底の長さを 2a, 2b，斜線と脚の交点座標を(x2, y2) とし，図のように記号を定め，方程式を解く。
include("julia-source.txt");
using SymPy
@syms y1::positive, x2::positive, y2::positive, r::positive;
@syms y1, x2, y2, r, a, b
(a, b) = (24, 120) .// 2
eq1 = (y1 - y2)*(b - a) - y1*(x2 - 12)
eq2 = distance(a, y1, b, 0, 0, y1 - r) - r^2
eq3 = distance(-b, 0, x2, y2, 0, y1 - r) - r^2
eq4 = distance(-b, 0, x2, y2, 0, r) - r^2;
eq2, eq3, eq4 は円の中心から斜線までの平方距離であるが，分数式になっておりそのままでは solve() で解けないので，分子=0 として解く。
@syms d
apart(eq2, d) |> numerator |> println
apart(eq3, d) |> numerator |> println
apart(eq4, d) |> numerator |> println
-r^2*y1^2 + 1152*r*y1 + 144*y1^2
-r^2*y2^2 - 2*r*x2^2*y1 - 240*r*x2*y1 + 120*r*x2*y2 - 7200*r*y1 + 7200*r*y2 + x2^2*y1^2 + 120*x2*y1^2 - 120*x2*y1*y2 + 3600*y1^2 - 7200*y1*y2 + 3600*y2^2
-r^2*y2^2 - 120*r*x2*y2 - 7200*r*y2 + 3600*y2^2
eq1 = (y1 - y2)*(b - a) - y1*(x2 - 12)
eq2 = -r^2*y1^2 + 1152*r*y1 + 144*y1^2
eq3 = -r^2*y2^2 - 2*r*x2^2*y1 - 240*r*x2*y1 + 120*r*x2*y2 - 7200*r*y1 + 7200*r*y2 + x2^2*y1^2 + 120*x2*y1^2 - 120*x2*y1*y2 + 3600*y1^2 - 7200*y1*y2 + 3600*y2^2
eq4 = -r^2*y2^2 - 120*r*x2*y2 - 7200*r*y2 + 3600*y2^2
solve([eq1, eq2, eq3, eq4], (y1, x2, y2, r))
9-element Vector{NTuple{4, Sym}}:
(-90, 180/7, -450/7, -20)
(-20, -60, -50, -60)
(0, x2, 0, r)
(20, -60, 50, 60)
(90, 180/7, 450/7, 20)
(-24*sqrt(5), 0, -30*sqrt(5), -12*sqrt(5))
(-24*sqrt(5), 60, 0, -12*sqrt(5))
(24*sqrt(5), 0, 30*sqrt(5), 12*sqrt(5))
(24*sqrt(5), 60, 0, 12*sqrt(5))
5 番目の組が適解である。
すなわち等円の半径が 20 であり，元の単位での直径は 20/10*2 = 4 寸である。
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b) = (60, 12)
(y1, x2, y2, r) = (90, 180/7, 450/7, 20)
@printf("y1 = %.6f; x2 = %.6f; y2 = %.6f; r = %.6f\n", y1, x2, y2, r)
plot([a, b, -b, -a, a], [0, y1, y1, 0, 0], color=:black, lw=0.5)
segment(-a, 0, x2, y2)
segment(a, 0, -x2, y2)
circle(0, y1 - r, r, :blue)
circle(0, r, r, :blue)
if more
point(0, y1 - r, " y1-r")
point(0, r, " r")
point(x2, y2, " (x2,y2)")
point(0, y1, " y1", :green, :left, :top)
point(a, 0, " a", :green, :left, :bottom)
point(b, y1, " (b, y1)")
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
