算額（その1630）
福島県三春町大字七草木字松山（旧七草木村）若草木神社 明治11年(1878)
～落書き帳「○△□」～ 927. 『算法天生法指南』巻之二（その５）
http://streetwasan.web.fc2.com/math20.04.25.2.html
キーワード：円2個，半円，長方形
#Julia, #Julia, #SymPy, #算額, #和算, #数学
長方形の中に大半円，中円，小円を容れる。長方形の短辺の長さが与えられたとき，中円と小円の直径を（単に数値として）かけ合わせた数を求める術を述べよ。
大円の半径と中心座標を r1, (x1, r1)
中円の半径と中心座標を r2, (r2, r2)
小円の半径と中心座標を r3, (x3, 2r1 - r3)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r1::positive, x1::positive, r2::positive,
r3::positive, x3::positive, 直平::positive,
直長::positive
eq1 = (x1 - r2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x3 - r2)^2 + (2r1 - r3 - r2)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (x1, r3, x3))[2] # 2 of 2
((2*r1^(3/2)*sqrt(r2) - 4*sqrt(r1)*r2^(3/2) + r1*r2 - 2*r2^2)/(r1 - 2*r2), r1^2/(4*r2), -sqrt(r1)*(r1 - 2*r2)/sqrt(r2) + r2)
res[2] |> println
r1^2/(4*r2)
r3 = r1^2/4r2 ゆえ，2r2 * 2r3 = r1^2 である。
つまり，中円の直径と小円の直径の積は大円の半径の 2 乗に等しい。
大円の直径（長方形の短辺）だけが与えられたのでは図は確定しない。
中円の直径も与えられると，小円の直径が確定する。
function draw(直平, r2, more)
r1 = 直平/2
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x1, r3, x3) = ((2*r1^(3/2)*sqrt(r2) - 4*sqrt(r1)*r2^(3/2) + r1*r2 - 2*r2^2)/(r1 - 2*r2), r1^2/(4*r2), -sqrt(r1)*(r1 - 2*r2)/sqrt(r2) + r2)
直長 = x1
str = @sprintf("大円径=%g，中円径=%g\n小円径=%g", 直平, 2r2, 2r3)
p = plot([0, 直長, 直長, 0, 0], [0, 0, 直平, 直平, 0], color=:green, lw=0.5)
circle(x1, r1, r1, beginangle=90, endangle=270)
circle(r2, r2, r2, :blue)
circle(x3, 2r1 - r3, r3, :magenta)
point(0, 直平, str, :black, :left, :bottom, delta=0.1, deltax=1, mark=false)
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(x1, r1, "大円:r1,(x1,r1)", :red, :right, :vcenter, deltax=-delta)
point(r2, r2, "中円:r2,(r2,r2)", :blue, :center, delta=-delta)
point(x3, 2r1 - r3, "小円:r3,(x3,2r1-r3)", :magenta, :center, delta=-delta)
else
plot!(xlims=(0, 16), ylims=(0, 12.5))
end
return p
end;
draw(12, 3.1, true)
