算額（その468）

算額（その468）

泉州堺大寺三村明神社 文政8年
関流 林助右衛門自弘門人
泉州堺 岸喜八郎忠義

安藤洋美: 泉州の算額，桃山学院大学総合研究所紀要，第 24 巻，第 3 号，p.119-144.
https://www.andrew.ac.jp/soken/assets/wr/sokenk104_4.pdf

長方形の中に楕円と正方形と円が入っている。長方形の長辺，短辺が与えられたとき，正方形の一変の長さと円の直径を求めよ。

楕円の長径，短径を a，c，長方形の長辺，短辺をそれぞれ 2b, 2a，楕円の長径円の半径を r とする。
正方形の左の頂点の座標を (x, y) とする。
長方形内の「四辺形」が正方形であるためには x = b - a でなければならない。
以上を踏まえて連立方程式を立て解を求める。なおそれぞれの方程式は独立なので，eq1, eq2, eq3 を順次解くこともできる（x, y, r それぞれの数式解を求めることもできる）。

include("julia-source.txt")

using SymPy

@syms a::positive, b::positive, c::positive, h::positive,
     k::positive, m::positive, x::positive, y::positive,
     r::positive;

(a, b, c) = (5, 7.5, 3)
eq1 = x^2/c^2 + y^2/a^2 - 1
eq2 = distance(x, y, b - y, 0, c + r, 0) - r^2
eq3 = b - x - a

res = solve([eq1, eq2, eq3], (x, y, r))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (2.50000000000000, 2.76385399196283, 0.759366262899728)

   x = 2.5;  y = 2.76385;  r = 0.759366
   正方形の一辺の長さ = 3.55517

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, c) = (5, 7.5, 3)
   (x, y, r) = res[1]
   @printf("x = %g;  y = %g;  r = %g\n", x, y, r)
   @printf("正方形の一辺の長さ = %g\n", sqrt(y^2 + (a - y)^2))
   plot([b, b, -b, -b, b], [-a, a, a, -a, -a], color=:black, lw=0.5)
   ellipse(0, 0, c, a, color=:green)
   plot!([b - y, b, x + y, x, b - y], [0, a - y, a, y, 0], color=:blue, lw=0.5)
   plot!([y - b, -b, -x - y, -x, y - b], [0, a - y, a, y, 0], color=:blue, lw=0.5)
   plot!([b - y, b, x + y, x, b - y], [0, y - a, -a, -y, 0], color=:blue, lw=0.5)
   plot!([y - b, -b, -x - y, -x, y - b], [0, y - a, -a, -y, 0], color=:blue, lw=0.5)
   circle(c + r, 0, r)
   circle(-c - r, 0, r)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
       point(x, y, " (x, y)", :green, :left, :vcenter)
       point(b - y, 0, "   b-y", :blue, :left, delta=-delta/2)
       point(c, 0, " c ", :green, :right, delta=-delta/2)
       point(c + r, 0, " c+r", :red, :center, delta=-delta/2)
       point(b, 0, "b ", :black, :right, :bottom, delta=delta/2)
       point(0, a, " a", :black, :left, :bottom, delta=delta/2)
   else
      plot!(showaxis=false)
   end
end;