算額（その1593）
長野県 諏訪神社 文化2年(1805)
深川英俊，トニー・ロスマン：聖なる数学：算額，p. 114，森北出版株式会社，2010年４月22日.
キーワード：円1個，長方形，対角線
#Julia, #SymPy, #算額, #和算, #数学
長方形の中に円と対角線を引く。長方形の長辺 AB と短辺 BC が与えられたとき，円によって切り取られる対角線の長さ PQ を求めよ。
円の半径を r
a = AB, b = BC = 2r
P, Q の座標を (x1, y1), (x2, y2)
θ = ∠BDC
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, b::positive, r::positive,
x1::positive, x2::positive, PQ::positive
r = b/2
tanθ = b/a
y1 = (a - x1)*tanθ
y2 = (a - x2)*tanθ
eq1 = (x1 - (a - r))^2 + (y1 - r)^2 - r^2
eq2 = (x2 - (a - r))^2 + (y2 - r)^2 - r^2
eq3 = sqrt((x1 - x2)^2 + (y1 - y2)^2) - PQ;
res = solve([eq1, eq2, eq3], (PQ, x1, x2))[1]; # 1 of 2
# PQ
res[1] |> println
sqrt(2)*sqrt(a)*b^(3/2)/sqrt(a^2 + b^2)
PQ の長さは sqrt((2a*b^3)/(a^2 + b^2)) である。
a = 185, b = 80 のとき，PQ = 68.2871764062511 である。
res[1](a => 185, b => 80).evalf() |> println
68.2871764062511
function draw(a, b, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x1, x2, PQ) = ((-sqrt(2)*a^(3/2)*b^(3/2) + a*(2*a^2 - a*b + b^2))/(2*(a^2 + b^2)), (sqrt(2)*a^(3/2)*b^(3/2) + a*(2*a^2 - a*b + b^2))/(2*(a^2 + b^2)), sqrt(2)*sqrt(a)*b^(3/2)/sqrt(a^2 + b^2))
r = b/2
tanθ = 2r/a
θ = atand(tanθ)
y1 = (a - x1)*tanθ
y2 = (a - x2)*tanθ
PQ2 = sqrt((x1 - x2)^2 + (y1 - y2)^2)
println("PQ = $PQ = $PQ2")
plot([0, a, a, 0, 0], [0, 0, 2r, 2r, 0], color=:magenta, lw=0.5)
circle(a - r, r, r, :blue)
segment(0, 2r, a, 0, :red)
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 - r, r, "r,(a-r,r)", :blue, :center, delta=-delta/2)
point(x1, y1, " Q:(x1,y1)", :blue, :left, :bottom, delta=delta/2)
point(x2, y2, "P:(x2,y2) ", :blue, :right, :vcenter)
point(a, b, "A:(a,b)", :magenta, :right, :bottom, delta=delta/2)
point(0, b, " B:(0,b)", :magenta, :left, :bottom, delta=delta/2)
point(0, 0, " C:(0,0)", :magenta, :left, delta=-delta/2)
point(a, 0, "D:(a,0)", :magenta, :right, delta=-delta/2)
circle(a, 0, a/20, beginangle=180-θ, endangle=180)
circle(a, 0, a/18, beginangle=180-θ, endangle=180)
point(a - a/17, 0, "θ", :red, :right, :bottom, delta=delta/2, mark=false)
ylims!(-5delta, b + 5delta)
end
end;
draw(5, 3, true)
