算額（その79）

算額（その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;