算額（その1397）

算額（その1397）

三四 武蔵国埼玉郡下忍村遍照院境内 金毘羅社（神楽堂） 天保11年(1840)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円1個，外円，楕円，弦1本
#Julia, #SymPy, #算額, #和算

外円の中に水平な弦を引き，上下に3個の合同な楕円を容れる。楕円は短径の端で外円に接し，長径が最も長いもの，すなわち外円は曲率円である。外円の直径が 5.05 寸のとき，短径はいかほどか。

外円の半径と中心座標を R, (0, 0)
下側の楕円の長半径と短半径，中心座標を a, b, (0, b - R)
楕円と接線の接点座標を　　（x0,y0), (x1, y1)
として以下の連立方程式を解く。

(x0, y0) を通る接線は右上の楕円の中心を通り，(x1, y1) を通る接線は右上の楕円の短径の端を通り直交する。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R, a, b, x0, y0, x1, y1
# a = sqrt(R/b)
eq1 = R - a^2/b
# eq01 = (R - 2b) - √Sym(2)b + sqrt((R - b)*(R - 2b))
eq2 = x0^2/a^2 + (y0 - b + R)^2/b^2 - 1
eq3 = -b^2*x0/(a^2*(y0 - b + R)) - y0/x0
eq4 = x1^2/a^2 + (y1 - b + R)^2/b^2 - 1
eq5 = -b^2*x1/(a^2*(y1 - b + R)) + (-y1 - (R - 2b)*y0/sqrt(x0^2 + y0^2))/(x1 + (R - 2b)*x0/sqrt(x0^2 + y0^2))
eq5 = eq5|> simplify |> numerator
eq6 = x0/y0*(x1 + (R - 2b)*x0/sqrt(x0^2 + y0^2))/(y1 + (R - 2b)*y0/sqrt(x0^2 + y0^2)) + 1
eq6 = eq6 |> simplify |> numerator |> expand |>simplify;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (a, b, x0, y0, x1, y1))

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (a, b, x0, y0, x1, y1) = u
   return [
       R - a^2/b,  # eq1
       -1 + (R - b + y0)^2/b^2 + x0^2/a^2,  # eq2
       -y0/x0 - b^2*x0/(a^2*(R - b + y0)),  # eq3
       -1 + (R - b + y1)^2/b^2 + x1^2/a^2,  # eq4
       a^2*(-y0*(R - 2*b) - y1*sqrt(x0^2 + y0^2))*(R - b + y1) - b^2*x1*(x0*(R - 2*b) + x1*sqrt(x0^2 + y0^2)),  # eq5
       R*x0^2 + R*y0^2 - 2*b*x0^2 - 2*b*y0^2 + x0*x1*sqrt(x0^2 + y0^2) + y0*y1*sqrt(x0^2 + y0^2),  # eq6
   ]
end;

R = 5.05/2
iniv = BigFloat[1.6, 1.1, -1.2, -0.9, 1.4, -1.1]
res = nls(H, ini=iniv)

   ([1.5898411600215787, 1.0010276887519838, -1.198768167169291, -0.86644302003174, 1.44704118366062, -1.1093330662559513], true)

外円の直径が 5.05 のとき，短径は 2.0020553775039676，長径は 3.1796823200431574 である。

全てのパラメータは，以下の通りである。

  a = 1.58984;  b = 1.00103;  x0 = -1.19877;  y0 = -0.866443;  x1 = 1.44704;  y1 = -1.10933

function draw(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho") 
   (a, b, x0, y0, x1, y1) = res[1]
   @printf("外円の直径が %g のとき，短径は %g，長径は %g である。\n", 2R, 2b, 2a)
   @printf("a = %g;  b = %g;  x0 = %g;  y0 = %g;  x1 = %g;  y1 = %g\n", a, b, x0, y0, x1, y1)
   plot()
   circle(0, 0, R)
   ellipse(0, b - R, a, b, color=:blue)
   y = 2b - R
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y, :green)
   θ = atand(y0/x0)
   xo = (R - b)*cosd(θ)
   yo = (R - b)*sind(θ)
   ellipse(xo, yo, a, b, φ=θ + 90, color=:blue)
   ellipse(-xo, yo, a, b, φ=90 - θ, color=:blue)
   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(x0, y0)
       point(x1, y1)
       abline(x0, y0, y0/x0, -R, R, :magenta)
       abline(x1, y1, -x0/y0, -0.5R, 0.9R, :magenta)
       point(xo, yo)
       point(0, b - R, "楕円:a,b,(0,b-R)", :blue, :center, delta=-delta)
       point(0, 2b - R, "2b-R", :blue, :center, delta=-delta)
       point(0, R, "R", :blue, :center, :bottom, delta=delta)
   end  
end;

draw(5.05/2, true)