算額（その1462）

算額（その1462）

七十五 群馬県吾妻町三島 馬頭観世音 嘉永4年(1851)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円3個，正三角形，矢，第二余弦定理
#Julia, #SymPy, #算額, #和算

算題の文字は判読できないということであるが，「正三角形の中に頂点から 3 ほんの斜線を引き，内部に外側の正三角形と相似である小さな正三角形を作り，隙間にできる相似な三角形に内接する，3 個の円を容れる。外側の正三角形の一辺の長さと円の直径が与えられたとき，内部の正三角形の一辺の長さはいかほどか。」ということでもあろう。

正三角形の一辺の長さを a，円の半径と中心座標を r, (x, y)
OA = a, OB = AC = b, BC = c, AB = c + b
とおき，以下の連立方程式を解く。

以下の get_b_c は，a, r を数値として与え，内部で方程式を解き，b, c を数値として返す関数である。
eq1 は ⊿OAB における第二余弦定理である。a, b, (b + c), ∠ABO = 120°
eq2 は ⊿OAB の面積をヘロンの公式で表したものである。

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

function get_b_c(a, r)
    @syms b::positive, c::positive, s::positive
    s = (a + c + (b + c))/2
    eq1 = a^2 - (c^2 + (b + c)^2 - 2c*(b + c)*cosd(Sym(120)))
    eq2 = sqrt(s*(s - a)*(s - c)*(s - (b + c))) - (a + c + (b + c))*r/2;
    res = solve([eq1, eq2], (b, c))[1]
    return float(res[1]), float(res[2])
end;

正三角形の一辺の長さが 10 寸，円の半径が 1.2 寸のとき，内部の正三角形の一辺の長さ c は 2.01343502545426 寸である。

get_b_c(5, 0.6)

    (1.6659502721190307, 2.01343502545426)

r の取りうる値の範囲は 0 < r < (a/2)\*tan(pi/12) = (a/2)\*(2 - √3) ≒ a*0.13397459621556135 である。
r = 0 のときは内部の正三角形は外部の正三角形と一致し，r = (a/2)\*(2 - √3) のときは内部の正三角形はなくなる。

以下は，図を描くためだけに必要なものである。

@syms a::positive, b::positive, c::positive,
      r::positive, s::positive, S::positive,
      x1::positive, y1::positive,
      x::positive, y::positive
eq3 = x1^2 + y1^2 - c^2
eq4 = (a - x1)^2 + y1^2 - (b + c)^2;
res2 = solve([eq3, eq4], (x1, y1))[1];

eq5 = dist2(0, 0, x1, y1, x, y, r)
eq6 = dist2(a, 0, x1, y1, x, y, r);
res3 = solve([eq5, eq6], (x, y))[3];

function draw(a, r, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (b, c) = get_b_c(a, r)
    @printf("a = %g;  r = %g;  b = %g;  c = %g\n", a, r, b, c)
    (x1, y1) = ((a^2 - b^2 - 2b*c)/(2a), sqrt((b^2 - a^2)*(a - b - 2c)*(a + b + 2c))/(2a))
    t = (a*y1 - r*sqrt(a^2 - 2a*x1 + 2x1^2 + 2y1^2 + 2sqrt(a^2*x1^2 + a^2*y1^2 - 2a*x1^3 - 2a*x1*y1^2 + x1^4 + 2x1^2*y1^2 + y1^4)))
    x = (a^2*r^2*y1 - a^2*y1^3 - 2a*r^2*x1*y1 + 4r^2*x1^2*y1 + 4r^2*y1^3 - 3y1*t^2 + t^3/a + t*(-a^2*r^2 + 3a^2*y1^2 - 4r^2*y1^2)/a)/(2r^2*y1*(-a + 2x1))
    y = t/a
    θ = atand(y1, x1)
    (x2, y2) = (b + c).*[cosd(θ), sind(θ)]
    θ2 = atand(y1, a - x1)
    (x3, y3) = [a - c*cosd(θ2), c*sind(θ2)]
    (Ox, Oy) = (a/2, √3a/6)
    θ3 = atand(√3a/6 - y, a/2 - x) + 180
    l = sqrt((a/2 - x)^2 + (√3a/6 - y)^2)
    (x4, y4) = l.*[cosd(θ3 + 120), sind(θ3 + 120)]
    (x5, y5) = l.*[cosd(θ3 - 120), sind(θ3 - 120)]
    p = plot([0, a, a/2, 0], [0, 0, √3a/2, 0], color=:green, lw=0.5, axis=false)
    circle(x, y, r)
    circle(Ox + x4, Oy + y4, r)
    circle(Ox + x5, Oy + y5, r)
    segment(a, 0, x1, y1, :blue)
    segment(0, 0, x2, y2, :blue)
    segment(a/2, √3a/2, x3, y3, :blue)
    #circle(a/2, √3a/6, √3b/3, :gray70)
    point(0.7a, 0.7a, @sprintf("a = %g\nr = %g", a, r),:black, mark=false)
    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(Ox, Oy)
        point(0, 0, "O:(0,0)", :blue, :center, delta=-delta)
        point(x1, y1, "B:(x1,y1) ", :blue, :right, :vcenter)
        point(x3, y3, "C:(x3,y3)", :blue, :right, delta=-delta)
        point(x, y, "r,(x,y)", :red, :center, delta=-delta)
        point(a, 0, "A:(a,0)", :green, :center, delta=-delta)
        plot!(xlims=(-5delta, a + 25delta), ylims=(-5delta, √3a/2 + 5delta))
    end
    return p
end;

draw(5, 0.6, true)