算額（その616）

算額（その616）

福島県白河市明神 境明神 万延元年(1860)
http://www.wasan.jp/fukusima/sakai1L2.html

団扇の中に小円 1 個と長方形が内接し，長方形の中に大円 2 個，小円 7 個が入っている。

小円の径を 1 としたとき，団扇の径はいかほどか。

長方形内の小円の入り方で長方形の長辺，短辺の長さおよび小円の中心座標は簡単に計算できる。

小円の半径を r，外円の半径と中心座標を R, (0, 0)，長方形の右上の頂点座標を (x1, y1) とする。
(x1, y1) は外円の円周上にあるので，等円の半径，外円の半径を r, R として以下の方程式をとき R = r*(4*sqrt(2) + 13)/4 を得る。

include("julia-source.txt")

using SymPy
@syms r::positive, R::positive, x5::positive, y5::negative
x1 = 2√Sym(2)r + r
y1 = R - 2r
eq = x1^2 + y1^2 - R^2
R = solve(eq, R)[1]
R |> println

   r*(4*sqrt(2) + 13)/4

r = 1/2 のとき，R = 2.3321067811865475 で，直径は 4.664213562373095 となる。答えでは 団扇の著系は 5 寸 となっているので，一致しない。

r = 1/2
R = r*(4*sqrt(2) + 13)/4
R, 2R

   (2.3321067811865475, 4.664213562373095)

以上で団扇の主要部を描くことができるが追加で，団扇の下部の竹の骨が見えている円弧部分の描画に必要な座標 (x5, y5) を計算する。

@syms r::positive, R::positive, x5::positive, y5::negative
y2 = R - 3r - √Sym(2)r
eq2 = x5^2 + (R + y5)^2 - (R + y2)^2/4
eq3 = x5^2 + y5^2 - R^2
solve([eq2, eq3], (x5, y5))

   1-element Vector{Tuple{Sym, Sym}}:
    (sqrt(3*R^4/4 - 3*R^3*r/2 - sqrt(2)*R^3*r/2 - 11*R^2*r^2/8 - 3*sqrt(2)*R^2*r^2/4 + 29*sqrt(2)*R*r^3/8 + 45*R*r^3/8 - 193*r^4/64 - 33*sqrt(2)*r^4/16)/R, (-4*R^2 - 12*R*r - 4*sqrt(2)*R*r + 6*sqrt(2)*r^2 + 11*r^2)/(8*R))

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   R = r*(4*sqrt(2) + 13)/4  # R は後々の制約条件から決まる
   println("R = $R, $(1 + sqrt(2))")
   println("2R = $(2R), $(2(1 + sqrt(2)))")
   a = 1 + sqrt(2)
   x1 = 2√2r + r
   y1 = R - 2r
   # (x1, y1) が外円の円周上にあるには x1^2 + y1^2 = R^2 でなければならない。
   y2 = R - 3r - √2r
   x2 = x1 - r - √2r
   y3 = y1 - 2a*r
   plot([x1, -x1, -x1, x1, x1], [y1, y1, y3, y3, y1], color=:blue, lw=0.5)
   circle(0, 0, R, :black)
   circle(0, R - r, r)
   circle(x2, y2 + √2r, r)
   circle(x2, y2 - √2r, r)
   circle(x2 - √2r, y2, r)
   circle(x2 + √2r, y2, r)
   circle(x2, y2, (√2 + 1)r, :green)
   circle(-x2, y2 + √2r, r)
   circle(-x2, y2 - √2r, r)
   circle(-x2 - √2r, y2, r)
   circle(-x2, y2, ( √2 + 1)r, :green)
   (x5, y5) = (sqrt(3*R^4/4 - 3*R^3*r/2 - sqrt(2)*R^3*r/2 - 11*R^2*r^2/8 - 3*sqrt(2)*R^2*r^2/4 + 29*sqrt(2)*R*r^3/8 + 45*R*r^3/8 - 193*r^4/64 - 33*sqrt(2)*r^4/16)/R, (-4*R^2 - 12*R*r - 4*sqrt(2)*R*r + 6*sqrt(2)*r^2 + 11*r^2)/(8*R))
   θ = atand(y5+R, x5)
   println(θ)
   circle(0, -R, (R + y2)/2, :black, beginangle=θ, endangle=181-θ)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(x1, y1, " (x1,y1)")
       point(x1, y2  - (√2 + 1)r, " (x1,y2)")
       point(x2, y2, " (x2,y2)", :green, :left, :vcenter)
       point(0, 0, " O0", :red, :left, delta=-delta/2)
       point(x2, y2 + √2r, " O1", :red, :left, :vcenter)
       point(0, y2, " O2", :red, :left, :vcenter)
       point(x2 + √2r, y2, " O3", :red, :left, :vcenter)
       point(0, R, " R", :black, :left, :bottom, delta=delta/2)
       point(x5, y5, " (x5,y5)", :black, :left, :vcenter)
   end
end;