算額（その471）

算額（その471）

宮城県角田市横倉 愛宕神社 明治15年(1882)1月
http://www.wasan.jp/miyagi/yokokuraatago.html

徳竹亜紀子，谷垣美保: 2021年度の算額調査，仙台高等専門学校名取キャンパス 研究紀要，第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf

二等辺三角形（圭）内に全円と菱形，および等円 2 個が入っている。全円の直径が 14 寸のとき，等円の直径を求めよ。

二等辺三角形の底辺の長さを 2a, 高さを b とする。また，菱形の右端の頂点座標を(c, r1) とする。
全円の半径と中心座標を r1，(0, r1)
等円の半径と中心座標を r2, (x2, 0)
とし，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms a::positive, b::positive, c::positive,
     r1::positive, r2::positive, x2::positive;
r1 = 14//2
eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (a - x2)r1 - r2*a
eq3 = (b - r1)a - b*c
eq4 = ((b - r1)a/r1)^2 - (b^2 + a^2)
eq5 = r1/(a - c) - b/a
eq5 = r1*a - (a - c)b
#eq5 = (2a +2sqrt(a^2 + b^2))r1/2 - a*b
#eq5 = distance(a, 0, 0, b, x2, r2) - r2^2
res = solve([eq1, eq2, eq3, eq4], (a, c, r2, x2))

   2-element Vector{NTuple{4, Sym}}:
    (sqrt(3*b - 2*sqrt(2)*sqrt(b^2 - 21*b + 98) - 28)*(7*b + 7*sqrt(2*b^2 - 42*b + 196) - 98)/(sqrt(b)*(b - 14)), sqrt(3*b - 2*sqrt(2)*sqrt(b^2 - 21*b + 98) - 28)*(7*b^2 + 7*b*sqrt(2*b^2 - 42*b + 196) - 147*b - 49*sqrt(2*b^2 - 42*b + 196) + 686)/(b^(3/2)*(b - 14)), (21*b - 14*sqrt(2*b^2 - 42*b + 196) - 196)/b, sqrt(588*b - 392*sqrt(2)*sqrt(b^2 - 21*b + 98) - 5488)/sqrt(b))
    (sqrt(3*b + 2*sqrt(2)*sqrt(b^2 - 21*b + 98) - 28)*(7*b - 7*sqrt(2*b^2 - 42*b + 196) - 98)/(sqrt(b)*(b - 14)), sqrt(3*b + 2*sqrt(2)*sqrt(b^2 - 21*b + 98) - 28)*(7*b^2 - 7*b*sqrt(2*b^2 - 42*b + 196) - 147*b + 49*sqrt(2*b^2 - 42*b + 196) + 686)/(b^(3/2)*(b - 14)), (21*b + 14*sqrt(2*b^2 - 42*b + 196) - 196)/b, sqrt(588*b + 392*sqrt(2)*sqrt(b^2 - 21*b + 98) - 5488)/sqrt(b))

2 組の解が得られるが，最初のものが適解である。しかし，解はすべて b が含まれている。つまり，b は任意の値を取りうるということである。

問の中の条件として「全円の直径が 14 寸のとき」とあるだけで b に関する条件がない。
図において，二等辺三角形の斜辺は全円の接線になるので，b が小さくなれば a は大きくなり，それに連れて等円の半径も大きくなる。「逆も真なり」である。

res[1][3] |> println

   (21*b - 14*sqrt(2*b^2 - 42*b + 196) - 196)/b

術では，「等円の直径は，全円の直径を 2 倍して 7 で割る」としている。等円の半径が 2 になるときの b は 784/31 = 25.2903225806452 である。

solve(res[1][3] - 2)[1].evalf() |> println

   25.2903225806452

res[1][3](b => 784/31) |> println

   2.00000000000000

等円の直径 = 4;  a = 10.4766;  b = 25.2903;  c = 7.57686;  r2= 2;  x2= 7.48331

たとえば b = 17 のときは，等円の直径は 7 になる。
等円の直径 = 7;  a = 19.799;  b = 16;  c = 11.1369;  r2= 3.5;  x2= 9.89949
これを術もどきでいうと「等円の直径は，全円の直径を 2 倍して 4 で割る」ことになってしまう。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 14//2
   b = 784/31  # 25.2903225806452
   (a, c, r2, x2) = (sqrt(3*b - 2*sqrt(2)*sqrt(b^2 - 21*b + 98) - 28)*(7*b + 7*sqrt(2*b^2 - 42*b + 196) - 98)/(sqrt(b)*(b - 14)), sqrt(3*b - 2*sqrt(2)*sqrt(b^2 - 21*b + 98) - 28)*(7*b^2 + 7*b*sqrt(2*b^2 - 42*b + 196) - 147*b - 49*sqrt(2*b^2 - 42*b + 196) + 686)/(b^(3/2)*(b - 14)), (21*b - 14*sqrt(2*b^2 - 42*b + 196) - 196)/b, sqrt(588*b - 392*sqrt(2)*sqrt(b^2 - 21*b + 98) - 5488)/sqrt(b))
   @printf("等円の直径 = %g;  a = %g;  b = %g;  c = %g;  r2= %g;  x2= %g\n", 2r2, a, b, c, r2, x2)
   plot([a, 0, -a, 0], [0, b, 0, 0], color=:black, lw=0.5)
   plot!([c, 0, -c, 0, c], [r1, 2r1, r1, 0, r1], color=:blue, lw=0.5)
   circle(0, r1, r1, :red)
   circle(x2, r2, r2, :magenta)
   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(a, 0, "a", :black, :left, :bottom, delta=delta/2)
       point(0, b, " b", :black, :left, :vcenter)
       point(c, r1, " (c,r1)", :blue, :left, :vcenter)
       point(0, r1, " 全円:r1,(0,r1)", :red, :left, :vcenter)
       point(x2, r2, "等円:r2,(x2,r2)", :magenta, :center, :top, delta=-delta)
       point(0, 2r1, " 2r1", :red, :left, :bottom, delta=delta/2)
       
   end
end;