算額（その708）
八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg
キーワード：円13個，外円，弦2本
#Julia, #SymPy, #算額, #和算
外円内に水平な 2 本の弦，甲円 1 個，乙円，丙円，丁円を 2 個ずつ，戊円 4 個，己円 1 個を容れる。
己円の直径が 33 寸のとき，乙円の直径はいかほどか。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, 6r5 - R)
乙円の半径と中心座標を r2, (x2, 2r5 - R)
丙円の半径と中心座標を r3, (x3, y + r3)
丁円の半径と中心座標を r4, (x4, R - 2r6 - r4)
戊円の半径と中心座標を r5, (0, r5 - R), (0, 3r5 - R), (0, 5r5 - R), (0, 7r5 - R)
弦と y 軸の交点座標を (0, y), (0, R - 2r6)
その他の変数の関係は
R = r6 + 4r5
y = 4r5 - R
r1 = 2r5
y4 = R - 2r6 - r4
とおき，以下の連立方程式を(r6 を変数のまま）解く。
include("julia-source.txt");
using SymPy
@syms r1::positive, r2::positive, x2::positive,
r3::positive, x3::positive,
r4::positive, x4::positive, y4::positive,
r5::positive, r6::positive
@syms r1, r2, x2, r3, x3, r4, x4, y4, r5, r6
R = r6 + 4r5
y = 4r5 - R
y4 = R - 2r6 - r4
r1 = 2r5
eq1 = x2^2 + (2r5 - R)^2 - (R - r2)^2
eq2 = x3^2 + (y + r3)^2 - (R - r3)^2
eq3 = x4^2 + y4^2 - (R - r4)^2
eq4 = x2^2 + r5^2 - (r2 + r5)^2
eq5 = x3^2 + (6r5 - R - y - r3)^2 - (r1 + r3)^2
eq6 = x4^2 + (6r5 - R - y - r4)^2 - (r1 + r4)^2
eq7 = (x3 - x4)^2 + (y4 - y - r3)^2 - (r3 + r4)^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r2, x2, r3, x3, r4, x4, r5));
10 組の解が得られるが，最初の 2 つが適解である。ただし，x2, x3, x4 は符号が負の解が得られる（y 軸で線対称なので適解である）。
res[1]
(28*r6/11, -2*sqrt(21)*sqrt(r6*(143*r6 + 73*sqrt(r6^2)))/33, r6*(10757*r6 + 11023*sqrt(r6^2))/(66*(59*r6 + 73*sqrt(r6^2))), -sqrt(2)*sqrt(r6^3*(10757*r6 + 11023*sqrt(r6^2)))/(33*r6), 2*r6*(59*r6 + 73*sqrt(r6^2))/(92*r6 + 73*sqrt(r6^2)), -2*sqrt(66)*sqrt(r6/(92*r6 + 73*sqrt(r6^2)))*(59*r6 + 73*sqrt(r6^2))/33, 59*r6/66 + 73*sqrt(r6^2)/66)
res[2]
(28*r6/11, 2*sqrt(21)*sqrt(r6*(143*r6 + 73*sqrt(r6^2)))/33, r6*(10757*r6 + 11023*sqrt(r6^2))/(66*(59*r6 + 73*sqrt(r6^2))), -sqrt(2)*sqrt(r6^3*(10757*r6 + 11023*sqrt(r6^2)))/(33*r6), 2*r6*(59*r6 + 73*sqrt(r6^2))/(92*r6 + 73*sqrt(r6^2)), -2*sqrt(66)*sqrt(r6/(92*r6 + 73*sqrt(r6^2)))*(59*r6 + 73*sqrt(r6^2))/33, 59*r6/66 + 73*sqrt(r6^2)/66)
乙円の半径(r2)は 28*r6/11 となり，己円の半径(r6)の 28/11 倍である。術でも「己径二十八段以十一個除（己円の直径を28倍して11で割る）」としている。
したがって，己円の直径が 33 寸のとき，乙円の直径は 84 寸である。
その他のパラメータは以下のとおりである。
R = 148.5; r1 = 66; r2 = 42; x2 = -67.3498; r3 = 41.25; x3 = -104.355; r4 = 26.4; x4 = -83.4841; y4 = 89.1; r5 = 33; r6 = 16.5
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
#(r2, x2, r3, x3, r4, x4, r5) =
# (42, 68, 41, 104, 26, 84, 33)
r6 = 33//2
(r2, x2, r3, x3, r4, x4, r5)= (28*r6/11, -2*sqrt(21)*sqrt(r6*(143*r6 + 73*sqrt(r6^2)))/33, r6*(10757*r6 + 11023*sqrt(r6^2))/(66*(59*r6 + 73*sqrt(r6^2))), -sqrt(2)*sqrt(r6^3*(10757*r6 + 11023*sqrt(r6^2)))/(33*r6), 2*r6*(59*r6 + 73*sqrt(r6^2))/(92*r6 + 73*sqrt(r6^2)), -2*sqrt(66)*sqrt(r6/(92*r6 + 73*sqrt(r6^2)))*(59*r6 + 73*sqrt(r6^2))/33, 59*r6/66 + 73*sqrt(r6^2)/66)
R = r6 + 4r5
y = 4r5 - R
r1 = 2r5
y4 = R - 2r6 - r4
@printf("乙円の直径 = %g; 己円の直径 = %g\n", 2r2, 2r6)
@printf("R = %g; r1 = %g; r2 = %g; x2 = %g; r3 = %g; x3 = %g; r4 = %g; x4 = %g; y4 = %g; r5 = %g; r6 = %g\n",
R, r1, r2, x2, r3, x3, r4, x4, y4, r5, r6)
plot()
circle(0, 0 , R, :blue)
circle(0, 6r5 - R, r1, :brown)
circle(x2, 2r5 - R, r2, :magenta)
circle(-x2, 2r5 - R, r2, :magenta)
circle(x3, y + r3, r3, :green)
circle(-x3, y + r3, r3, :green)
circle(x4, y4, r4, :orange)
circle(-x4, y4, r4, :orange)
x1 = sqrt(R^2 - y^2)
segment(-x1, y, x1, y, :green)
circle(0, 7r5 - R, r5)
circle(0, 5r5 - R, r5)
circle(0, 3r5 - R, r5)
circle(0, r5 - R, r5)
y20 = R - 2r6
x20 = sqrt(R^2 - y20^2)
segment(-x20, y20, x20, y20, :green)
circle(0, R - r6, r6, :gray)
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(0, 6r5 - R, "甲円:r1\n(0,6r5-R)", :brown, :center, delta=-delta)
point(-x2, 2r5 - R, "乙円:r2\n(x2,2r5-R)", :magenta, :center, delta=-delta/2)
point(-x3, y + r3, "丙円:r3\n(x3,y+r3)", :green, :center, delta=-delta/2)
point(-x4, y4, "丁円:r4\n(x4,y4)", :orange, :center, delta=-delta/2)
point(0, r5 - R, "戊円:r5\n(0,r5-R)", :red, :center, delta=-delta/2)
point.(0, [3, 5, 7].*r5 .- R, "", :red)
point(0, R - r6, "己円:r6,(0,R-r6)", :black, :center, :bottom, delta=delta/2)
point(0, y, " y", :green, :left, delta=-delta/2)
point(0, R - 2r6, " R-2r6", :green, :left, :bottom, delta=delta/2)
end
end;
