算額（その252）
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（234)
長野県北佐久郡軽井沢町峠 熊野神社 明治5年(1872)
円内の弦の上に等円 2 個と菱形が入っている。菱形の長い方の対角線の長さが 8 寸，等円の径が 4 寸，のとき図の「矢」の長さはいかほどか。
外円の半径，中心座標を r0, (0, r1) とする。
弦の位置の y 座標を y とする。「矢」の長さは y である。
右側の等円の半径，中心座標を r1, (r1, y + r1) とする。r1 = 2 である。
菱形の短い方の対角線の長さを 2b とおく。b = 4 である。
菱形の短い方の対角線の長さを 2a とおく。菱形の右の頂点は (4, 2r0 - a) である。
以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms r0::positive, r1::positive, a::positive, b::positive, y::positive;
r1 = 2
b = 4
eq1 = b^2 +(r0 - a)^2 - r0^2
eq2 = r1^2 + (r0 - r1 - y)^2 - (r0 - r1)^2
eq3 = distance(0, 2r0 - 2a, b, 2r0 - a, r1, y + r1) - r1^2;
eq3 = (distance(0, 2r0 - 2a, b, 2r0 - a, r1, y + r1) - r1^2)*(a^2 + b^2)/b |> expand |> simplify;
res = solve([eq1, eq2, eq3], (r0, a, y))
3-element Vector{Tuple{Sym, Sym, Sym}}:
(9/2, 9/2 - sqrt(17)/2, 1)
(9/2, sqrt(17)/2 + 9/2, 1)
(9*sqrt(5)/5, 2*sqrt(5), 1 + 3*sqrt(5)/5)
(9/2, 9/2 - sqrt(17)/2, 1) |> println
(9/2, sqrt(17)/2 + 9/2, 1) |> println
(9*sqrt(5)/5, 2*sqrt(5), 1 + 3*sqrt(5)/5) |> println
(4.5, 2.4384471871911697, 1)
(4.5, 6.561552812808831, 1)
(4.024922359499621, 4.47213595499958, 2.341640786499874)
a < b なので，最初の解が適切である。
このとき「矢」は 1，もとの単位では 1 寸である。
r0 = 4.5000000; a = 2.4384472; y（矢） = 1.0000000
SymPy の solve() では，r1, b を変数のままで r0, a, y を求めることはできない (計算時間が長く一向に終わらない）。
術では，等円^2 / 4(菱長 - 等円) である。
func(r1, b) = 2(r1^2//4(b - r1))
func(2, 4), func(3, 7), func(2, 7), func(1, 6)
(1//1, 9//8, 2//5, 1//10)
using Plots
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r0, a, y) = res[1]
@printf("r0 = %.7f; a = %.7f; y（矢） = %.7f\n", r0, a, y)
r1 = 2
x = 4
plot([x, 0, -x, 0, x], [2r0 - a, 2r0, 2r0 - a, 2r0 - 2a, 2r0 - a], color=:black, lw=0.5)
circle(0, r0, r0)
circle(r1, y + r1, r1, :blue)
circle(-r1, y + r1, r1, :blue)
x1 = sqrt(r0^2 - (r0 - y)^2)
segment(-x1, y, x1, y)
if more == true
point(0, 2r0, " 2r0", :red, :left, :bottom)
segment(-4, 2r0 - a, 4, 2r0 - a)
point(0, 2r0 - a, " 2r0-a", :black)
point(0, 2r0 - 2a, " 2r0-2a", :black)
point(0, r0, " r0", :red)
point(0, y, " y", :black)
point(r1, y + r1, "(r1,y+r1)", :blue)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
