算額（その252）

算額（その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;