算額（その79）
福島県田村郡田母神村（現・郡山市田村町）の馬頭観世音
今有如図大円, 内設弦, 容乾坤二円, 欲使黒積至多。
乾円径若干, 問得坤円径術如何。
外円の中に，乾円，坤円が入っている。黒い部分の面積が最大になるとき，坤円の径を乾円の径で表わせ。
90度回転させて，上半分を考える。外円の半径を 1 とし，2円の接点の x 座標を a とする。
using SymPy
@syms x::real, a::positive;
乾円方程式 = sqrt(1 - x^2)
乾円面積 = integrate(乾円方程式, (x, -1, a))
坤円面積 = ((a + 1)/2)^2*PI/2
黒面積 = 乾円面積 - 坤円面積
f = diff(黒面積, a) # 導関数
using Plots
plot(f, aspectratio=:none, size=(400, 200))
hline!([0])
vline!([(16 - pi^2)/(pi^2 + 16)])
a0 = solve(f)[1] # 導関数 = 0 を解く
a0 |> println
(16 - pi^2)/(pi^2 + 16)
乾円径 = 1 + a0
坤円径 = 1 - a0
坤円径 / 乾円径 |> simplify |> println
pi^2/16
using Plots
function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
θ = beginangle:0.1:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
if fill
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5,
seriestype=:shape, fillcolor=color)
else
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end
end;
function polygon(x, y, color=:black)
plot!(x, y, seriestype=:shape, fillcolor=color)
end;
function point(x, y, string="", color=:green, position=:left, vertical=:top; fontsize=10, mark=true)
mark && scatter!([x], [y], color=color, markerstrokewidth=0)
annotate!(x, y, text(string, fontsize, vertical, position, color))
end;
function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
a = (16 - pi^2)/(16 + pi^2)
θ0 = acos(a)/pi*180
θ = θ0:0.1:360-θ0
x = cosd.(θ)
y = sind.(θ)
append!(x, a)
append!(y, 0)
plot()
# segment(a, sqrt(1-a^2), a, -sqrt(1-a^2))
if more
polygon(y, -x, :black)
circle(0, 0, 1, :black)
circle(0, -(a-1)/2, (a+1)/2, :white, fill=true)
circle(0, -(a+1)/2, (1-a)/2, :red)
point(0, -(a-1)/2, "乾円", :blue, :center, mark=false)
point(0, -(a+1)/2, "坤円", :red, :center, mark=false)
else
polygon(x, y, :black)
circle(0, 0, 1, :black)
circle((a-1)/2, 0, (a+1)/2, :white, fill=true)
circle((a+1)/2, 0, (1-a)/2, :red)
point(a, 0, "a", :black)
plot!(ylims=(0, 1.05))
end
end;
