算額（その1511）

算額（その1511）

百四 岩手県大船渡市田茂山 根城八幡宮 天保12年(1841)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

今有如図 03018
https://w.atwiki.jp/sangaku/pages/206.html

キーワード：円9個，正三角形
#Julia, #SymPy, #算額, #和算

直角三角形の中に斜線を引き，大円 1 個，中円 1 個，小円 2 個を容れる。鈎，股が 3 寸，4 寸のとき，中円の直径はいかほどか。

注：あまりきれいではないが中円，小円の中心は水平線上にある。

正三角形の一辺の長さを 2a
大円の半径と中心座標を r1,(r1, r1)
中円の半径と中心座標を r2, (0, y2), (2r2, y2)
小円の半径と中心座標を r3, (r3, y3), (3r3, y3)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms a, r1, r2, y2, r3, y3
eq1 = (2r2 - r1)^2 + (y2 - r1)^2 - (r1 + r2)^2
eq2 = (3r3 - 2r2)^2 + (y3 - y2)^2 - (r2 + r3)^2
eq3 = dist(a, 0, 0, √Sym(3)*a, r1, r1) - r1^2
eq4 = dist(a, 0, 0, √Sym(3)*a, 2r2, y2) - r2^2
eq5 = dist(a, 0, 0, √Sym(3)*a, 3r3, y3) - r3^2;
# res = solve([eq1, eq2, eq3, eq4, eq5], (a, r1, r2, y2, y3))

function H(u)
    (a, r1, r2, y2, y3) = u
    return [
        (2r2 - r1)^2 + (y2 - r1)^2 - (r1 + r2)^2,
        (3r3 - 2r2)^2 + (y3 - y2)^2 - (r2 + r3)^2,
        dist(a, 0, 0, √3*a, r1, r1) - r1^2,
        dist(a, 0, 0, √3*a, 2r2, y2) - r2^2,
        dist(a, 0, 0, √3*a, 3r3, y3) - r3^2
    ]
end;
r3 = 4/2
iniv = BigFloat[23.5, 8.6, 3.6, 20.7, 26.2]
res = nls(H, ini=iniv)

    ([23.453792062852013, 8.58468371008167, 3.6391447593398145, 20.738502725786578, 26.23085463760209], true)

大円の直径は res[1][2]*2 = 17.16936742016334 である。

術は小円の直径の 4.25 倍とのみ書いているが，不明。4.25*4 = 17

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r3 = 4/2  # 0.08933
    (a, r1, r2, y2, y3) = [1, 0.36, 0.15239, 0.88, 1.12]
    (a, r1, r2, y2, y3) = res[1]
    plot([a, 0, -a, a], [0, √3a, 0, 0], color=:magenta, lw=0.5)
    circle2(r1, r1, r1)
    circle(0, y2, r2, :blue)
    circle2(2r2, y2, r2, :blue)
    circle2(r3, y3, r3, :green)
    circle2(3r3, y3, r3, :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(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta)
        point(2r2, y2, "中円:r2,(2r2,y2)", :blue, :center, delta=-delta)
        point(3r3, y3, "小円:r3,(3r3,y3)", :black, :center, delta=-delta)
        point(a, 0, "a", :magenta, :left, :bottom, deltax=delta/2)
        point(0, √3a, "√3a", :magenta, :left, :bottom, deltax=delta/2)
    end
end;

draw(true)

算額（その1510）

算額（その1510）

九十九 岩手県江刺市 雨宝堂 現在中善観音堂 文政10年(1827)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

今有如図 03007
https://w.atwiki.jp/sangaku/pages/205.html

キーワード：円4個，直角三角形，斜線1本
#Julia, #SymPy, #算額, #和算

直角三角形の中に斜線を引き，大円 1 個，中円 1 個，小円 2 個を容れる。鈎，股が 3 寸，4 寸のとき，中円の直径はいかほどか。

鈎，股をそのまま「鈎」，「股」
大円の半径と中心座標を r1, (x1, r1)
中円の半径と中心座標を r2, (股 - r2, y2)
小円の半径と中心座標を r3, (x3, r3), (x32, y32)
斜線と斜辺の交点座標を (x0, y0)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms 鈎, 股, r1, x1, r2, y2, r3, x3, x32, y32, x0, y0
eq1 = dist(0, 0, 股, 鈎, x1, r1) - r1^2
eq2 = dist(0, 0, 股, 鈎, 股 - r2, y2) - r2^2
eq3 = dist(0, 0, 股, 鈎, x32, y32) - r3^2
eq4 = dist(x0, y0, 股, 0, x1, r1) - r1^2
eq6 = dist(x0, y0, 股, 0, x32, y32) - r3^2
eq7 = dist(x0, y0, 股, 0, 股 - r2, y2) - r2^2
eq8 = x0^2 + y0^2 - 股^2
eq9 = (2股 + sqrt((股 - x0)^2 + y0^2))*r1 - 股*y0
eq11 = r3/(股 - x3) - r1/(股 - x1)
eq12 = (x3 - x1)^2 + (r1 - r3)^2 - (r1 + r3)^2;

function H(u)
    (r1, x1, r2, y2, r3, x3, x32, y32, x0, y0) = u
    return [
        dist(0, 0, 股, 鈎, x1, r1) - r1^2,
        dist(0, 0, 股, 鈎, 股 - r2, y2) - r2^2,
        dist(0, 0, 股, 鈎, x32, y32) - r3^2,
        dist(x0, y0, 股, 0, x1, r1) - r1^2,
        dist(x0, y0, 股, 0, x32, y32) - r3^2,
        dist(x0, y0, 股, 0, 股 - r2, y2) - r2^2,
        x0^2 + y0^2 - 股^2,
        (2股 + sqrt((股 - x0)^2 + y0^2))*r1 - 股*y0,
        r3/(股 - x3) - r1/(股 - x1),
        (x3 - x1)^2 + (r1 - r3)^2 - (r1 + r3)^2
    ]
end;
(鈎, 股) = (3, 4)
iniv = BigFloat[0.91753, 2.72942, 0.36849, 2.26201, 0.2366, 3.6581, 3.07768, 2.00873, 3.2, 2.4]
res = nls(H, ini=iniv)

    ([0.9116963119775494, 2.7350889359326485, 0.3675444679663241, 2.2649110640673515, 0.2389194451114356, 3.668516977010949, 3.0781652486756204, 2.0099746301174206, 3.2, 2.4], true)

中円の直径は res[1][3]*2 = 0.7350889359326482 である。

function draw(r, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (鈎, 股) = (3, 4)
    (r1, x1, r2, y2, r3, x3, x32, y32, x0, y0) = [0.91753, 2.72942, 0.36849, 2.26201, 0.2366, 3.6581, 3.07768, 2.00873, 3.2, 2.4]
    (r1, x1, r2, y2, r3, x3, x32, y32, x0, y0) = res[1]
    plot([0, 股, 股, 0], [0, 0, 鈎, 0], color=:magenta, lw=0.5)
    circle(x1, r1, r1)
    circle(股 - r2, y2, r2, :blue)
    circle(x3, r3, r3, :green)
    circle(x32, y32, r3, :green)
    segment(股, 0, x0, y0)
    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(x1, r1, "大円:r1,(x1,r1)", :red, :center, delta=-delta)
        point(股 - r2, y2, "中円:r2\n(股-r2,y2)", :blue, :center, :bottom, delta=delta)
        point(x3, r3, "小円:r3,(x3,r3)", :green, :right, :vcenter, deltax=-8delta)
        point(x32, y32, "小円:r3,(x32,y32)", :green, :right, :bottom, delta=2delta, deltax=-6delta)
        point(股, 0, " 股", :magenta, :left, :bottom, delta=delta, deltax=-delta/2)
        point(股, 鈎, "(股,鈎)", :magenta, :left, :bottom, delta=delta, deltax=-delta/2)
    end
end;

draw(2, true)

算額（その1509）

算額（その1509）

六十八 岩手県一関市川崎町薄衣諏訪前 浪分神社 文政5年(1822)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

今有如図 03004
https://w.atwiki.jp/sangaku/pages/202.html

キーワード：円6個，外円，弦4本
#Julia, #SymPy, #算額, #和算

外円の中に弦を 4 本引き，区画された領域に等円 5 個を容れる。等円の直径が与えられたとき，外円の直径を求める術を述べよ。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (0, R - r), (x1, y1), (x2, y2)
弦の端点座標を (x01, sqrt(R^2 - x01^2)), (x02, sqrt(R^2 - x02^2)), (0, -R)
とおき，以下の連立方程式を解く。

R, r1 を記号のままにして SymPy では簡単に解くことができない。数値を代入して解けば，数値解は求まる。

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

using SymPy
@syms R::positive, r::positive, x1::positive, y1::positive,
      x2::positive, y2::negative,  x01::positive, x02::positive
eq1 = x1^2 + y1^2 - (R - r)^2
eq2 = x2^2 + y2^2 - (R - r)^2
eq3 = dist(0, -R, x01, sqrt(R^2 - x01^2), 0, R - r) - r^2
eq4 = dist(0, -R, x01, sqrt(R^2 - x01^2), x1, y1) - r^2
eq5 = dist(0, -R, x02, sqrt(R^2 - x02^2), x1, y1) - r^2
eq6 = dist(0, -R, x02, sqrt(R^2 - x02^2), x2, y2) - r^2
eq7 = (sqrt(R^2 - x02^2) + R)/x02 * y2/x2 + 1;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (R, x1, y1, x2, y2, x01, x02))  

function H(u)
    (R, x1, y1, x2, y2, x01, x02) = u
    return [
        x1^2 + y1^2 - (R - r)^2,
        x2^2 + y2^2 - (R - r)^2,
        dist(0, -R, x01, sqrt(R^2 - x01^2), 0, R - r) - r^2,
        dist(0, -R, x01, sqrt(R^2 - x01^2), x1, y1) - r^2,
        dist(0, -R, x02, sqrt(R^2 - x02^2), x1, y1) - r^2,
        dist(0, -R, x02, sqrt(R^2 - x02^2), x2, y2) - r^2,
        (sqrt(R^2 - x02^2) + R)/x02 * y2/x2 + 1
    ]
end;
r = 2
iniv = BigFloat[10, 5, 5, 7, -3, 3, 8]
res = nls(H, ini=iniv)

    ([7.509090269827044, 3.8101684797908897, 3.9790315098916964, 4.870544338179145, -2.5744656631880116, 2.2798713379711115, 6.204719439604394], true)

等円の直径が 4 のとき，外円の直径は 2*7.509090269827044 = 15.018180539654088 である。

function draw(r, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (R, x1, y1, x2, y2, x01, x02) = res[1]
    plot()
    circle(0, 0, R, :green)
    circle(0, R - r, r)
    circle2(x1, y1, r)
    circle2(x2, y2, r)
    plot!([-x01, 0, x01], [sqrt(R^2 - x01^2), -R, sqrt(R^2 - x01^2)], color=:blue, lw=0.5)
    plot!([-x02, 0, x02], [sqrt(R^2 - x02^2), -R, sqrt(R^2 - x02^2)], color=:blue, lw=0.5)
    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, R, "R", :green, :center, :bottom, delta=delta)
        point(x1, y1, "r,(x1,y1)", :red, :center, delta=-delta)
        point(x2, y2, "r,(x2,y2)", :red, :center, delta=-delta)
        point(x01, sqrt(R^2 - x01^2), "(x01,y01)", :blue, :left, :bottom, delta=delta)
        point(x02, sqrt(R^2 - x02^2), "(x02,y02)", :blue, :left, :bottom, delta=delta)
        point(0, R - r, " R-r", :blue, :left, :vcenter)
    end
end;

draw(2, true)

術は，「外円の直径を x，等円の直径を d として 以下の d の 5 次式を解くとしている。
-d^5 + 4*d^4*x - 14*d^3*x^2 + 52*d^2*x^3 - 73*d*x^4 + 16*x^5

using SymPy
@syms x, d
A = ((16x - 73d)*x + 52d^2)*x - 14d^3
eq = (A*x + 4d^4)*x - d^5 |> expand
eq |> println

    -d^5 + 4*d^4*x - 14*d^3*x^2 + 52*d^2*x^3 - 73*d*x^4 + 16*x^5

d = 4 のときの解を求める。5 個の解のうち，3 番目のものが適解 15.0181805396541 である。

res = solve(eq(d => 4));
res[3].evalf() |> println

    15.0181805396541

 

算額（その1508）

算額（その1508）

茨城県水戸市文京 笠原神社 明治10年(1877)
https://w.atwiki.jp/sangaku/pages/201.html

キーワード：円6個，外円，弦
#Julia, #SymPy, #算額, #和算

外円の中に弦を引き，甲円 1 個，乙円 4 個を容れる。外円の直径が 12.5 寸，甲円の直径が 10.5 寸のとき，乙円の直径はいかほどか。

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

R, r1 を記号のままにして SymPy では簡単に解くことができない。数値を代入して解けば，数値解は求まる。

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, x2::positive, y2::negative, x22::positive
R = 12.5/2
r1 = 10.5/2
eq1 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + y2^2 - (R - r2)^2
eq3 = x22^2 + (r1 - r2)^2 - (r1 + r2)^2
eq4 = (x2 - x22)^2 + (y2 - (R - 2r1 + r2))^2 - 4r2^2;
res = solve([eq1, eq2, eq3, eq4], (r2, x2, y2, x22))[1]  # 1 of 2

    (0.750000000000000, 4.96078370824611, -2.37500000000000, 3.96862696659689)

外円の直径が 12.5 寸，甲円の直径が 10.5 寸のとき，乙円の直径は 1.5 寸である。

function draw(R, r1, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r2, x2, y2, x22) = (0.750000000000000, 4.96078370824611, -2.37500000000000, 3.96862696659689)
    y = R - 2r1
    x = sqrt(R^2 - y^2)
    plot()
    circle(0, 0, R, :green)
    circle(0, R - r1, r1)
    circle2(x2, y2, r2, :blue)
    circle2(x22, y + r2, r2, :blue)
    segment(-x, y, x, y)
    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, 0, "外円:R,(0,0)", :green, :center, delta=-delta)
        point(0, R - r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta)
        point(x2, y2, "乙円:r2,(x2,y2)", :blue, :center, delta=-delta)
        point(x22, y + r2, "乙円:r2,(x22,y+r2)", :blue, :center, delta=-delta)
        point(0, R, "R", :green, :center, :bottom, delta=delta)
        point(0, y, " y=R-2r", :black, :center, :bottom, delta=delta)
    end
end;

draw(12.5/2, 10.5/2, true)

古武そば・さぬきうどん 阿讃亭

高松市塩江町 古武そば・さぬきうどん 阿讃亭
前に言及した，もうすぐ徳島県という手前の 2 番目の阿讃亭。
どちらかというとそばが売りの一般店。ちょっと単価は高めである。常連さんが多い。営業日に注意。
最近のハマりで，しっぽくうどんを頼んだ。かなりの細麺で，コシはある。

算額（その1507）

算額（その1507）

八十九 陸前高田市小友町 常膳寺観音堂 天保12年(1841)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03019
https://w.atwiki.jp/sangaku/pages/208.html

キーワード：円2個，直角三角形，正方形，斜線1本
#Julia, #SymPy, #算額, #和算

直角三角形の中に甲円，乙円，正方形，斜線を容れる。鈎，股が与えられたとき乙円の直径を求める術を述べよ。

鈎，股をそのまま，「鈎」，「股」
正方形の一辺の長さを a
斜線と正方形の一辺の交点座標を (b, 0)
甲円の半径と中心座標を r1, (a + r1, r1)
乙円の半径と中心座標を r2, (a - r2, a - r2)
とおき，以下の方程式を解く。

一度に解けないので，逐次求めていく。

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

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive, r1::positive, r2::positive, a::positive, b::positive
eq1 = (股 - a) + a - sqrt((股 - a)^2 + a^2) - 2r1
eq2 = dist2(0, 鈎, b, 0, a + r1, r1, r1)
eq3 = dist2(0, 鈎, b, 0, a - r2, a - r2, r2)
eq4 =  a/(股 - a) - 鈎/股;
# res = solve([eq1, eq4, eq5, eq6], (r2, r1, a, b))

# a
ans_a = solve(eq4, a)[1]
ans_a |> println

    股*鈎/(股 + 鈎)

# r1
eq11 = eq1(a => ans_a)
ans_r1 = solve(eq11, r1)[1]#(sqrt(鈎^2 + 股^2) => 弦)
ans_r1 |> println

    股*(股 + 鈎 - sqrt(股^2 + 鈎^2))/(2*(股 + 鈎))

# b
eq12 = eq2(a => ans_a, r1 => ans_r1) |> simplify
ans_b = solve(eq12, b)[2](sqrt(鈎^2 + 股^2) => 弦) |> simplify
ans_b |> println

    股*鈎^2*(-弦 + 股 + 2*鈎)/(弦*股^2 + 弦*股*鈎 - 股^3 - 股^2*鈎 + 股*鈎^2 + 鈎^3)

# r2
eq13 = eq3(a => ans_a, r1 => ans_r1, b => ans_b)
ans_r2 = solve(eq13, r2)[2];  # 1 of 2
ans_r2 |> println

    ((股 + 鈎)^3*(-弦 + 股 + 2*鈎)*(弦*股 - 股^2 + 鈎^2)*sqrt(弦^2*股^4 + 2*弦^2*股^3*鈎 + 2*弦^2*股^2*鈎^2 - 2*弦*股^5 - 4*弦*股^4*鈎 - 2*弦*股^3*鈎^2 + 2*弦*股*鈎^4 + 股^6 + 2*股^5*鈎 + 3*股^2*鈎^4 + 2*股*鈎^5 + 鈎^6) + (弦*股^2 - 股^3 + 3*股*鈎^2 + 鈎^3)*(弦^2*股^4 + 3*弦^2*股^3*鈎 + 3*弦^2*股^2*鈎^2 + 弦^2*股*鈎^3 - 2*弦*股^5 - 8*弦*股^4*鈎 - 11*弦*股^3*鈎^2 - 5*弦*股^2*鈎^3 + 弦*股*鈎^4 + 弦*鈎^5 + 股^6 + 5*股^5*鈎 + 8*股^4*鈎^2 + 2*股^3*鈎^3 - 7*股^2*鈎^4 - 7*股*鈎^5 - 2*鈎^6))*(弦*股^2 + 弦*股*鈎 + 弦*鈎^2 - 股^3 - 股^2*鈎 - 鈎^3)/(2*(股 + 鈎)^2*(-弦 + 股 + 2*鈎)*(弦*股 - 股^2 + 鈎^2)*(弦^2*股^4 + 3*弦^2*股^3*鈎 + 3*弦^2*股^2*鈎^2 + 弦^2*股*鈎^3 - 2*弦*股^5 - 8*弦*股^4*鈎 - 11*弦*股^3*鈎^2 - 5*弦*股^2*鈎^3 + 弦*股*鈎^4 + 弦*鈎^5 + 股^6 + 5*股^5*鈎 + 8*股^4*鈎^2 + 2*股^3*鈎^3 - 7*股^2*鈎^4 - 7*股*鈎^5 - 2*鈎^6))

計算式はかなり長く，SymPy でも簡約化しきれない。
弦を前もって計算しておき，以下により小円の半径を求める。

鈎 = 3, 股 = 4, 弦 = 5 のとき，小円の半径は 0.303296703296703 である。

ans_r2(弦 => sqrt(鈎^2 + 股^2))(鈎 => 3, 股 => 4) |> println
ans_r2(弦 => sqrt(鈎^2 + 股^2))(鈎 => 3, 股 => 4).evalf() |> println

    138/455
    0.303296703296703

function draw(鈎, 股, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    弦 = sqrt(鈎^2 + 股^2)
    term1 = 弦^2*股^4 + 3*弦^2*股^3*鈎 + 3*弦^2*股^2*鈎^2 + 弦^2*股*鈎^3 - 2*弦*股^5 - 8*弦*股^4*鈎 - 11*弦*股^3*鈎^2 - 5*弦*股^2*鈎^3 + 弦*股*鈎^4 + 弦*鈎^5 + 股^6 + 5*股^5*鈎 + 8*股^4*鈎^2 + 2*股^3*鈎^3 - 7*股^2*鈎^4 - 7*股*鈎^5 - 2*鈎^6
    term2 = 弦*股 - 股^2 + 鈎^2
    term3 = -弦 + 股 + 2*鈎
    term4 = 股 + 鈎
    num = (term4^3*term3*term2*sqrt(弦^2*股^4 + 2*弦^2*股^3*鈎 + 2*弦^2*股^2*鈎^2 - 2*弦*股^5 - 4*弦*股^4*鈎 - 2*弦*股^3*鈎^2 + 2*弦*股*鈎^4 + 股^6 + 2*股^5*鈎 + 3*股^2*鈎^4 + 2*股*鈎^5 + 鈎^6) + (弦*股^2 - 股^3 + 3*股*鈎^2 + 鈎^3)*term1)*(弦*股^2 + 弦*股*鈎 + 弦*鈎^2 - 股^3 - 股^2*鈎 - 鈎^3)
    den = 2*term4^2*term3*term2*term1
    r1 = 股*(-弦 + 股 + 鈎)/(2*(股 + 鈎))
    b = 股*鈎^2*(-弦 + 股 + 2*鈎)/(弦*股^2 + 弦*股*鈎 - 股^3 - 股^2*鈎 + 股*鈎^2 + 鈎^3)
    a = 股*鈎/(股 + 鈎)
    r2 = num/den
    println(r2)
    plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:magenta, lw=0.5)
    plot!([a, a, 0], [0, a, a], color=:magenta, lw=0.5)
    segment(0, 鈎, b, 0)
    circle(a + r1, r1, r1, :magenta)
    circle(a - r2, a - r2, r2, :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(a + r1, r1, "甲円:r1,(a+r1,r1)", :red, :center, delta=-delta)
        point(a - r2, a - r2, "乙円:r2,(a-r2,a-r2)", :red, :center, delta=-delta)
        point(a, 0, "a", :black, :center, delta=-delta)
        point(b, 0, "b", :black, :center, delta=-delta)
        point(0, 鈎, "鈎", :magenta,  :left, :bottom, delta=delta)
        point(股, 0, "股", :magenta,  :left, :bottom, delta=delta)
        ylims!(-5delta, 鈎 +5delta)
    end
end;

draw(3, 4, true)

 

算額（その1506）

算額（その1506）

八十八 陸前高田市小友町 常膳寺観音堂 天保12年(1841)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03019
https://w.atwiki.jp/sangaku/pages/208.html

キーワード：円3個，半円3個，正三角形
#Julia, #SymPy, #算額, #和算

この算額は，「問」，「答」，「術」とも読み取ることができないそうだ。図形のみ認識できるようなので，いくつかは再現できそうである。

ただし，山村の図は不正確で，円弧は半円であるという条件がなければ解けない。「今有如図」を見て算額（その1052）から大円を取り除いたものと同じであることに気づいた。

正三角形の中に，半円，円を 3 個ずつ容れる。円の直径が 1 寸のとき，正三角形の一辺の長さはいかほどか。

半円の半径は正三角形の一辺の長さの √3/2 倍で，その中心は正三角形の辺の中点である。
正三角形の一辺の長さを 2a
半円の半径と中心座標を R, (0, 0), (a/2, √3a/2), (-a/2, √3a/2)
円の半径と中心座標を r, (0, r)
とおき，以下の方程式を解く。

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

using SymPy
@syms a::positive, R::positive, r1::positive, x1::positive, r::positive
R = √Sym(3)a/2
eq3 = (a/2)^2 + (√Sym(3)a/2 - r)^2 - (R + r)^2

res = solve(eq3, a)[1];
res |> println

    8*sqrt(3)*r

正三角形の一辺の長さは，円の半径の 8√3 倍である。
円の直径が 1 寸のとき，正三角形の一辺の長さは 13.856406460551018 寸である。

function draw(r, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    a = 8√3r
    R = √3a/2
    @printf("r = %g;  a = %g;  R = %g\n",  r, a, R)
    plot([a, 0, -a, a], [0, √3a, 0, 0], color=:magenta, lw=0.5)
    circle(0, 0, R, :blue,  beginangle=0, endangle=180)
    circle(a/2, √3a/2, R, :blue, beginangle=120, endangle= 300)
    circle(-a/2, √3a/2, R, :blue, beginangle=-120, endangle= 60)
    circle(0, r, r, :green)
    circle2(a/2 - r*cosd(30), √3a/2 - r*sind(30), r, :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, r, " 円:r,(0,r)", :black, :left, :vcenter)
        point(a/2, √3a/2, "半円:R\n(a/2,√3a/2)", :blue, :left, :vcenter, deltax=2delta)
        point(0, R, " R", :blue, :center, :bottom, delta=delta/2)
        point(0, √3a, " √3a", :blue, :center, :bottom, delta=delta/2)
        point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
        point(a, 0, " a", :magenta, :left, :bottom, delta=delta/2)
    end
end;

draw(1/2, true)

算額（その1505）

算額（その1505）

七十一 岩手県一関市川崎町薄衣諏訪前 浪分神社 明治35年(1902)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03054
https://w.atwiki.jp/sangaku/pages/138.html

キーワード：円10個，長方形
#Julia, #SymPy, #算額, #和算

外円と長方形が交差している。区画された領域に甲円 2 個，乙円 2 個，丙円 5 個を容れる。
(1) 甲円径，乙円径，丙円径の差はそれぞれ 1 寸
(2) 「外円径の 3 乗と丙円径の積」と「甲円径と乙円径の積」の差は 337
(3) 外円径と丙円径の差は 6 寸
のとき，外円の直径はいかほどか。



条件が過剰なので，(2) が無くても解ける。

外円，甲円，乙円，丙円の直径を D, d1, d2, d3 とおき，以下の連立方程式を解く。

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

using SymPy

@syms D, d1, d2, d3
eq1 = d1 - d2 - 1
eq2 = d2 - d3 - 1
# eq2 = D^3*d3 - d1*d2 - 337
eq3 = D - d3 - 6
eq4 = (2d1 + d3) - (2d2 + 3d3)
res = solve([eq1, eq2, eq3, eq4], (D, d1, d2, d3));

res |> println

    Dict{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}(d1 => 3, d2 => 2, D => 7, d3 => 1)

外円の直径は 7 寸である。

function draw(r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (R, r1, r2, r3) = (7, 3, 2, 1)./2
    a = 2r1 + 3r3
    b = 2r2 + r3
    plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:magenta, lw=0.5)
    circle(0, 0, R, :blue)
    circle(0, 0, r3)
    circle2(2r3 + 2r1, 0, r3)
    circle22(0, 2r3 + 2r2, r3)
    circle2(r3 + r1, 0, r1, :orange)
    circle22(0, r3 + r2, r2, :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(r3 + r1, 0, "甲円", :orange, :center, :vcenter, mark=false)
        point(0, r3 + r2, "乙円", :green, :center, :vcenter, mark=false)
        point(0, 0, "丙円", :red, :center, :vcenter, mark=false)
    end
end;

draw(1/2, true)

算額（その1042）改訂版

算額（その1042）改訂版

算額（その1042）は山村の描いた不適切な図に基づいていたので，改訂版を書いた。

八十六 室根村室根山 室根神社 明治32年(1899)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03082
https://w.atwiki.jp/sangaku/pages/233.html

キーワード：円12個，正方形
#Julia, #SymPy, #算額, #和算

2 個の正方形が交差している。区分された領域に甲円，乙円を 2 個ずつと丙円を 8 個描く。甲円の直径が 3 寸のとき，乙円の径はいかほどか。

算額（その1042）で，『甲円の直径が 3 寸のとき，「答曰乙円径二寸」とあるが，算額の図を見るだけでそんな解はない。そんな図は描けない。山村も，術を鸚鵡返ししているだけで，答，術がおかしいことを指摘していない。』と書いたが，そもそも山村の描いた図のほうがおかしいので，答，術がおかしいのではない。

甲円の半径と中心座標を r1, (r1, 0)
乙円の半径と中心座標を r2, (0, 2r1 - r2)
丙円の半径と中心座標を r3, (x3, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive, x3::positive
eq1 = r1^2 + (2r1 - r2)^2 - (r1 + r2)^2
eq2 = dist2(2√Sym(2)*r1, 0, 0, 2√Sym(2)*r1, 2r1 + r3, 0, r3)
res = solve([eq1, eq2], (r2, r3))[1];

# r2
res[1] |> println

    2*r1/3

# r3
res[2] |> simplify |> println

    2*r1*(3 - 2*sqrt(2))

乙円の直径は甲円の直径の 2/3 倍である。
答，術ともに正しい。

function draw(r1, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r2 = 2r1/3
    r3 = 2r1*(3 - 2√2)
    @printf("甲円の直径が %g のとき，乙円の直径は %g である。\n", 2r1, 2r2)
    @printf("r1 = %g;  r2 = %g;  r3 = %g\n", r1, r2, r3)
    a = 2r1
    x = y = 2√2r1
    plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:blue, lw=0.5)
    plot!([x, 0, -x, 0, x], [0, y, 0, -y, 0], color=:orange, lw=0.5)
    circle2(r1, 0, r1)
    circle22(0, a - r2, r2, :magenta)
    circle4(a - r3, a - r3, r3, :green)
    circle22(0, a + r3, r3, :green)
    circle2(a + r3, 0, r3, :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(r1, 0, "甲円:r1,(r1,0)", :red, :center, delta=-delta/2)
        point(0, a - r2, " 乙円:r2\n(0,a-r2)", :magenta, :center, delta=-delta/2)
        point(0, a + r3, " 丙円:r3", :black, :left, :vcenter)
        point(a, a, "(2r1,2r1)", :blue, :right, :bottom, delta=delta/2)
        point(x, 0, " x", :orange, :left, :bottom, delta=delta/2)
        point(0, y, " y", :orange, :left, :bottom, delta=delta/2)
    end
end;

draw(3/2, true)

算額（その1504）

算額（その1504）

参考文献

1).  七十一 岩手県一関市川崎町薄衣諏訪前 浪分神社 明治35年(1902)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

2).  今有如図 03054
https://w.atwiki.jp/sangaku/pages/138.html

キーワード：球13個，円錐
#Julia, #SymPy, #算額, #和算

円錐の中に大球 1 個，小球 12 個を容れる。大球は円錐の側面に内接し，小球と外接する。小球は円錐の側面に内接し，大球に外接し，互いに外接し合っている。
小球の直径が 1 寸のとき，円錐の底面の直径はいかほどか。

算額（その1495）の類題である。

左図は横から見た図，右図は上から見た図 青は大球，赤は小球，オレンジ色は円錐の外縁（右の図では底面の円）。小球の中心は灰色の円周上にある。

円錐の底面の半径と座標を a, (0, 0, 0)
円錐の高さを h
大球の半径と中心座標を r1, (0, 0, r1)
小球の半径と中心座標を r2, (b, 0, r2)
とおき，以下の連立方程式を r2 を既知として解き (a, r1, h) を求める。

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

using SymPy

@syms a, b, h, r1, r2
b = r2/sind(Sym(30)/2)
eq1 = r1/(h - r1) - a/sqrt(a^2 + h^2)
eq2 = 2sqrt(r1*r2) - b
eq3 = r2/(a - b) - r1/a;
res = solve([eq1, eq2, eq3], (a, r1, h))[1]

    (r2*(sqrt(6) + 2*sqrt(2)), r2*(sqrt(3) + 2), 4*r2*(sqrt(3) + 2))

底面の円の半径 a は，小球の半径 r2 の (√6 + 2√2) 倍である。

小球の直径が 1 寸のとき，底面の円の直径は 5.27791686752936 寸である。

ちなみに円錐の高さは 7.46410161513775 寸，大球の直径は 3.73205080756888 寸である。

res[2](r2 => 1).evalf() |> println

    3.73205080756888

術は sqrt(√192 + 14) であるが，二重根号を外すと上と同じく √6 + 2√2 である。

function draw(r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (a, r1, h) = (r2*(sqrt(6) + 2*sqrt(2)), r2*(sqrt(3) + 2), 4*r2*(sqrt(3) + 2))
    b = r2/sind(30/2)
    @printf("小球の直径が %g のとき，底面の円の直径は %g， 高さは %g， 大球の直径は %g である。\n", 2r2, 2a, h, 2r1)
    p1 = plot([a, 0, -a, a], [0, h, 0, 0], color=:green, lw=0.5)
    circle(0, r1, r1)
    circle(b, r2, r2, :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)
    end
    p2 = plot()
    rotate(b, 0, r2, angle=30, :blue)
    circle(0, 0, b, :gray80)
    circle(0, 0, r1)
    circle(0, 0, a, :green)
    plot(p1, p2)
end;

draw(1/2, true)

算額（その1503）

算額（その1503）

六十三 花泉町金沢 大門神社・大門観世音菩薩堂 明治13年(1880)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03064
https://w.atwiki.jp/sangaku/pages/292.html

キーワード：円2個，円弧，斜線
#Julia, #SymPy, #算額, #和算

大円の一部である円弧（弓形）の中に小円を 2 個と界斜を容れる。大円の径が 8 寸，矢が 2 寸，小円の直径が 1 寸のとき，界斜はいかほどか。

山村の図は右側の小円と界斜の端点の位置が違う。

「今有如図」の仮定で解を導き，答と一致することを確認した。

大円の半径と中心座標を R, (0, 0)
矢をそのまま「矢」
小円の半径と中心座標を r, (x, R - 矢 + r)
界斜の端点を (x1, R - 矢), (x2, sqrt(R^2 - x2^2))
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, 矢::positive, x::positive, r::positive,
      x1::positive, x2::negative
y = R - 矢
eq1 = x^2 + (y + r)^2 - (R - r)^2
eq2 = dist(x2, sqrt(R^2 - x2^2), x1, y, x, y + r) - r^2
eq3 = (sqrt(R^2 - x2^2) - y)/(x1 - x2) - r/x1;
# res = solve([eq1, eq2, eq3], (x, x1, x2))

function H(u)
    (x, x1, x2) = u
    return [
        x^2 + (y + r)^2 - (R - r)^2,
        dist(x2, sqrt(R^2 - x2^2), x1, y, x, y + r) - r^2,
        (sqrt(R^2 - x2^2) - y)/(x1 - x2) - r/x1
    ]
end;
R = 8//2
矢 = 2
y = R - 矢
r = 1//2
iniv = BigFloat[2.5, 2.5, -2.6]
res = nls(H, ini=iniv)

    ([2.449489742783178, 2.3979157616563596, -2.597742075127723], true)

界斜は sqrt((x1 - x2)^2 + (sqrt(R^2 - x2^2) - (R - 矢))^2) により求めることができる。
大円の直径が 8 寸，矢が 2 寸，円の直径が 1 寸のとき，5.103103630798287 寸である。

(x, x1, x2) = res[1]
sqrt((x1 - x2)^2 + (sqrt(R^2 - x2^2) - (R - 矢))^2)

    5.103103630798287

function draw(r1, r2, more=false)
    pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = 8//2
    矢 = 2
    r = 1//2
    y = R - 矢
    (x, x1, x2) = res[1]
    界斜 = sqrt((x1 - x2)^2 + (sqrt(R^2 - x2^2) - (R - 矢))^2)
    println("界斜 = $界斜")
    plot()
    circle(0, 0, R, :purple)
    segment(x2, sqrt(R^2 - x2^2), x1, R - 矢, :orange)
    circle2(x, y + r, r, :blue)
    segment(-sqrt(R^2 - y^2), y, sqrt(R^2 - y^2), y)
    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(x2, sqrt(R^2 - x2^2), "(x2,sqrt(R^2-x2^2))", :green, :left, :bottom, deltax=3delta)
        point(x1, y, "(x1,y)", :green, :center, delta=-delta)
        point(R, 0, " R", :purple, :left, :bottom, delta=delta)
        point(x, y + r, "(x,y+r)", :blue, :center, delta=-delta)
        ylims!(0, R)
    end
end;

draw(12/2, 8/2, true)

算額（その1502）

算額（その1502）

六十三 花泉町金沢 大門神社・大門観世音菩薩堂 明治13年(1880)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03064
https://w.atwiki.jp/sangaku/pages/292.html

キーワード：円8個，外円，弦
#Julia, #SymPy, #算額, #和算

外円の中に水平な弦と甲円，乙円を容れ，隙間に丙円 3 個，丁円 2 個を容れる。甲円，乙円の直径がそれぞれ 12 寸，8 寸のとき，丁円の直径はいかほどか。

「今有如図」の図が正しい。これに基づいて得た解は，答，術に一致する。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, r1 - R)
乙円の半径と中心座標を r2, (0, R - r2)
丙円の半径と中心座標を r3, (0, R - r3), (x3, R - 3r3)
丁円の半径と中心座標を r4, (x4, R - 2r3 + r4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, r2::positive, r3::positive, x3::positive, r4::positive, x4::positive
R = r1 + r2
eq1 = x3^2 + (3r3 - r2)^2 - (r2 + r3)^2
eq2 = x4^2 + (r2 - 2r3 + r4)^2 - (r2 + r4)^2
eq3 = x3^2 + (2R - 3r3 - r1)^2 - (r1 + r3)^2
eq4 = x4^2 + (R - 2r3 + r4)^2 - (R - r4)^2
res = solve([eq1, eq2, eq3, eq4], (r3, x3, r4, x4))[1]

    (r2*(r1 + r2)/(2*r1 + r2), 2*sqrt(2)*sqrt(r1)*r2*sqrt(r1 + r2)/(2*r1 + r2), r1*r2/(2*r1 + r2), 2*sqrt(2)*sqrt(r1)*r2*sqrt(r1 + r2)/(2*r1 + r2))

丁円の半径 r4 は，甲円，乙円の半径を r1, r2 とすると r1*r2/(2*r1 + r2) である。
甲円，乙円の直径がそれぞれ 12 寸，8 寸のとき，丁円の直径は 12*8/(2*12 + 8) = 3 寸である。

術は，「丁円径 = 甲円径/(2甲円径/乙円径 + 1) 」で同じ結果を得る。

function draw(r1, r2, more=false)
    pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r3, x3, r4, x4) = (r2*(r1 + r2)/(2*r1 + r2), 2*sqrt(2)*sqrt(r1)*r2*sqrt(r1 + r2)/(2*r1 + r2), r1*r2/(2*r1 + r2), 2*sqrt(2)*sqrt(r1)*r2*sqrt(r1 + r2)/(2*r1 + r2))
    R = r1 + r2
    y = R - 2r3
    x = sqrt(R^2 - y^2)
    @printf("甲円，乙円の直径がそれぞれ %g, %g のとき，丁円の直径は %g である。", 2r1, 2r2, 2r4)
    plot()
    circle(0, 0, R, :purple)
    segment(-x, y, x, y, :orange)
    circle(0, R - r3, r3, :blue)
    circle2(x3, R - 3r3, r3, :blue)
    circle(0, R - r2, r2, :green)
    circle(0, r1 - R, r1, :magenta)
    circle2(x4, R - 2r3 + r4, r4)
    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, r1 - R, "甲円:r1,(0,r1-R)", :magenta, :center, delta=-delta)
        point(0, R - r2, "乙円:r2,(0,R-r2)", :green, :center, delta=-4delta)
        point(0, R - r3, "丙円:r3\n(0,R-r3)", :blue, :center, :bottom, delta=delta)
        point(x3, R - 3r3, "丙円:r3\n(x3,R-3r3)", :blue, :center, delta=-delta)
        point(x4, R - 2r3 + r4, "丁円:r4,(x4,R-2r3+r4)", :red, :center, delta=-delta)
        point(0, R, "R", :purple, :center, :bottom, delta=delta)
    end
end;

draw(12/2, 8/2, true)

算額（その1501）

算額（その1501）

五十八 岩手県花泉町 花泉天満宮 明治3年(1870)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03056
https://w.atwiki.jp/sangaku/pages/167.html

キーワード：円7個，外円，円弧，弦2本，扇形
#Julia, #SymPy, #算額, #和算

外円の中に扇形を容れ，甲円 1 個，乙円 3 個，丙円 2 個を容れる。甲円の直径が 41 寸のとき，丙円の直径はいかほどか。

これも，山村の図は明らかにおかしい。接すべきものが接していないのは手描きのずなのでやむを得ないかもしれないが，扇型の上にある黄色の 3 個の円は同じ乙円ではない。中央にある円と左右の円の大きさは普通に考えれば同じではない。同じだとすれば真ん中の円が外円と扇形に接するならば，左右の円は外円と扇形に同時に接することはできない。
ということで，この算額図も「今有如図」の方が正しい。黄色の真ん中の円と下の赤い 2 この円が乙円，左右の黄色の 2 個が丙円である。そのように仮定して求めた解は，答，術に一致する。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, 0); r1 = R - 2r2
乙円の半径と中心座標を r2, (0, R - r2), (x2, y2)
丙円の半径と中心座標を r3, (x3, y3)
扇型の半径と中心座標を (2R - 2r2), (0, -R)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, r2::positive, r3::positive, x3::positive, y3::positive
R = r1 + 2r2
eq1 = x3^2 + y3^2 - (R - r3)^2
eq2 = x3^2 + (y3 + R)^2 - (r3 + (2R - 2r2))^2
eq3 = x3^2 + (R - r2 - y3)^2 - (r2 + r3)^2
eq4 = r1^2 + (R - r2)^2 - R^2
res = solve([eq1, eq2, eq3, eq4], (r2, r3, x3, y3))[1]

    (r1/3, 12*r1/41, 8*sqrt(10)*r1/41, 151*r1/123)

丙円の半径 r3 は，甲円の半径 r1 の 12/41 倍である。
甲円の直径が 41 寸のとき，丙円の直径は 12 寸である。

算額の問に答えるのはここまでである。
以下では，図を描くためのパラメータを求める。

扇型と外円の交点座標を (x0, y0)
扇型の半径と外円に接する乙円の中心座標を (x2, y2)
として加え，8 元連立方程式を解く。

@syms x0, y0, x2, y2
eq5 = x0^2 + y0^2 - R^2
eq6 = x0^2 + (y0 + R)^2 - (2R - 2r2)^2
eq7 = dist2(0, -R, x0, y0, x2, y2, r3)
eq8 = x2^2 + y2^2 - (R - r2)^2
eq7 = y2/x2*(y0 + R)/x0 + 1;

res2 = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (r2, r3, x3, y3, x0, y0, x2, y2))
res2[4]

    (r1/3, 12*r1/41, 8*sqrt(10)*r1/41, 151*r1/123, 8*r1/5, 7*r1/15, 16*r1/15, -4*r1/5)

甲円の直径が 41 のとき，全パラメータは以下のとおりである。

r2 = 6.83333;  R = 34.1667;  r1 = 20.5;  r3 = 6;  x3 = 12.6491;  y3 = 25.1667;  x2 = 21.8667;  y2 = -16.4

function draw(r1, more=false)
    pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r2, r3, x3, y3, x0, y0, x2, y2) = (r1/3, 12*r1/41, 8*sqrt(10)*r1/41, 151*r1/123, 8*r1/5, 7*r1/15, 16*r1/15, -4*r1/5)
    R = r1 + 2r2
    @printf("r2 = %g;  R = %g;  r1 = %g;  r3 = %g;  x3 = %g;  y3 = %g;  x2 = %g;  y2 = %g\n", r2, R, r1, r3, x3, y3, x2, y2)
    plot()
    circle(0, 0, R, :green)
    circle(0, R - r2, r2)
    circle2(x3, y3, r3, :blue)
    circle(0, 0, r1, :magenta)
    θ = atand(y0 + R, x0)
    circle(0, -R, 2R - 2r2, :orange, beginangle=θ, endangle=180 - θ)
    plot!([-x0, 0, x0], [y0, -R, y0], color=:orange, lw=0.5)
    circle2(x2, y2, r2, :red)
    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, 0, "外円:R", :green, :center, :bottom, delta=delta)
        point(0, 0, "甲円:r1", :magenta, :center, delta=-delta)
        point(0, R - r2, "乙円:r2\n(0,R-r2)", :red, :center, delta=-delta)
        point(x2, y2, "乙円:r2\n(x2,y2)", :red, :center, delta=-delta)
        point(x3, y3, "丙円:r3\n(x3,y3)", :blue, :center, delta=-delta)
        point(0, R, "R", :green, :center, :bottom, delta=delta)
        point(x0, y0, "(x0,y0)  ", :orange, :right, :vcenter)
    end
end;

draw(41/2, true)

算額（その1500）

算額（その1500）

五十八 岩手県花泉町 花泉天満宮 明治3年(1870)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03056
https://w.atwiki.jp/sangaku/pages/167.html

キーワード：円10個
#Julia, #SymPy, #算額, #和算

交差する 2 個の甲円の 2 つの交点を通る乙円を描き隙間に丙円 3 個，丁円 4 個を容れる。丙円の直径が 4 寸のとき，丁円の直径はいかほどか。

山村は丙円を 6 個，丁円を 8 個描いているが，問との不一致をどのように思っていたのか?術の解説にもその点には触れておらず，実際に問を解いたとは思えない。山村のような図は描けないこともないが，連立方程式の本数が多くなるだろうと，手をつけないままにしていた。
これも，「今有如図」が正しく，その解釈に基づいて得られる解は答，術とも一致する。

甲円の半径と中心座標を r1, (r1 - r3, 0), (r3 - r1, 0)
乙円の半径と中心座標を r2, (0, 0)
丙円の半径と中心座標を r3, (0, 0), (r2 + r3)
丁円の半径と中心座標を r4, (x4, y4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive, r4::positive, x4::positive, y4::positive
eq1 = x4^2 + y4^2 - (r2 + r4)^2
eq2 = (x4 - (r1 - r3))^2 + y4^2 - (r1 - r4)^2
eq3 = (r1 - r3 - x4)^2 + y4^2 - (r3 + r4)^2
eq4 = (2r3 + r2) - (2r1 - r3)
eq5 = (r1 - r3)^2 + r2^2 - r1^2
res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, r4, x4, y4))[1]

    (5*r3/2, 2*r3, 3*r3/4, 9*r3/4, sqrt(10)*r3/2)

丁円の半径 r4 は，丙円の半径 r3 の 3/4 倍である。
丙円の直径が 4 寸のとき，丁円の直径は 3 寸である。

function draw(r3, more=false)
    pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, r2, r4, x4, y4) = (5*r3/2, 2*r3, 3*r3/4, 9*r3/4, sqrt(10)*r3/2)
    @printf("r1 = %g;  r2 = %g;  r3 = %g;  r4 = %g;  x4 = %g;  y4 = %g\n", r1, r2, r3, r4, x4, y4)
    plot()
    circle2(r1 - r3, 0, r1)
    circle(0, 0, r3, :blue)
    circle2(2r1 - 2r3, 0, r3, :blue)
    circle(0, 0, r2, :green)
    circle4(x4, y4, r4, :magenta)
    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(r2 + r3, 0, "丙円:r3\n(r2+r3,0)", :blue, :center, delta=-delta)
        point(x4, y4, "丁円:r4\n(x4,y4)", :magenta, :center, delta=-delta)
        point(r2, 0, "r2 ", :green, :right, delta=-delta)
        point(r3, 0, " r3", :blue, :left, delta=-delta)
        point(r1 - r3, 0, "甲円", :red, :center, :bottom, delta=delta)
        point(0, 0, "乙円:r2", :green, :center, :bottom, delta=delta)
        point(0, 0, "丙円:r3", :blue, :center, delta=-delta)
    end
end;

draw(4/2, true)

 

算額（その1499）

算額（その1499）

四十七 岩手県一関市平沢 平沢白山神社 慶応2年(1866)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03056
https://w.atwiki.jp/sangaku/pages/330.html

キーワード：円3個，円弧2個
#Julia, #SymPy, #算額, #和算

直径が同じ 2 個の円弧が交差してできる月形の中に等円 3 個を容れる。等円の直径が 873 寸のとき，円弧の直径はいかほどか。
注：左側の円弧は右側の円弧の中心を通る。

円弧の半径と中心座標を R, (0, 0), (-R, 0)
等円の半径と中心座標を r, (R - r, 0), (x, y)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r::positive, x::positive, y::positive
eq1 = x^2 + y^2 - (R - r)^2
eq2 = (R - r - x)^2 + y^2 - 4r^2
eq3 = (R + x)^2 + y^2 - (R + r)^2
res = solve([eq1, eq2, eq3], (R, x, y))[1]

    (r*(sqrt(57) + 9)/6, r*(15 - sqrt(57))/12, r*sqrt(-1/8 + 3*sqrt(57)/8))

円弧の半径 R は，等円の半径 r の (√57) + 9)/6 倍 である。
等円の直径が 873 寸のとき，円弧の半径は 873*(√57 + 9)/6 = 2408.000910331894 寸である。

function draw(r, more=false)
    pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (R, x, y) = (r*(sqrt(57) + 9)/6, r*(15 - sqrt(57))/12, r*sqrt(-1/8 + 3*sqrt(57)/8))
    @printf("等円の直径が %gのとき，円弧の直径は %g である。\n", 2r, 2R)
    plot()
    circle(0, 0, R, beginangle=240, endangle=480)
    circle(-R, 0, R, beginangle=300, endangle=420)
    circle22(x, y, r, :blue)
    circle(R - r, 0, r, :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(R, 0, " R", :red, :left, :vcenter)
        point(R - r, 0, "等円:r,(R-r,0)", :blue, :center, delta=-delta/2)
        point(x, y, "等円:r,(x,y)", :blue, :center, delta=-delta/2)
    end
end;

draw(873/2, true)