算額（その149）

算額（その149）

岐阜県 田代神社 天保12年(1841)8月
http://www.wasan.jp/gifu/tasiro.html

第4問: 正方形内に四分円弧を描き，その中に円弧に接する青円，黃円，黒円を入れる。また，左上，右下に赤円と白円を含み，四分円に接する楕円を描く。黒円が与えられたとき楕円の長軸の長さを求めよ。

正方形の一辺の長さを 1 とし，以下のように記号を定め連立方程式を解く。
青円の半径を r1，中心は (1/2, 1/2)
黃円の半径を r2，中心を (x2, x2)
黒円の半径を r3，中心を (x3, x3)
左上楕円内の赤円の半径を r4, 中心は (r4, 1 - r4 - 2r5)
左上楕円内，下の白円の半径を r5，中心は (r4, 1 - 2r4 - 3r5) なお，r5 = r4 / 2
楕円と四分円の接点座標を (x0, y0)

using SymPy
@syms r1::positive, r2::positive, x2::positive, r3::positive, x3::positive,
       r4::positive, r5::positive, x0::positive, y0::positive;
eq1 = (1 - r1)^2 - 1//2  # 青円の半径
eq2 = (1 - x2)^2 + x2^2 - (1-r2)^2  # 黃円が四分円に内接する
eq3 = 2(1//2 - x2)^2 - (r1 + r2)^2  # 青円と黃円が外接する
eq4 = (1 - x3)^2 + x3^2 - (1 - r3)^2  # 黒円が四分円に内接する
eq5 = 2(x2 - x3)^2 - (r2 + r3)^2;  # 黃円と黒円が外接する
r5 = r4/2
a = r4
b = r4 + 2r5
y4 = 1 - b
eq7 = b^2*(x0 - a)/(a^2*(1 - b - y0)) + (x0 - 1) / y0  # （x0, y0） における円と楕円の接線の傾きが等しい
eq8 = (x0 - r4)^2/a^2 + (y0 - y4)^2/b^2 - 1  # (x0, y0) が楕円に周上にある
eq9 = (x0 - 1)^2 + y0^2 - 1;  # (x0, y0) が円周上にある

res = solve([eq1, eq2, eq3, eq4, eq5, eq7, eq8, eq9], (r1, r2, x2, r3, x3, r4, x0, y0))

   12-element Vector{NTuple{8, Sym}}:
    (1 - sqrt(2)/2, -1/7 + 3*sqrt(2)/14, 3/14 + 3*sqrt(2)/7, -1/287 + 17*sqrt(2)/574, 263/574 + 102*sqrt(2)/287, 1, 1, 1)
    (1 - sqrt(2)/2, -1/7 + 3*sqrt(2)/14, 3/14 + 3*sqrt(2)/7, -1/287 + 17*sqrt(2)/574, 263/574 + 102*sqrt(2)/287, -sqrt(3)/2 - sqrt(2)/2 + sqrt(6)/2 + 3/2, sqrt(6)/3 + 1, sqrt(3)/3)
    (1 - sqrt(2)/2, -1/7 + 3*sqrt(2)/14, 3/14 + 3*sqrt(2)/7, -1/287 + 17*sqrt(2)/574, 263/574 + 102*sqrt(2)/287, -sqrt(6)/2 - sqrt(3)/2 + sqrt(2)/2 + 3/2, 1 - sqrt(6)/3, sqrt(3)/3)
    (1 - sqrt(2)/2, -1/7 + 3*sqrt(2)/14, 11/14 - 3*sqrt(2)/7, -1/287 + 17*sqrt(2)/574, 311/574 - 102*sqrt(2)/287, 1, 1, 1)
    (1 - sqrt(2)/2, -1/7 + 3*sqrt(2)/14, 11/14 - 3*sqrt(2)/7, -1/287 + 17*sqrt(2)/574, 311/574 - 102*sqrt(2)/287, -sqrt(3)/2 - sqrt(2)/2 + sqrt(6)/2 + 3/2, sqrt(6)/3 + 1, sqrt(3)/3)
    (1 - sqrt(2)/2, -1/7 + 3*sqrt(2)/14, 11/14 - 3*sqrt(2)/7, -1/287 + 17*sqrt(2)/574, 311/574 - 102*sqrt(2)/287, -sqrt(6)/2 - sqrt(3)/2 + sqrt(2)/2 + 3/2, 1 - sqrt(6)/3, sqrt(3)/3)
    (1 - sqrt(2)/2, -1/7 + 3*sqrt(2)/14, -sqrt(2)*sqrt(11 - 6*sqrt(2))/7 + 1/2, 1 - sqrt(2)/2, 1/2, 1, 1, 1)
    (1 - sqrt(2)/2, -1/7 + 3*sqrt(2)/14, -sqrt(2)*sqrt(11 - 6*sqrt(2))/7 + 1/2, 1 - sqrt(2)/2, 1/2, -sqrt(3)/2 - sqrt(2)/2 + sqrt(6)/2 + 3/2, sqrt(6)/3 + 1, sqrt(3)/3)
    (1 - sqrt(2)/2, -1/7 + 3*sqrt(2)/14, -sqrt(2)*sqrt(11 - 6*sqrt(2))/7 + 1/2, 1 - sqrt(2)/2, 1/2, -sqrt(6)/2 - sqrt(3)/2 + sqrt(2)/2 + 3/2, 1 - sqrt(6)/3, sqrt(3)/3)
    (1 - sqrt(2)/2, -1/7 + 3*sqrt(2)/14, sqrt(2)*sqrt(11 - 6*sqrt(2))/7 + 1/2, 1 - sqrt(2)/2, 1/2, 1, 1, 1)
    (1 - sqrt(2)/2, -1/7 + 3*sqrt(2)/14, sqrt(2)*sqrt(11 - 6*sqrt(2))/7 + 1/2, 1 - sqrt(2)/2, 1/2, -sqrt(3)/2 - sqrt(2)/2 + sqrt(6)/2 + 3/2, sqrt(6)/3 + 1, sqrt(3)/3)
    (1 - sqrt(2)/2, -1/7 + 3*sqrt(2)/14, sqrt(2)*sqrt(11 - 6*sqrt(2))/7 + 1/2, 1 - sqrt(2)/2, 1/2, -sqrt(6)/2 - sqrt(3)/2 + sqrt(2)/2 + 3/2, 1 - sqrt(6)/3, sqrt(3)/3)

using Printf
for i = 1:12
   (r1, r2, x2, r3, x3, r4, x0, y0) = res[i]
   if x3 < x2 < 0.5 && r4 < 0.5 && x0 < 0.5
       println("i = $i")
       println("r1 = $(res[i][1])")
       println("r2 = $(res[i][2])")
       println("r3 = $(res[i][3])")
       println("r4 = $(res[i][4])")
       println("x2 = $(res[i][5])")
       println("x3 = $(res[i][6])")
       println("x0 = $(res[i][7])")
       println("y0 = $(res[i][8])")
   end
end

   i = 6
   r1 = 1 - sqrt(2)/2
   r2 = -1/7 + 3*sqrt(2)/14
   r3 = 11/14 - 3*sqrt(2)/7
   r4 = -1/287 + 17*sqrt(2)/574
   x2 = 311/574 - 102*sqrt(2)/287
   x3 = -sqrt(6)/2 - sqrt(3)/2 + sqrt(2)/2 + 3/2
   x0 = 1 - sqrt(6)/3
   y0 = sqrt(3)/3

青径は 2r1 = 2 - √2
黄径は 2r2 = (3√2 - 2) / 7
黒径は 2r3 = (17√2 - 2) / 287
赤径は 2r4 = 3 + √2 - √3 - √6
白径は 2r5 = r4 = (3 + √2 - √3 - √6) / 2
楕円の長軸は 2a + 4b = 2r4 + 4r5 = 4r4 = 2(3 + √2 - √3 - √6)
長軸 / 黒径 = 40 + 53√2 - 36√3 - 19√6

つまり，
長軸 = (40 + 53√2 - 36√3 - 19√6) × 黒径 ≒ 6.05918 × 黒径

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function ellipse(ox, oy, ra, rb; φ=0, beginangle=0, endangle=360,
                    color=:black, lty=:solid, lwd=0.5, fcolor="")
"""
(ox, oy) を中心，ra, rb を半径とする楕円（楕円弧）。
fcolor を指定すると塗りつぶし。
"""
   θ = beginangle:0.1:endangle
   if φ == 0
       if fcolor == ""
           plot!(ra .* cosd.(θ) .+ ox, rb .* sind.(θ) .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd)
       else
           plot!(ra .* cosd.(θ) .+ ox, rb .* sind.(θ) .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd,
                 seriestype=:shape, fillcolor=fcolor)
       end
   else
       x = ra .* cosd.(θ)
       y = rb .* sind.(θ)
       cosine = cosd(φ)
       sine = sind(φ)
       if fcolor == ""
           plot!(cosine .* x .- sine .* y .+ ox,
                 sine .* x .+ cosine .* y .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd)
       else
           plot!(cosine .* x .- sine .* y .+ ox,
                 sine .* x .+ cosine .* y .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd,
                 seriestype=:shape, fillcolor=fcolor)
       end
   end
end

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
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 rect(x1, y1, x2, y2, color=:pink)
   plot!([x1, x2, x2, x1, x1], [y1, y1, y2,  y2, y1], color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, x2, r3, x3, r4, x0, y0) = res[6]
   r5 = r4 / 2
   @printf("r1 = %.5f; r2 = %.5f; r3 = %.5f; r4 = %.5f; r5 = %.5f; x2 = %.5f; x3 = %.5f; x0 = %.5f; y0 = %.5f\n",
       r1, r2, r3, r4, r5, x2, x3, x0, y0)
   plot([0, 1, 1, 0, 0], [0, 0, 1, 1, 0], color=:black, lw=0.5)
   circle(1, 0, 1, :magenta, beginangle=90, endangle=180)
   circle(0, 1, 1, :magenta, beginangle=270, endangle=360)
   circle(1/2, 1/2, r1, :blue)
   circle(x2, x2, r2, :orange)
   circle(1 - x2, 1 - x2, r2, :orange)
   circle(x3, x3, r3, :black)
   circle(1 - x3, 1 - x3, r3, :black)
   y4 = 1 - r4 - 2r5
   ellipse(r4, y4, r4, r4 + 2r5)
   circle(r4, y4, r4, :red)
   circle(r4, y4 - r4 - r5, r5, :gray)
   circle(r4, 1 - r5, r5, :gray)
   ellipse(1 - r4, r4 + 2r5, r4, r4 + 2r5)
   circle(1 - r4, r4 + 2r5, r4, :red)
   circle(1 - r4, r5, r5, :gray)
   circle(1 - r4, 2r4 + 3r5, r5, :gray)

   if more == true
       point(1/2, 1/2, "青")
       point(x2, x2, "黃")
       point(x3, x3, "黒")
       point(r4, y4, "赤")
       point(r4, y4 - r4 - r5, "白")
       point(x0, y0, "(x0,y0)")
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

   r1 = 0.29289; r2 = 0.16019; r3 = 0.03840; r4 = 0.11634; r5 = 0.05817; x2 = 0.17962; x3 = 0.03920; x0 = 0.18350; y0 = 0.57735