算額（その771）

算額（その771）

神壁算法 天明5年乙巳 關流四傅藤田貞資六弟子 筑後州久留米 佐田甚助方之
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

和算問題あれこれ2 令和６年４月の問題-No.３
https://gunmawasan.web.fc2.com/k-n-mondai.html

円弧の弦の上に，底辺が弦の上で頂点が円弧上にある甲，乙，丙の正三角形が載っている。一辺の長さがそれぞれ 206 寸 1 分，229 寸，183 寸 2 分のとき，弦の長さはいかほどか。

正三角形の底辺の 2 つの x 座標，頂点の x 座標を，左から（丙，乙，甲の順で），a, b, c, d, e, f, g とする。
弦と y 軸の交点の y 座標を h とする。
頂点の y 座標 は h + √3(b - a)，h + √3(d - c)，h + √3(f - e) とする。
A = g - e，B = e - c，C = c - a とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     c::positive, d::positive,
     e::positive, f::positive,
     g::negative, h::positive,
     R::positive
@syms a, b, c, d, e, f, g, h, R, A, B, C
s3 = sqrt(Sym(3))
eq1 = g - e - A
eq2 = e - c - B
eq3 = c - a - C
eq4 = b - (a + c)/2
eq5 = d - (c + e)/2
eq6 = f - (e + g)/2
eq7 = (h + (b - a)s3)^2 + b^2 - R^2
eq8 = (h + (d - c)s3)^2 + d^2 - R^2
eq9 = (h + (f - e)s3)^2 + f^2 - R^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9],
   (R, h, a, b, c, d, e, f, g))

   2-element Vector{NTuple{9, Sym{PyCall.PyObject}}}:
    (-sqrt(3)*sqrt((A^2 - A*B + B^2)*(B^2 - B*C + C^2)*(A^2 + A*B - A*C + B^2 + B*C + C^2))/(3*(A*C - B^2)), sqrt(3)*(-A^2*B - A^2*C - A*B*C - A*C^2 + B^3 - B*C^2)/(6*(A*C - B^2)), (A^2*B - A^2*C - A*B*C - A*C^2 + B^3 + 2*B^2*C - B*C^2)/(2*(A*C - B^2)), (A^2*B - A^2*C - A*B*C + B^3 + B^2*C - B*C^2)/(2*(A*C - B^2)), (B - C)*(A^2 - A*C + B^2 + B*C)/(2*(A*C - B^2)), (A - C)*(A*B - A*C + B*C)/(2*(A*C - B^2)), (A - B)*(A*B - A*C + B^2 + C^2)/(2*(A*C - B^2)), (A^2*B - A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)), (A^2*B + A^2*C - 2*A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)))
    (sqrt(3)*sqrt((A^2 - A*B + B^2)*(B^2 - B*C + C^2)*(A^2 + A*B - A*C + B^2 + B*C + C^2))/(3*(A*C - B^2)), sqrt(3)*(-A^2*B - A^2*C - A*B*C - A*C^2 + B^3 - B*C^2)/(6*(A*C - B^2)), (A^2*B - A^2*C - A*B*C - A*C^2 + B^3 + 2*B^2*C - B*C^2)/(2*(A*C - B^2)), (A^2*B - A^2*C - A*B*C + B^3 + B^2*C - B*C^2)/(2*(A*C - B^2)), (B - C)*(A^2 - A*C + B^2 + B*C)/(2*(A*C - B^2)), (A - C)*(A*B - A*C + B*C)/(2*(A*C - B^2)), (A - B)*(A*B - A*C + B^2 + C^2)/(2*(A*C - B^2)), (A^2*B - A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)), (A^2*B + A^2*C - 2*A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)))

2 組の解が得られるが，2 番目のものが適解である。

弦の長さは 2sqrt(R^2 - h^2) である。
正三角形の一辺の長さがそれぞれ 206 寸 1 分，229 寸，183 寸 2 分のとき，弦の長さは 1029.80004005632 寸である。

2sqrt(res[2][1]^2 - res[2][2]^2)(A => 206.1, B => 229.0, C => 183.2) |> println

   1029.80004005632

その他のパラメータは以下のとおりである。
R = 764.582;  h = 565.211;  a = -337.775;  b = -246.175;  c = -154.575;  d = -40.075;  e = 74.425;  f = 177.475;  g = 280.525

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (A, B, C) = (206.1, 229.0, 183.2)
   (R, h, a, b, c, d, e, f, g) = (-sqrt(3)*sqrt((A^2 - A*B + B^2)*(B^2 - B*C + C^2)*(A^2 + A*B - A*C + B^2 + B*C + C^2))/(3*(A*C - B^2)), sqrt(3)*(-A^2*B - A^2*C - A*B*C - A*C^2 + B^3 - B*C^2)/(6*(A*C - B^2)), (A^2*B - A^2*C - A*B*C - A*C^2 + B^3 + 2*B^2*C - B*C^2)/(2*(A*C - B^2)), (A^2*B - A^2*C - A*B*C + B^3 + B^2*C - B*C^2)/(2*(A*C - B^2)), (B - C)*(A^2 - A*C + B^2 + B*C)/(2*(A*C - B^2)), (A - C)*(A*B - A*C + B*C)/(2*(A*C - B^2)), (A - B)*(A*B - A*C + B^2 + C^2)/(2*(A*C - B^2)), (A^2*B - A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)), (A^2*B + A^2*C - 2*A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)))
   @printf("R = %g;  h = %g;  a = %g;  b = %g;  c = %g;  d = %g;  e = %g;  f = %g;  g = %g\n", R, h, a, b, c, d, e, f, g)
   s3 = √3
   plot([a, b, c, d, e, f, g], h .+ [0, (b - a)s3, 0, (d - c)s3, 0, (f - e)s3, 0], color=:blue, lw=0.5)
   θ = asind(h/R)
   circle(0, 0, R, beginangle=θ, endangle=180 - θ, n=1000)
   x = R*cosd(θ)
   segment(-x, h, x, h, :red)
   plot!([-x, 0, x], [h, 0, h], color=:gray50, lw=0.5)
   @printf("弦の長さ = %g\n", 2R*cosd(θ))
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point.([a, b, c, d, e, f, g], h, Char["abcdefg"...], :blue, :center, delta=-delta)
       point(x, h, "(x,h) ", :red, :right, delta=-delta)
       point(x/2+4delta, h/2, "R", :gray50, :center, :vcenter, mark=false)
       point.([b, d, f], h + 75, ["丙三角", "乙三角", "甲三角"], :blue, :center, delta=-delta, mark=false)
   end
end;