算額（その1503）
六十三 花泉町金沢 大門神社・大門観世音菩薩堂 明治13年(1880)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
今有如図 03064
https://w.atwiki.jp/sangaku/pages/292.html
キーワード：円8個，外円，弦
#Julia, #SymPy, #算額, #和算
大円の一部である円弧（弓形）の中に小円を 2 個と界斜を容れる。大円の径が 8 寸，矢が 2 寸，小円の直径が 1 寸のとき，界斜はいかほどか。
山村の図は右側の小円と界斜の端点の位置が違う。
「今有如図」の仮定で解を導き，答と一致することを確認した。
大円の半径と中心座標を R, (0, 0)
矢をそのまま「矢」
小円の半径と中心座標を r, (x, R - 矢 + r)
界斜の端点を (x1, R - 矢), (x2, sqrt(R^2 - x2^2))
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms R::positive, 矢::positive, x::positive, r::positive,
x1::positive, x2::negative
y = R - 矢
eq1 = x^2 + (y + r)^2 - (R - r)^2
eq2 = dist(x2, sqrt(R^2 - x2^2), x1, y, x, y + r) - r^2
eq3 = (sqrt(R^2 - x2^2) - y)/(x1 - x2) - r/x1;
# res = solve([eq1, eq2, eq3], (x, x1, x2))
function H(u)
(x, x1, x2) = u
return [
x^2 + (y + r)^2 - (R - r)^2,
dist(x2, sqrt(R^2 - x2^2), x1, y, x, y + r) - r^2,
(sqrt(R^2 - x2^2) - y)/(x1 - x2) - r/x1
]
end;
R = 8//2
矢 = 2
y = R - 矢
r = 1//2
iniv = BigFloat[2.5, 2.5, -2.6]
res = nls(H, ini=iniv)
([2.449489742783178, 2.3979157616563596, -2.597742075127723], true)
界斜は sqrt((x1 - x2)^2 + (sqrt(R^2 - x2^2) - (R - 矢))^2) により求めることができる。
大円の直径が 8 寸，矢が 2 寸，円の直径が 1 寸のとき，5.103103630798287 寸である。
(x, x1, x2) = res[1]
sqrt((x1 - x2)^2 + (sqrt(R^2 - x2^2) - (R - 矢))^2)
5.103103630798287
function draw(r1, r2, more=false)
pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
R = 8//2
矢 = 2
r = 1//2
y = R - 矢
(x, x1, x2) = res[1]
界斜 = sqrt((x1 - x2)^2 + (sqrt(R^2 - x2^2) - (R - 矢))^2)
println("界斜 = $界斜")
plot()
circle(0, 0, R, :purple)
segment(x2, sqrt(R^2 - x2^2), x1, R - 矢, :orange)
circle2(x, y + r, r, :blue)
segment(-sqrt(R^2 - y^2), y, sqrt(R^2 - y^2), y)
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(x2, sqrt(R^2 - x2^2), "(x2,sqrt(R^2-x2^2))", :green, :left, :bottom, deltax=3delta)
point(x1, y, "(x1,y)", :green, :center, delta=-delta)
point(R, 0, " R", :purple, :left, :bottom, delta=delta)
point(x, y + r, "(x,y+r)", :blue, :center, delta=-delta)
ylims!(0, R)
end
end;
draw(12/2, 8/2, true)
