算額（その601）

算額（その601）

長野県東御市祢津 長明寺・大日如来堂 嘉永6年(1853)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html

外円内に弦，矢，大円，小円を入れる。大円，小円の直径がそれぞれ 16 寸，12 寸，矢の長さが 15 寸のとき，外円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (x + r1, y + r1)
小円の半径と中心座標を r2, (x - r2, y + r2)
矢の長さを「矢」，弦と y 軸の交点座標を (x, y)
として以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R, 矢, r1, r2, x, y
eq1 = (x + r1)^2 + (y + r1)^2 - (R - r1)^2
eq2 = (x - r2)^2 + (y + r2)^2 - (R - r2)^2
eq3 = x^2 + (y + 矢)^2 - R^2;
res = solve([eq1, eq2, eq3], (R, x, y))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-(r1^2 - 2*r1*矢 + r2^2 - 2*r2*矢 + 2*矢^2)/(4*(r1 + r2 - 矢)), -(r1 - r2)*(r1 + r2 - 2*矢)/(4*(r1 + r2 - 矢)), (-r1*r2 - r1*矢 - r2*矢 + 矢^2)/(2*r1 + 2*r2 - 2*矢))
    (-(2*r1^2*r2^2 - 2*r1^2*r2*矢 + r1^2*矢^2 - 2*r1*r2^2*矢 + r2^2*矢^2)/(4*r1*r2*矢), -(r1 - r2)*(2*r1*r2 - r1*矢 - r2*矢)/(4*r1*r2), (r1*r2 - r1*矢 - r2*矢 - 矢^2)/(2*矢))

2 組の解が得られるが，最初のものが適解である。

(r1, r2) = (16, 12) ./ 2
矢 = 15
(-(r1^2 - 2*r1*矢 + r2^2 - 2*r2*矢 + 2*矢^2)/(4*(r1 + r2 - 矢)), -(r1 - r2)*(r1 + r2 - 2*矢)/(4*(r1 + r2 - 矢)), (-r1*r2 - r1*矢 - r2*矢 + 矢^2)/(2*r1 + 2*r2 - 2*矢)),
(-(2*r1^2*r2^2 - 2*r1^2*r2*矢 + r1^2*矢^2 - 2*r1*r2^2*矢 + r2^2*矢^2)/(4*r1*r2*矢), -(r1 - r2)*(2*r1*r2 - r1*矢 - r2*矢)/(4*r1*r2), (r1*r2 - r1*矢 - r2*矢 - 矢^2)/(2*矢))

   ((32.5, -8.0, 16.5), (-2.4125, 1.1875, -12.9))

外円の直径は 65 寸である。

その他のパラメータは以下の通りである。
r1 = 8;  r2 = 6;  矢 = 15;  R = 32.5;  x = -8;  y = 16.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (16, 12) .// 2
   矢 = 15
   (R, x, y) = (-(r1^2 - 2*r1*矢 + r2^2 - 2*r2*矢 + 2*矢^2)/(4*(r1 + r2 - 矢)), -(r1 - r2)*(r1 + r2 - 2*矢)/(4*(r1 + r2 - 矢)), (-r1*r2 - r1*矢 - r2*矢 + 矢^2)/(2*r1 + 2*r2 - 2*矢))
   @printf("外円の直径 = %g;  r1 = %g;  r2 = %g;  矢 = %g;  R = %g;  x = %g;  y = %g\n", 2R, r1, r2, 矢, R, x, y)
   plot()
   circle(0, 0, R, :blue)
   circle(x + r1, y + r1, r1, :green)
   circle(x - r2, y + r2, r2, :red)
   segment(-sqrt(R^2 - y^2), y, sqrt(R^2 - y^2), y)
   segment(x, y, x, y + 矢, :brown, lw=2)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x + r1, y + r1, " 大円:r1\n (x+r1,y+r1)", :black, :left, :vcenter)
       point(x - r2, y + r2, "小円:r2\n(x-r2,y+r2)", :black, :center, delta=-delta/2)
       point(x, y + 矢/2, " 矢", :brown, mark=false)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
       point(0, y, " y", :black, :left, :top, delta=-delta/2)
       point(x, y, "(x,y)", :black, :center, :top, delta=-delta/2)
   end
end;