算額（その1427）

算額（その1427）

九十七 大船渡市猪川町 雨宝堂 現雨宝山竜宝院 文政7年(1824)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円4個，正三角形，斜線2本
#Julia, #SymPy, #算額, #和算

正三角形の中に斜線 2 本を描き，大円 1 個，中円 1 個，小円 2 個を容れる。小円の直径が 49 寸のとき中円の直径はいかほどか。

正三角形の一辺の長さを 2a
斜線と正三角形の底辺の交点座標を (b, 0)
大円の半径と中心座標を r1, (0, √3a/3); r1 = √3a/3
小円の半径と中心座標を r2, (x2, r2)
中円の半径と中心座標を r3, (0, r3)
とおき，以下の連立方程式を解く

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

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive, x2::positive, r3::positive
r1 = a/√Sym(3)
eq1 = r2/(a - x2) - r1/a
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq3 = dist2(0, √Sym(3)a, b, 0, 0, r3, r3)
eq4 = dist2(0, √Sym(3)a, b, 0, x2, r2, r2)
res = solve([eq1, eq2, eq3, eq4], (a, b, x2, r3))[1];

# a
res[1] |> println

    3*sqrt(3)*r2

# b
res[2] |> println

    11*sqrt(3)*r2/7

# x2
res[3] |> println

    2*sqrt(3)*r2

# r3
res[4] |> println

    99*r2/49

中円の半径 r3 は，小円の半径 r2 の 99/49 倍である。
小円の直径が 49 寸のとき，中円の直径は 99 寸である。

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

    r2 = 24.5;  a = 127.306;  b = 66.684;  r1 = 73.5;  x2 = 84.8705;  r3 = 49.5

function draw(r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (a, b, x2, r3) = (3*sqrt(3)*r2, 11*sqrt(3)*r2/7, 2*sqrt(3)*r2, 99*r2/49)
    r1 = a/√3
    @printf("r2 = %g;  a = %g;  b = %g;  r1 = %g;  x2 = %g;  r3 = %g\n", r2, a, b, r1, x2, r3)
    plot([a, 0, -a, a], [0, √3a, 0, 0], color=:green, lw=0.5)
    circle(0, √3a/3, r1)
    circle2(x2, r2, r2, :blue)
    plot!([b, 0, -b], [0, √3a, 0], color=:magenta, lw=0.5)
    circle(0, r3, r3, :orange)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
        vline!([0], color=:gray80, lw=0.5)
        hline!([0], color=:gray80, lw=0.5)
        point(a, 0, "a", :green, :left, :bottom, delta=delta/2)
        point(b, 0, "b ", :green, :right, :bottom, delta=delta/2)
        point(0, √3a, " √3a", :green, :left, :vcenter)
        point(0, r1, "大円:r1,(0,r1)", :red, :center, delta=-delta)
        point(0, r3, "中円:r3,(0,r3)", :orange, :center, delta=-delta)
        point(x2, r2, "小円:r2\n(x2,r2)", :blue, :center, delta=-delta)
    end
end;

draw(49/2, true)

算額（その1426）

算額（その1426）

八十九 陸前高田市小友町 常膳寺観音堂 天保13年(1842)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円4個，半円2個，四分円1個
#Julia, #SymPy, #算額, #和算

正方形の中に半円 2 個，四分円 1 個，大円 1 個，中円 2 個，小円 1 個を容れる。小円の直径が 1 寸のとき，中円の直径はいかほどか。

正方形の一辺の長さを R
四分円の半径を R
半円の半径を R/2
大円の半径と中心座標を r1, (x1, y1); y1 = x1/2
中円の半径と中心座標を r2, (x21, y21), (x22, y22)
小円の半径と中心座標を r3, (R/2, y3)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms R::positive, r1::positive, x1::positive, y1::positive,
      r2::positive, x21::positive, y21::positive,
      x22::positive, y22::positive, r3::positive, y3::positive
R2 = R/2
y1 = x1/2
eq1 = R2^2 + y3^2 - (R - r3)^2
eq2 = R2^2 + (y3 - R2)^2 - (R2 + r3)^2
eq3 = (R - x1)^2 + (R2 - y1)^2 - (R2 - r1)^2
eq4 = (R - x21)^2 + (R2 - y21)^2 - (R2 - r2)^2
eq5 = (R - x22)^2 + (R2 - y22)^2 - (R2 - r2)^2
eq6 = x21^2 + y21^2 - (R - r2)^2
eq7 = x22^2 + y22^2 - (R - r2)^2
eq8 = (x21 - x1)^2 + (y1 - y21)^2 - (r1 + r2)^2
eq9 = (x22 - x1)^2 + (y1 - y22)^2 - (r1 + r2)^2;

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)
    (R, r1, x1, r2, x21, y21, x22, y22, y3) = u
    return [
        R^2/4 + y3^2 - (R - r3)^2,  # eq1
        R^2/4 + (-R/2 + y3)^2 - (R/2 + r3)^2,  # eq2
        -(R/2 - r1)^2 + (R/2 - x1/2)^2 + (R - x1)^2,  # eq3
        -(R/2 - r2)^2 + (R/2 - y21)^2 + (R - x21)^2,  # eq4
        -(R/2 - r2)^2 + (R/2 - y22)^2 + (R - x22)^2,  # eq5
        x21^2 + y21^2 - (R - r2)^2,  # eq6
        x22^2 + y22^2 - (R - r2)^2,  # eq7
        -(r1 + r2)^2 + (-x1 + x21)^2 + (x1/2 - y21)^2,  # eq8
        -(r1 + r2)^2 + (-x1 + x22)^2 + (x1/2 - y22)^2,  # eq9
    ]
end;

r3 = 1/2
iniv = BigFloat[6.8, 1.3, 4.9, 0.7, 6.1, 0.81, 4.3, 4.4, 5.3]
res = nls(H, ini=iniv)

    ([6.82842712474619, 1.2810685816644012, 4.920484251655344, 0.7126517245944497, 6.059977212374935, 0.8242481001480593, 4.2953848075434085, 4.353432909811113, 5.32842712474619], true)

小円の直径が 1 寸のとき，中円（茶色で示した円）の直径は 2r2 = 2*0.7033666754658107 = 1.4067333509316213 である。

算額の「答」，「術」による中円の直径は 1.775275591337676 寸である。この数値では，図に示した灰色の円になり，明らかに不適切である。

function draw(r3, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (R, r1, x1, r2, x21, y21, x22, y22, y3) = [6.82842712474619, 1.2810685816644012, 4.920484251655344, 0.7126517245944497, 6.059977212374935, 0.8242481001480593, 4.2953848075434085, 4.353432909811113, 5.32842712474619]
    R2 = R/2
    y1 = x1/2
    @printf("小円の直径が %g のとき，中円の直径は %.15g である。\n", 2r3, 2r2)
    #    @printf("r2 = %g;  R = %g;  r1 = %g;  x1 = %g;  x0 = %g;  y0 = %g;  y = %g\n", r2, R, r1, x1, x0, y0, y)
    plot([0, R, R, 0, 0], [0, 0, R, R, 0], color=:green, lw=0.5)
    circle(0, R2, R2, beginangle=-90, endangle=90)
    circle(R, R2, R2, beginangle=90, endangle=270)
    circle(0, 0, R, :blue, beginangle=0, endangle=90)
    circle(x1, y1, r1, :magenta)
    circle(x21, y21, r2, :brown)
    circle(x22, y22, r2, :brown)
    false_r2 = (√45 - 1) * (√3 + 2)/12/2
    println("算額の「答」，「術」による不適切な中円径 = $(2false_r2)")
    circle(x21, y21, false_r2, :gray60, lw=0.3)
    circle(x22, y22, false_r2, :gray60, lw=0.3)
    circle(R2, y3, r3, :orange)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
        vline!([0], color=:gray80, lw=0.5)
        hline!([0], color=:gray80, lw=0.5)
        point(x1, y1, "大円:r1,(x1,y1)", :magenta, :center, delta=-delta)
        point(x21, y21, "中円:r2\n(x21,y21)", :brown, :center, delta=-delta)
        point(x22, y22, "中円:r2\n(x22,y22)", :brown, :center, delta=-delta)
        point(R2, y3, "小円:r3,(R2,y3)", :orange, :right, delta=-delta, deltax=-6delta)
        point(R, 0, " R", :green, :left, :bottom, delta=delta/2)
        point(0, R, " R", :green, :left, :bottom, delta=delta/2)
    end
end;

draw(1/2, true)

算額（その1425）

算額（その1425）

八十九 陸前高田市小友町 常膳寺観音堂 天保13年(1842)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円6個，外円，弦3本
#Julia, #SymPy, #算額, #和算

全円の中に 3 本の弦を引き，区画された領域に甲円 3 個，乙円 2 個を容れる。乙円の直径が 1 寸のとき，甲円の直径はいかほどか。

全円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1), (x1, R - 3r1)
乙円の半径と中心座標を r2, (x2, R - 2r1 + r2)
とおき，以下の連立方程式を解く。

注1：左右にある乙円と甲円の x 座標値は同じになる。x2 = x1
注2：eq1, eq2 は甲円，乙円が全円に内接するという条件を記述したものである。
注3：コメントアウトした eq3, eq4は斜線と甲円，乙円の中心の距離を記述したものであるが，連立方程式を解くことができないほど複雑になる。図に灰色で示した 2 つの円は甲円と乙円と合同なものである。この位置にある甲円，乙円についての式が有効な eq3, eq4 である。これを使えば，連立方程式は容易に解くことができる

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

using SymPy
@syms R::positive, r1::positive, r2::positive, x1::positive, x2, x0::positive
x2 = x1
eq1 = x1^2 + (R - 3r1)^2 - (R - r1)^2
eq2 = x2^2 + (R - 2r1 + r2)^2 - (R - r2)^2
# eq3 = dist2(0, -R, x0, sqrt(R^2 - x0^2), x1, R - 3r1, r1)
# eq4 = dist2(0, -R, x0, sqrt(R^2 - x0^2), x2, R - 2r1 + r2, r2)
eq3 = 2R*r1 - sqrt(x0^2 + (R - sqrt(R^2 - x0^2))^2)*(2R - r1)
eq4 = r2*(2R - r1) - r1*R
res = solve([eq1, eq2, eq3, eq4], (R, r1, x1, x0))[1];

# R
res[1] |> println

    r2*(2 + sqrt(5))

# r1
res[2] |> println

    r2*(1 + sqrt(5))/2

# x1
res[3] |> println

    r2*sqrt(2 + 2*sqrt(5))

# x0
res[4] |> sympy.sqrtdenest |> println

    2*r2*sqrt(-8 + 4*sqrt(5))

甲円の半径 r1 は，乙円の半径 r2 の (1 + √5)/2 倍である。
乙円の直径が 1 寸のとき，甲円の直径は (1 + √5)/2 = 1.618033988749895 寸である。

全てのパラメータは以下のとおりである

   r2 = 0.5;  R = 2.11803;  r1 = 0.809017;  x1 = 1.27202;  x0 = 0.971737;  y0 = 1.88197;  y = 0.5

function draw(r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = r2*(√5 + 2)
    r1 = r2*(√5 + 1)/2
    x1 = r2*sqrt(2√5 + 2)
    x0 = 2r2*sqrt(4√5 - 8)
    x2 = x1
    y0 = sqrt(R^2 - x0^2)
    y = R - 2r1
    @printf("乙円の直径が %g のとき，甲円の直径は %g である。\n", 2r2, 2r1)
    @printf("r2 = %g;  R = %g;  r1 = %g;  x1 = %g;  x0 = %g;  y0 = %g;  y = %g\n", r2, R, r1, x1, x0, y0, y)
    plot([x0, 0, -x0], [y0, -R, y0], color=:green, lw=0.5)
    circle(0, 0, R)
    circle(0, R - r1, r1)
    circle2(x1, y - r1, r1)
    circle2(x2, y + r2, r2, :blue)
    circle(0, y - r2, r2, :gray80)
    circle(0, r1 - R, r1, :gray80)
    segment(sqrt(R^2 - y^2), y, -sqrt(R^2 - y^2), y, :magenta)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
        vline!([0], color=:gray80, lw=0.5)
        hline!([0], color=:gray80, lw=0.5)
        point(0, R, "R", :magenta, :center, :bottom, delta=delta/2)
        point(x1, y + r2, "乙円:r2\n(x2,y+r2)", :blue, :center, delta=-delta/2)
        point(0, R - r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta/2)
        point(x1, y - r1, "甲円:r1,(x1,y-r1)", :red, :center, delta=-delta/2)
        point(x0, sqrt(R^2 - x0^2), "(x0,sqrt(R^2-x0^2) ", :green, :left, :bottom, delta=delta/2)
        point(0, y, " y", :magenta, :left, :bottom, delta=delta/2)
    end
end;

draw(1/2, true)

さぬきのセルフうどん 愉楽家（ゆらくや） 林店

高松市林町 さぬきのセルフうどん 愉楽家（ゆらくや） 林店
ツルツル、やわらか、こしがある
評判通り，麺の量は多め。とり天もでかい

お店の横に，お土産用さぬきうどんの自販機がある。バーコード決済もできる。

算額（その1424）

算額（その1424）

八十九 陸前高田市小友町 常膳寺観音堂 天保13年(1842)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円12個，外円，1/3円3個
#Julia, #SymPy, #算額, #和算

全円の中に正三角形と弦を容れ，区画された領域に 4 個の等円を容れる。等円の等円の直径が 1 寸のとき，全円の直径はいかほどか。

全円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x1, y1), (x2, y2)
弦と y 軸の交点座標を (0, y)
y1 = y + r,  y2 = y - r
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, y::positive, x1::positive, y1::positive, x2::positive, y2::positive
y1 = y + r
y2 = y - r
x0 = R*cosd(Sym(30))
y0 = -R*sind(Sym(30))
eq1 = x1^2 + y1^2 - (R - r)^2 |> expand
eq2 = x2^2 + y2^2 - (R - r)^2 |> expand
eq3 = dist2(0, R, x0, y0, x1, y1, r)
# eq4 = dist2(0, R, x0, y0, x2, y2, r)  # これは eq3 と従属
eq4 = (y1 - y2) - √Sym(3)*(x2 - x1);  # これにすべき

res = solve([eq1, eq2, eq3, eq4], (R, y, x1, x2))[1];

# R
res[1] |> println

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

# y
res[2] |> println

    r*(sqrt(13) + 6)/6

# x1
res[3] |> println

    r*(sqrt(39) + 4*sqrt(3))/6

# x2
res[4] |> factor |> println

    r*(sqrt(39) + 8*sqrt(3))/6

全円の半径 R は，等円の半径の 2(3+√13)/3 倍である。
等円の直径が 1 寸のとき，全円の直径は 2(3+√13)/3 = 4.403700850309327 寸である。

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

   R = 2.20185042515466;  y = 0.800462606288666;  x1 = 1.09776676905616;  x2 = 1.67511703824578

function draw(r, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = 2*r*(3 + √13)/3
    y = r*(√13 + 6)/6
    x1 = r*(√39 + 4√3)/6
    x2 = r*(√39 + 8√3)/6
    x = sqrt(R^2 - y^2)
    y1 = y + r
    y2 = y - r
    x0 = R*cosd(30)
    y0 = -R*sind(30)
    @printf("R = %.15g;  y = %.15g;  x1 = %.15g;  x2 = %.15g\n", R, y, x1, x2)
    plot([-x0, 0, x0, -x0], [y0, R, y0, y0], color=:green, lw=0.5)
    circle(0, 0, R)
    circle2(x1, y1, r, :blue)
    circle2(x2, y2, r, :blue)
    segment(-x, y, x, y, :magenta)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
        vline!([0], color=:gray80, lw=0.5)
        hline!([0], color=:gray80, lw=0.5)
        point(0, R, "R", :magenta, :center, :bottom, delta=delta/2)
        point(x1, y1, "等円:r,(x1,y1)", :blue, :center, delta=-delta/2)
        point(x2, y2, "等円:r,(x2,y2)", :blue, :center, delta=-delta/2)
        point(x0, y0, "(x0,y0) ", :green, :right, :bottom, delta=delta/2)
        point(0, y, " y", :magenta, :left, :bottom, delta=delta/2)
    end
end;

draw(1/2, true)