算額（その1547）
神壁算法 播州垂井 住吉大明神社 寛政12年(1800)
藤田貞資門人 播州小野 高瀬恒右衛門信之
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf
キーワード：円6個，面積
#Julia, #SymPy, #算額, #和算，#数学
大円 3 個が交わっており，小円 3 個を容れる。黒積が 3329.29のとき，小円の直径はいかほどか。
大円の半径と中心座標を r1, (0, r1), (√3r1, -r1/2)
小円の半径と中心座標を r2, (x2, x2/√3)
とおく。
まず，r2 が与えられたとき，r1 を求める。
include("julia-source.txt");
# # julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
@syms r1::positive, x1::positive, y1::negative, r2::positive, x2::positive, y2::negative
x1 = √Sym(3)r1/2
y1 = -r1/2
y2 = -x2/√Sym(3)
eq1 = x1^2 + (2r1 - r2 - y1)^2 - (r1 + r2)^2
eq2 = x2^2 + (r1 - y2)^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r1, x2))[1]
(7*r2/6, 2*sqrt(3)*r2/3)
大円の半径は小円の半径の 7/6 倍である。
黒積は，図の灰色部分の 6 倍である。
灰色部分の面積は，(大円の面積/2 - 小円の面積/2) - 2(大円の面積/6 - 正三角形OBCの面積) である。
@syms r1, r2, 黒積, d
r1 = 7r2/6
eq = 6*((PI*r1^2/2 - PI*r2^2/2) - 2(PI*r1^2/6 - r1^2*√Sym(3)/4)) - 黒積
eq |> println
-59*pi*r2^2/36 + 49*sqrt(3)*r2^2/12 - 黒積
方程式を解き，小円の半径 r2 = 6*sqrt(黒積)/sqrt(-59*pi + 147*sqrt(3)) である。
黒積が 3329.29 のとき，小円の直径は 12*sqrt(黒積)/sqrt(-59*pi + 147*sqrt(3)) = 83.20006146161401 である。
ans_r2 = solve(eq, r2)[2]
ans_r2 |> println
6*sqrt(黒積)/sqrt(-59*pi + 147*sqrt(3))
# r2
黒積 = 3329.29
12*sqrt(黒積)/sqrt(-59*pi + 147*sqrt(3))
83.20006146161401
「術」は「sqrt(3329.29/(sqrt(64827) - 59*π))*12 = 83.20006146161398」，上の結果と一致する。
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = 83.20006146161401/2
(r1, x2) = (7*r2/6, 2*sqrt(3)*r2/3)
y2 = -x2/√3
x1 = √3r1
y1 = -r1/2
円周率 = 3.16
円周率 = π
S = 6(円周率*r1^2/2 - 円周率*r2^2/2 - 2(円周率*r1^2/6 - r1^2*√3/4))
println("r1 = $r1, r2 = $r2 のとき，黒積は $S")
plot()
rotate(0, r1, r1)
rotate(0, 2r1 - r2, r2, :blue)
θ = -30:0.1:90
x = r1.*cosd.(θ)
y = r1.*sind.(θ) .+ r1
θ = 90:-0.1:-90
append!(x, r2.*cosd.(θ))
append!(y, r2.*sind.(θ) .+ (2r1 - r2))
θ = 150:-0.1:90
append!(x, r1.*cosd.(θ) .+ √3r1/2)
append!(y, r1.*sind.(θ) .- r1/2)
plot!(x, y, color=:gray80, seriestype=:shape, fillcolor=:gray80, lw=0.5)
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, r1, "大円:r1,(0,r1)", :red, :center, delta=-delta/2)
point(0, 2r1 - r2, "小円:r2,(0,2r1-r2)", :blue, :center, :bottom, delta=delta)
point(0, 0, "O", :green, :center, delta=-1.5delta)
point(r1*sind(30), r1 - r1*cosd(30), "A", :green, :left, delta=-delta)
point(r1*sind(60), r1 - r1*cosd(60), "B", :green, :left, delta=-delta)
point(0, r1, "C ", :green, :right, :bottom, delta=-delta)
end
end;
draw(true)
