算額（その1649）
三十一 一関市舞川 観福寺内地蔵堂前額 明治34年(1901)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
今有如図 03085
https://w.atwiki.jp/sangaku/pages/157.html
キーワード：円2個，楕円，正三角形，曲率円
#Julia, #SymPy, #算額, #和算
楕円の中に曲率円を 2 個設け，それらの円と外接するように正三角形を容れる。楕円の長径と短径が与えられたとき，正三角形の一辺の長さを求める術を述べよ。
注：問題文では「設極等円隔斜洩三角」とあるが，山村の図に「斜」はないし，円が曲率円になっていない。たとえ曲率円であっても，山村の図では正三角形の一辺の長さは「楕円の長径から曲率円の直径の 2 倍を引いたもの」という，つまらない解しか得られない。結局，今回も「今有如図」の図に基づいて解く。
楕円の長半径，短半径，中心座標を a, b, (0, 0)
曲率円の半径と中心座標を r, (a - r, 0), (r - a, 0)
正三角形の一辺の長さと頂点の座標を c, (2r - a, y1), (2r - a, y1 + c), (x0, y0)
とおき，以下の連立方程式を解く。
正三角形の頂点 (x0, y0) は楕円の周上にある。
SymPy の能力では解析解を求めることができないので，特定の a, b, における数値解を求める。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
function getc(a, b)
@syms r::positive,
x0::positive, y0::positive, y1::negative,
c::positive;
eq1 = r - b^2/a
eq2 = x0^2/a^2 + y0^2/b^2 - 1
eq3 = dist(2r - a, y1, x0, y0, a - r, 0) - r^2
eq4 = √Sym(3)*(y0 - y1) -(x0 - 2r + a)
eq5 = √Sym(3)*(y0 - y1 - c) + (x0 - 2r + a)
res = solve([eq1, eq2, eq3, eq4, eq5], (r, x0, y0, y1, c))[1]
@printf("楕円の長径，短径が %g，%g のとき，正三角形の一辺の長さは %.15g である。\n", 2a, 2b, res[5])
return res
end;
楕円の長径，短径が 18，10 のとき，正三角形の一辺の長さは 11.5725008614986 である。
山村の「術」の解説は以下のようなっており，11.447328388628943 という，微妙に違う解を与える。
function jutu(長径, 短径)
天 = 長径^2 - 短径^2
地 = sqrt(天 - 短径^2)
人 = sqrt(天)
A = √3長径
B = 短径^2/A
C = 地 + 人 - B
println(C/2)
end;
jutu(18, 10)
11.447328388628943
function draw(a, b, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r, x0, y0, y1, c) = float.(getc(a, b))
plot([x0, 2r - a, 2r - a, x0], [y0, y1 + c, y1, y0], color=:green, lw = 0.5)
circle2(a - r, 0, r)
ellipse(0, 0, a, b, 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(0, b, "b", :blue, :center, :bottom, delta=delta/2, deltax=delta/2)
point(a, 0, "a", :blue, :left, :bottom, delta=delta/2, deltax=delta/2)
point(a - r, 0, "a-r", :red, :center, :bottom, delta=delta/2)
point(2r - a, 0, "2r-a", :blue, :right, :bottom, delta=delta/2, deltax=-delta/2)
point(2r - a, y1, "(2r-a,y1)", :green, :center, delta=-delta/2, deltax=-delta/2)
point(2r - a, y1 + c, "(2r-a,y1+c)", :green, :right, delta=-delta/2, deltax=-delta/2)
point(x0, y0, "(x0,y0)", :blue, :left, :bottom, delta=delta/2)
end
end;
draw(9, 5, true)
