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