算額（その1407）

算額（その1407）

百四十三 群馬県榛名町榛名山 榛名神社 明治33年(1900)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：楕円，直角三角形
#Julia, #SymPy, #算額, #和算

直角三角形の中に斜線を 2 本描き，区画された領域に楕円を容れる。鈎が 3 寸，股が 4 寸のとき，楕円の長径，短径はいかほどか。

注：問には明記されておらず，図も正確ではないが，斜線と鈎，股の交点は鈎，股の中点である。

鈎，股，楕円の長径，短径をそのまま変数「鈎」，「股」，「長径」，「短径」とする。

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

三角形に内接する楕円の長径・短径に関する「算法助術」の公式97を援用する。
三角形の底辺, 斜辺1, 斜辺2, 高さを与えて，長径，短径の関係式を返す関数 equation を定義する。

@syms 底辺::positive, 斜辺1::positive, 斜辺2::positive, 高さ::positive, 長径::positive, 短径::positive
equation(底辺, 斜辺1, 斜辺2, 高さ) = -(斜辺1^2 - 斜辺2^2)^2*高さ*短径^2 + (斜辺1^2 - 斜辺2^2)^2*短径^3 + 底辺^4*高さ*(2§高さ - 短径)^2 - 底辺^4*短径*(2高さ - 短径)^2 - 底辺^2*高さ*長径^2*(2高さ - 短径)^2;

同じ楕円が 2 つの直角三角形に内接しているので，以下の連立方程式を解く。

@syms 鈎, 股
eq1 = equation(股, sqrt(股^2 + (鈎/2)^2), 鈎/2, 鈎/2)
eq2 = equation(股/2, sqrt((股/2)^2 + 鈎^2), 鈎, 鈎)
res = solve([eq1, eq2], (長径, 短径))[4];

# 長径
p = res[1] |> simplify |> sympy.sqrtdenest |> simplify
p |> println
p(鈎 => 3, 股 => 4).evalf() |> println

    股*(3 - sqrt(5))/2
    1.52786404500042

# 短径
q = res[2] |> simplify
q |> println
q(鈎 => 3, 股 => 4).evalf() |> println

    鈎*(3 - sqrt(5))/2
    1.14589803375032

長径は股の (3 - √5)/2 倍，短径は鈎の (3 - √5)/2 倍である。
鈎，股が 3 寸，4 寸のとき，長径は 4(3 - √5)/2 = 1.5278640450004204，短径は 3(3 - √5)/2 = 1.1458980337503153 である。

function draw(鈎, 股, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (a, b) = (股, 鈎) .* ((3 - √5)/ 2)  ./ 2
    @printf("鈎，股が %g，%g のとき，楕円の長径，短径は %g，%g である。\n", 鈎, 股, 2a, 2b)
    plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
    ellipse(a, b, a, b, color=:red)
    segment(0, 鈎/2, 股, 0, :green)
    segment(0, 鈎, 股/2, 0, :green)    
    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, 鈎, "鈎 ", :blue, :right, :vcenter)
        point(0, 鈎/2, "鈎/2 ", :blue, :right, :vcenter)
        point(股, 0, "股", :blue, :center, delta=-delta)
        point(股/2, 0, "股/2", :blue, :center, delta=-delta)
        point(a, b, "楕円:a,b,(a,b)", :red, :center, delta=-delta)
        plot!(xlims=(-10delta, 股 + 5delta), ylims=(-7delta, 鈎 + 5delta))
    end
end;

draw(3, 4, true)