算額（その1649）

算額（その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)

算額（その1648）

算額（その1648）

三十二 一関市舞川 観福寺内地蔵堂 明治34年(1901)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03086
https://w.atwiki.jp/sangaku/pages/158.html

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

水平線上に合同な直角三角形（緑）を縦・横に配し，甲乙の 2 本の斜線（紫）を引いてそれぞれの直角三角形の面積を等分割する。分割された区画に大円 1 個，小円 2 個を容れる。小円の直径が与えられたとき，大円の直径を求める術を述べよ。

山村の図は細部が不適切なので（算額へのリスペクトがまったくない），今回も「今有如図」が示す図に基づいて解を求める。その結果，術と同じ解が得られた。山村の術の解説には誤りがあることもわかった。

直角三角形の直角を挟む二辺のうち，短い方を「鈎」，長い方を「股」とする。
大円の半径と中心座標を r1, (x1, y1)
小円の半径と中心座標を r2, (股/2, r2), (股 + 鈎 - r2, 股/2)
とおき，以下の連立法手式を解く。
鈎（または股）が既知として，順次パラメータを確定する手順を示し，最後に一挙に連立方程式を解く手順を示す。

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

using SymPy

@syms 鈎::positive, 股::positive, l::positive,
      r1::positive, x1::positive, y1::positive,
      r2::positive;

まず最初に，図のように甲斜，乙斜（2 本の斜線）がそれぞれの直角三角形の面積を等分割するということは，甲斜は 2 つの合同な直角三角形の弦（斜辺）の中点 A, B を結ぶ直線が左側の直角三角形の直角の頂点 O を通らなければならない。すなわち，鈎，股（直角を挟む二辺）の間に特定の関係が成り立っている必要がある。

ここでは仮に，鈎が事前に決められている場合について考える。

eq1 = 鈎/股 - (股/2)/(股 + 鈎/2);

# 股
ans_股 = solve(eq1, 股)[1]
ans_股 |> println

    鈎*(1 + sqrt(2))

股は鈎の (1 + √2) 倍でなければならない。

# 鈎
ans_鈎 = solve(eq1, 鈎)[1]
ans_鈎 |> println

    股*(-1 + sqrt(2))

逆に，股が事前に決められている場合には 鈎は股の (√2 - 1)倍でなければならない。

2変数の間の関係式（方程式）は同じで，どちらを未知数として解くかが違うだけである。

次に，小円が内接している 2 つの二等辺三角形は合同であり，小円の半径は弦の長さによって決まる（言うまでもなく，弦は鈎，股によって決まる）。

# r2
弦 = sqrt(鈎^2 + 股^2)
eq2 =  鈎*股/4 - (股 + 弦/2 + 弦/2)*r2/2
ans_r2 = solve(eq2, r2)[1]
ans_r2 |> println

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

また，大円が内接している三角形は，直角二等辺三角形であり，内接円の半径 r1 は，等辺 = 弦/2，斜辺 = √2等辺 により決定される。
2r1 = 等辺 + 等辺 - √2等辺

# r1
等辺 = 弦/2
eq3 = 2r1 - (等辺 + 等辺 - √Sym(2)等辺)
ans_r1 = solve(eq3, r1)[1]
ans_r1 |> println

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

変数は 鈎，股，r1, r2 の 4 変数で，方程式は eq1, eq2, eq3 の 3 本である。
変数のうちの 1 つが既知で，残りの 3 つが未知として連立方程式を立てて一挙に解を求めることができる。

上の例では，鈎が既知で，股，r1, r2 を求める場合は以下のようになる。

res1 = solve([eq1, eq2, eq3], (股, r1, r2))[1]

    (鈎*(1 + sqrt(2)), 鈎*(-1 + sqrt(2))*sqrt(sqrt(2) + 2)/2, 鈎*(1 + sqrt(2))/(2*(1 + sqrt(2) + sqrt(2)*sqrt(sqrt(2) + 2))))

同じ連立方程式で，r2 が 既知で，r1，鈎，股を求める場合は以下のようになる。
算額では，鈎，股を求めることは要求していないが，既知の r2 から未知の r1 を求める術になっている。

res2 = solve([eq1, eq2, eq3], (r1, 鈎, 股))[1]

    (r2*(-2 - sqrt(sqrt(2) + 2) + sqrt(2*sqrt(2) + 4) + 2*sqrt(2)), 2*r2*(-sqrt(2*sqrt(2) + 4) + 1 + 2*sqrt(sqrt(2) + 2)), 2*r2*(1 + sqrt(2) + sqrt(2*sqrt(2) + 4)))

# r1
res2[1] |> factor

大円の半径 r1 は，小円の半径 r2 の (-1 + sqrt(2))*(sqrt(sqrt(2) + 2) + 2) = 1.59379398947637 倍である。
小円の直径が 1 のとき，大円の直径は 1.59379398947637 である。

# r1
res2[1] |> factor |> println

    r2*(-1 + sqrt(2))*(sqrt(sqrt(2) + 2) + 2)

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

    1.59379398947637

「術」は以下のように「大円径は小円径の 1.5937939894763695 倍」としている。上に述べた結果と同じである。

なお，山村は術の解説で，「方斜率 = √3」としている。方斜率は √2 なので，とんでもない間違いである。山村の本は，いろんな間違いがたくさんある，稀代の悪書である。

# 術
小円径 = 1
方斜率 = √2
天 = 2 - √2
大円径 = (√天 + 天*方斜率)*小円径

    1.5937939894763695

図を描くためだけに大円の中心座標 (x1, y1) の解析解を求めるのは，労多くして益少しなので数値解を求める。

eq1 = 鈎/股 - (股/2)/(股 + 鈎/2);
弦 = sqrt(鈎^2 + 股^2)
eq2 =  鈎*股/4 - (股 + 弦/2 + 弦/2)*r2/2
等辺 = 弦/2
eq3 = 2r1 - (等辺 + 等辺 - √Sym(2)等辺)
eq4 = dist(0, 鈎, 股, 0, x1, y1) - r1^2
eq5 = dist(股, 0, 股 + 鈎, 股, x1, y1) - r1^2;

function H(u)
    (r1, 鈎, 股, x1, y1) = u
    return [
        -股/(2*(股 + 鈎/2)) + 鈎/股,  # eq1
        -r2*(股 + sqrt(股^2 + 鈎^2))/2 + 股*鈎/4,  # eq2
        2*r1 - sqrt(股^2 + 鈎^2) + sqrt(2)*sqrt(股^2 + 鈎^2)/2,  # eq3
        -r1^2 + (x1 - 股*(x1*股 - 鈎*(y1 - 鈎))/(股^2 + 鈎^2))^2 + (y1 - 鈎 + 鈎*(x1*股 - 鈎*(y1 - 鈎))/(股^2 + 鈎^2))^2,  # eq4
        -r1^2 + (y1 - 股*(y1*股 + 鈎*(x1 - 股))/(股^2 + 鈎^2))^2 + (x1 - 股 - 鈎*(y1*股 + 鈎*(x1 - 股))/(股^2 + 鈎^2))^2,  # eq5
    ]
end;
r2 = 1 # 0.11769
iniv = BigFloat[1.59, 4.16, 10.05, 9.19, 2.08]
res = nls(H, ini=iniv)

    ([1.5937939894763693, 4.164784400584788, 10.054678984251696, 9.192123892710637, 2.082392200292394], true)

function draw(r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, 鈎, 股, x1, y1) = r2 .* [1.5937939894763693, 4.164784400584788, 10.054678984251696, 9.192123892710637, 2.082392200292394]
    plot([0, 股 + 鈎, 股 + 鈎, 股, 0, 0],
         [0, 0, 股, 0, 鈎, 0], color=:green, lw = 0.5)
    plot!([0, 股 + 鈎/2, 股 + 鈎],
         [0, 股/2, 0], color=:magenta, lw = 0.5)
    circle(x1, y1, r1)
    circle(股/2, r2, r2, :blue)
    circle(股 + 鈎 - r2, 股/2, 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(股/2, 鈎/2, "A", :green, :center, :bottom, delta=delta)
        point(股 + 鈎/2, 股/2, " B", :green, :left, :vcenter)
        point(0, 0, "O", :green, :center, delta=-2delta)
        point(0, 鈎, "鈎 ", :green, :right, :vcenter)
        point(股, 0, "股", :green, :center, delta=-2delta)
        point(股 + 鈎, 0, "股+鈎", :green, :center, delta=-2delta)
        point(股 + 鈎, 股, "(股+鈎,股) ", :green, :right, :vcenter)
        point(x1, y1, "大円:r1\n(x1,y1)", :red, :center, delta=-delta)
        point(股/2, r2, "小円:r2,(股/2,r2)", :blue, :center, delta=-delta, deltax=-10delta)
        point(股 + 鈎 - r2, 股/2, "小円:r2,(股+鈎-r2,股2)", :blue, :right, :bottom, delta=delta, deltax=-8delta)
        plot!(xlims=(-8delta, 股 + 鈎 + 5delta), ylims=(-8delta, 股 + 5delta))
    end
end;

draw(1/2, true)

算額（その1404）改訂版

算額（その1404）　改訂版

算額（その1404）は，依拠した図がでたらめなものであったので，改訂版を書いた。

三十二 一関市舞川 観福寺内地蔵堂 明治34年(1901)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

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


キーワード：円2個，直角三角形4個，等脚台形
#Julia, #SymPy, #算額, #和算

問の原文は以下のとおりである。
「今有如図半円内設等円二個従其親所洩二斜容等円一個其等円径若干問得半円径術如何」

山村の図は正しいものではない。「今有如図」の図に従って解を求めると，術に一致する解が得られた。

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

斜線は，(x0, y0), (x, 0), (0, y01) を通る。

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

@syms R::positive, r::positive, x::positive, x0::positive, y0::positive
eq1 = dist2(x, 0, x0, y0, 0, -r, r)
eq2 = x^2 + r^2 - (R - r)^2
eq3 = (x0 - x)^2 + (y0 - r)^2 - r^2
eq4 = x0^2 + y0^2 - R^2
res = solve([eq1, eq2, eq3, eq4], (R, x, x0, y0))[1]

    (3*r, sqrt(3)*r, 3*sqrt(3)*r/2, 3*r/2)

円弧の半径 R は，等円の半径 r の 3 倍である。
等円の直径が 1 寸のとき，円弧直径は 3 寸である。
x = 0.866025，x0 = 1.29904, y0 = 0.75

function draw(r, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (R, x, x0, y0) = (3*r, sqrt(3)*r, 3*sqrt(3)*r/2, 3*r/2)
    @printf("等円の直径が %g のとき，円弧の直径は %g である。ただし，x = %g，x0 = %g, y0 = %g である。\n", 2r, 2R, x, x0, y0)
    plot()
    circle(0, 0, R, beginangle=0, endangle=180)
    circle2(x, r, r, :blue)
    circle(0, -r, r, :blue)
    (x01, y01) = float.(intersection(x0, y0, x, 0, -x0, y0, -x, 0))
    println((x01, y01))
    segment(x0, y0, 0, y01)
    segment(-x0, y0, 0, y01)
    segment(-R, 0, R, 0, :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(x0, y0, "(x0,y0)", :blue, :right, :vcenter, deltax=-delta)
        point(x, 0, "x", :blue, :center, delta=-delta)
        point(R, 0, "R", :red, :center, delta=-delta)
        point(x, r, "等円:r,(x,r)", :blue, :center, delta=-delta, deltax=-2delta)
        point(0, -r, "等円:r,(0,-r)", :blue, :center, delta=-delta, deltax=-2delta)
        point(0, y01, "y01", :black, :left, :vcenter, deltax=delta)
    end
end;

draw(1/2, true)

 

算額（その1403） 改訂版

算額（その1403）　改訂版

算額（その1403）は，依拠した図がでたらめなものであったので，改訂版を書いた。

三十二 一関市舞川 観福寺内地蔵堂後額 明治34年(1901)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03086
https://w.atwiki.jp/sangaku/pages/158.html

キーワード：円2個，等脚台形，斜線
#Julia, #SymPy, #算額, #和算, #数学

合同な直角三角形を4個組み合わせて等脚台形を作り，大円 1 個，小円 1 個を容れる。小円の直径が与えられたときに大円の直径を得る術を述べよ。

山村の図は正しいものではない。「今有如図」の図に従って解を求めると，術に一致する解が得られた。

直角三角形の直角を挟む 2 辺の短い方，長い方の長さをそれぞれ a, h
大円の半径と中心座標を r1, (2a - r1, r1)
小円の半径と中心座標を r2, (a + r2, h - r2)
とおき，以下の連立方程式を解く。

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

using SymPy

@syms a::positive, h::positive,
      r1::positive, r2::positive
eq1 = dist2(0, 0, a, h, 2a - r1, r1, r1)
eq2 = dist2(3a, 0, 2a, h, a + r2, h - r2, r2)
eq3 = ((a + r2) - (2a - r1))^2 + (h - r2 - r1)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3], (r1, a, h))[5]

    (-5*r2 + sqrt(41)*r2, -3*r2/2 + sqrt(41)*r2/2, r2*(3 + sqrt(41))/2)

大円の半径 r1 は，小円の半径 r2 の (√41 - 5) = 1.4031242374328485 倍である。
小円の直径が 1 寸のとき，大円の直径は 1.4031242374328485 寸である。

function draw(r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, a, h) = (-5*r2 + sqrt(41)*r2, -3*r2/2 + sqrt(41)*r2/2, r2*(3 + sqrt(41))/2)
    plot([0, 3a, 2a, a, 0], [0, 0, h, h, 0], color=:green, lw=0.5)
    plot!([a, a, 2a, 2a], [0, h, 0, h], color=:green, lw=0.5)
    circle(2a - r1, r1, r1)
    circle(a + r2, h - 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(2a - r1, r1, "大円:r1\n(2a-r1,r1)", :red, :center, delta=-delta, deltax=2delta)
        point(a + r2, h - r2, "小円:r2\n(a+r2,h-r2)", :blue, :center, delta=-delta)
        point(a, h, "(a,h)", :green, :center, :bottom, delta=delta)
        point(2a, h, "(2a,h)", :green, :center, :bottom, delta=delta)
        point(3a, 0, "3a", :green, :right, :bottom, delta=delta/2, deltax=-delta)
    end
end;

draw(1/2, true)

算額（その1647）

算額（その1647）

三四 武蔵国埼玉郡下忍村遍照院境内 金毘羅社（神楽堂） 天保11年(1840)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円7個，楕円3個，正方形
#Julia, #Julia, #SymPy, #算額, #和算, #数学

正方形の中に，交差する大楕円 2 個を設け，四隅に丙円 4 個，中央部に小楕円，更にその中に甲円 1 個，乙円 2 個を容れる。乙円とは丙円はそれぞれの楕円の長径端において楕円と 1 点で接する最大の円である（曲率円）。小楕円の短径が 1 寸のとき，大楕円の短径はいかほどか。

大楕円の長半径，短半径，中心座標を a, b, (0, 0)
小楕円の長半径，短半径，中心座標を c, d, (0, 0)
甲円の半径と中心座標を r1, (0, 0); r1 = b
乙円の半径と中心座標を r2, (b - r2, 0)
丙円の半径と中心座標を r4, ((c - r4)/√2, (c - r4)/√2)
大楕円と小楕円の接点座標を A:(x0, y0)
接点から長軸へおろした垂線の脚を B，原点を O としたとき AB = y1, 0B = x1
とおき，以下の連立方程式を解く。

小楕円の中の乙円は曲率円なので，r2 = b^2/a
また，a = b + 2r2 である。

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

using SymPy

@syms a::positive, b::positive, r2::positive,
      c::positive, d::positive, r4::positive,
      x0::positive, y0::positive, x1::positive, y1::positive

eq01 = r2 - b^2/a
eq02 = a - (b + 2r2)
solve([eq01, eq02], (a, r2))[1]

    (2*b, b/2)

a = 2b, r2 = b/2 である。

eq03 = r4 - d^2/c
eq04 = c - (d + 2r4)
solve([eq03, eq04], (c, r4))[1]

    (2*d, d/2)

c = 2d, r4 = d/2 である。

大楕円，小楕円は相似である。

y1 = (x0 + y0)/√Sym(2)
eq0 = y0^2 - (4b^2 - x0^2)/4
eq4 = (x0^2 + y0^2) - ((x0 + y0)^2/2 + x1^2)
eq5 = 3√Sym(2)x0*y0/(4y0 + x0) - 3x1/4
eq6 = x1^2 + 4y1^2 - 4d^2
res = solve([eq0, eq4, eq5, eq6], (d, x0, y0, x1))[2]  # 2 of 2

    (b*sqrt(3*sqrt(41) + 25)/4, b*sqrt(205 - sqrt(41))*(41*sqrt(2) + 13*sqrt(82))/1312, b*sqrt(205 - sqrt(41))*(-3*sqrt(82) + 41*sqrt(2))/1312, b*sqrt(8405 - 41*sqrt(41))/82)

# d
res[1] |> println
res[1].evalf() |> println

    b*sqrt(3*sqrt(41) + 25)/4
    1.66225322815709*b

大楕円の短半径 d は，小楕円の短半径 b の sqrt(3√41+ 25)/4 = 1.66225322815709 倍」である。
小楕円の短径が 1 寸のとき，大楕円の短径は 1.66225322815709 寸である。

術は「大楕円の短径 d は，小楕円の短径 b の sqrt(√369 + 25)/4 倍」である。
sqrt(√369 + 25)/4 = sqrt(3√41 + 25)/4 なので，同じである。

function draw(b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (d, x0, y0, x1) = (b*sqrt(3*sqrt(41) + 25)/4,
        b*sqrt(205 - sqrt(41))*(41*sqrt(2) + 13*sqrt(82))/1312,
        b*sqrt(205 - sqrt(41))*(-3*sqrt(82) + 41*sqrt(2))/1312,
        b*sqrt(8405 - 41*sqrt(41))/82)
    f = d/b
    a = 2b
    c = 2d
    r1 = b
    @printf("factor = %g\n", f)
    @printf("a = %g;  b = %g;  c = %g;  d = %g\n", a, b, c, d)
    r2 = b^2/a  # 曲率円
    r3 = r2*f
    r4 = d^2/c  # 曲率円
    s = sqrt(2(c^2 + d^2))  # 算法助術−94
    plot(s/2 .*[1, 1, -1, -1, 1], s/2 .*[-1, 1, 1, -1, -1], color=:orange, lw=0.5)
    ellipse(0, 0, a, b, color=:red)
    ellipse(0, 0, a, b, color=:pink, φ=90)
    circle(0, 0, b, :green)
    circle2(a - r2, 0, r2, :brown)
    ellipse(0, 0, c, d, color=:blue, lw=0.5, φ=45)
    ellipse(0, 0, c, d, color=:blue, lw=0.5, φ=135)
    ellipse(0, 0, c, d, color=:gray80, lw=0.5, φ=0)
    circle(0, 0, d, :lightgreen)
    circle4((c - r4)/√2, (c - r4)/√2, r4, :deeppink)
    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, "甲円:r1,(0,0)", :green, :center, :bottom, delta=2delta)
        point(a - r2, 0, "乙円:r2\n(a-r2,0)", :brown, :center, delta=-delta/2)
        point((c - r4)/√2, (c - r4)/√2, "丙円:r4\n((c-r4)/√2,\n(c-r4)/√2)", :magenta, :center, delta=-delta/2)
        point(x0, y0, "(x0,y0)", :red, :left, :bottom, delta=delta, deltax=-delta/2)
        point(0, b, "b", :red, :center, :bottom, delta=delta/2)
        point(0, d, "d", :red, :center, :bottom, delta=delta/2)
        point(a, 0, "a", :red, :right, :vcenter, deltax=-delta/2)
        point(c, 0, "c", :red, :right, :vcenter, deltax=-delta/2)
        segment(0, 0, -s/2, -s/2, :gray80)
        point(0, 0, "O", :black, :left, delta=-delta)
        segment(0, 0, -x0, y0)
        point(-x0, y0, "A", :black, :right, :vcenter, deltax=-delta)
        segment(0, 0, -x1/√2, -x1/√2)
        point(-x1/√2, -x1/√2, "B", :black, :right, :vcenter, deltax=-delta)
        segment(-x1/√2, -x1/√2, -x0, y0)
    end
end;

draw(1, true)

算額（その1646）

算額（その1646）

三十二 一関市舞川 観福寺内地蔵堂後額 明治34年(1901)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03086
https://w.atwiki.jp/sangaku/pages/158.html

キーワード：円11個，楕円5個
#Julia, #SymPy, #算額, #和算, #数学

小円を中心として楕円を 5 個配置する。楕円の内外に大円を 10個置く。小円の直径が与えられたときに大円の直径を求める術をのべよ。

図の解釈としては「今有如図」を採った。

楕円の長半径と短半径と中心座標を a, b, (0, r2 + b)
大円の半径と中心座標を r1, (0, r2 + b)
小円の半径と中心座標を r2, (0, 0)
楕円と右隣の楕円の接点座標を (x0, y0)
とおき，以下の連立方程式を解く。

なお，方程式中に tan(54°) = sqrt(2*sqrt(5)/5 + 1) の定数が出てくるが，これをそのまま用いると式が複雑になりすぎて SymPy では解がもとまらない。t = tan(54°) として，t をつかって方程式を解く。

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,
      x0::positive, y0::positive
@syms a, b, r1, r2, x0, y0, t
# t = tand(Sym(54))
a = r1
eq1 = x0^2/a^2 + (y0 - (r2 + b))^2/b^2 - 1
eq2 = -b^2*x0/(a^2*(y0 - (r2 + b))) - t
eq3 = (r2 + b)/2r1 - t;
eq4 = t - y0/x0;

res = solve([eq1, eq2, eq3, eq4], (r1, b, x0, y0))[2]  # 2 of 2

    (r2*(4*sqrt(3) + 7)/(t*(sqrt(3) + 2)), r2*(3 + 2*sqrt(3)), r2*(sqrt(3) + 2)/(2*t), r2*(sqrt(3) + 2)/2)

大円の半径は以下のようになる。

# r1
res[1] |> println

    r2*(4*sqrt(3) + 7)/(t*(sqrt(3) + 2))

t に tan(54°) を代入すると以下のようになり，SymPy では直接的には簡約化できない。

res[1](t => tand(Sym(54)))

SymPy の力を借りながら手作業で簡約化して最終的に以下を得る。

大円の半径は，小円の半径の sqrt(2√5 + 5)*(√5 - 2)*(√3 + 2) = 2.71149362837554 倍である。
たとえば，小円の直径が 1 のとき，大円の直径は 2.71149362837554 である。

# r1
r1 = r2*sqrt(2√5 + 5)*(√5 - 2)*(√3 + 2)
r1 |> println

    2.71149362837554*r2

r1(r2 => 1)

    2.71149362837554

術は以下のようであるが，小円径 = 1 のとき，大円径 = 3.804226065180614 になってしまう。
術がおかしいのか，算額の図形の解釈が間違えているのかわからない。
一つの可能性としては山村の図のように 10 個の楕円の中心が同一円の円周上にあるような図形であるが，ちょっと立式が難しいか？

法 = sqrt(sqrt(0.8) + 1)
小円径 = 1
大円径 = (sqrt(5) + 3)*小円径/法

    3.804226065180614

function draw(r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    t = tand(54)
    (r1, b, x0, y0) = (r2*(4*sqrt(3) + 7)/(t*(sqrt(3) + 2)), r2*(3 + 2*sqrt(3)), r2*(sqrt(3) + 2)/(2*t), r2*(sqrt(3) + 2)/2)
    r1 = r2*sqrt(2√5 + 5)*(√5 - 2)*(√3 + 2)
    a = r1
    @printf("r2 = %g;  r1 = %g;  a = %g;  b = %g;  x0 = %g;  y0 = %g\n", r2, r1, a, b, x0, y0)
    plot()
    circle(0, 0, r2, :green)
    l = (r2 + b)
    l2 = l/cosd(36)
    for i = 0:4
        θ = 90 + 72i
        ellipse(l*cosd(θ), l*sind(θ), a, b, color=:blue, φ=θ+90)
        circle(l*cosd(θ), l*sind(θ), r1, :red)
        circle(l2*cosd(θ + 36), l2*sind(θ + 36), r1, :red)
        i == 4 &&  segment(0, 0, l2*cosd(θ + 36), l2*sind(θ + 36), :black)
    end
    x = l2.*cosd.(54:72:414)
    y = l2.*sind.(54:72:414)
    plot!(x, y, color=:gray80, 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(x0, y0, " (x0,y0)", :blue, :left, :vcenter)
        point(0, r2 + b, "大円:r1\n(0,r2+b)", :red, :center, delta=-delta/2)
    end
end;

draw(1/2, true)

算額（その1645）

算額（その1645）

三十二 一関市舞川 観福寺内地蔵堂後額 明治34年(1901)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

大円の中に斜線を設け，中小の象限を 1 個ずつ描く。それぞれの象限長が与えられたとき，大円の直径を求める術を述べよ。
注：象限の半径は斜線を共有する。

キーワード：円1個，四分円2個，斜線
#Julia, #SymPy, #算額, #和算, #数学

山村の図は「象限」というのが何を指すのか理解できていない。
「象限」とは「四分円」である。中象限（赤），小象限（青）で，∠CAD = ∠DBE = 90°
象限長は四分円の半径で，中象限長 = AC = AD，小象限長 = BD = BE
∠CAD = 90° なので，CF は円の中心を通る直径である。
∠DBE = 90° なので，BD = BE = BF である。
△ACF は直角三角形で CF^2 = AC^2 + AF^2 = AC^2 + (AD + 2BD)
直径^2 = 中象限長^2 + (中象限長 + 2小象限長)^2
直径を 2R, 中象限長を r1，小象限長を r2 とすれば
(2R)^2 = r1^2 + (r1 + 2r2)^2
である。
r1, r2 が既知のとき R は
R = sqrt(r1^2 + (r1 + 2r2)^2)/2
である。

using SymPy
@syms r1::positive, r2::positive, R::positive
eq1 = r1^2 + (r1 + 2r2)^2 - (2R)^2
ans_R = solve(eq1, R)[1]
ans_R |> println

    sqrt(2*r1^2 + 4*r1*r2 + 4*r2^2)/2

たとえば，中象限長が 80，小象限長が 35 のとき，大円の直径は 170 である。

2*ans_R(r1 => 80, r2 => 35) |> println

    170

術は以下のようになっている。

@syms 小象限長, 中象限長, A, 大円径
中象限長 = 80
小象限長 = 35
A = (2小象限長 + 中象限長)^2 + 中象限長^2
大円径 = sqrt(A)

    170.0

 

算額（その1644）

算額（その1644）

三十二 一関市舞川 観福寺内地蔵堂後額 明治34年(1901)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円11個，外円，直角三角形4個
#Julia, #SymPy, #算額, #和算, #数学

全円の中に合同な直角三角形を 4 個，大円 4 個，中円 2 個，小円 4 個を容れる。小円の直径が与えられたとき，中円の直径を得る術を述べよ。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (r1, 0), (0, √3r1)
中円の半径と中心座標を r2, (R - r2, 0)
とおき，以下の方程式を順に解いていく。

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,
      r3::positive, x3::positive, y3::positive,
      l::positive;

# r1
eq1 = dist2(0, R, R/√Sym(3), 0, r1, 0, r1)
ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println

    R*(2 - sqrt(3))

# r2
eq2 = r2/(R - r2 - R/√Sym(3)) - 1//2
eq2 = 2r2 - (R - r2 - R/√Sym(3))
ans_r2 = solve(eq2, r2)[1]
ans_r2 |> println

    R*(3 - sqrt(3))/9

# (x0, y0)
eq3 = (R/√Sym(3) + l*cosd(Sym(30)))^2 + (l*sind(Sym(30)))^2 - R^2
ans_l = solve(eq3, l)[1]
ans_l |> println
x0 = R/√Sym(3) + ans_l*cosd(Sym(30)) |> simplify
y0 = ans_l*sind(Sym(30))
(x0, y0) |> println

    R*(-3 + sqrt(33))/6
    (sqrt(3)*R*(1 + sqrt(33))/12, R*(-3 + sqrt(33))/12)

# r3, x3, y3
eq4 = dist2(0, R, x0, y0, x3, y3, r3)
eq5 = y3/x3*(y0 - R)/x0 + 1
eq5 = y3/x3 + x0/(y0 - R)
eq6 = x3^2 + y3^2 - (R - r3)^2
(ans_r3, ans_x3, ans_y3) = solve([eq4, eq5, eq6], (r3, x3, y3))[1]
# r3
ans_r3 |> println

    R*(-43*sqrt(14814 - 1794*sqrt(33)) - 15*sqrt(54318 - 6578*sqrt(33)) + 6144)/12288

中円の直径は，小円の直径の (R*(3 - sqrt(3))/9) / (R*(-43*sqrt(14814 - 1794*sqrt(33)) - 15*sqrt(54318 - 6578*sqrt(33)) + 6144)/12288) = 1.30332288183929 倍である。

p = ans_r2/ans_r3
p |> println

    4096*(3 - sqrt(3))/(3*(-43*sqrt(14814 - 1794*sqrt(33)) - 15*sqrt(54318 - 6578*sqrt(33)) + 6144))

p.evalf() |> println

    1.30332288183929

山村の術の解説は以下のように，「中円径 = 1.1882963312409864 * 小円径」となっているが，比は小さすぎる。
また，無理数は √3，√5 しか出てこない。

@syms 天, A, 小円径, 中円径
小円径 = 1
天 = √3
A = sqrt((√5天 - 3)*2)
# 中円径 = (天 - 1)*2/((天 - A)*3)*小円径
# 中円径/小円径
(天 - 1)*2/((天 - A)*3) |> println

    1.1882963312409864

function draw(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r1 = R*(2 - sqrt(3))
    plot([0, R/√3, 0, -R/√3, 0], [R, 0, -R, 0, R], color=:blue, lw=0.5)
    circle(0, 0, R, :green)
    circle2(r1, 0, r1)
    circle22(0, √3r1, r1)
    l = R*(-3 + sqrt(33))/6
    (x0, y0) = (R/√3 + l*cosd(30), l*sind(30))
    plot!([0, x0, R/√3], [R, y0, 0], color=:blue, lw=0.5)
    plot!([0, x0, R/√3], -[R, y0, 0], color=:blue, lw=0.5)
    plot!(-[0, x0, R/√3], [R, y0, 0], color=:blue, lw=0.5)
    plot!(-[0, x0, R/√3], -[R, y0, 0], color=:blue, lw=0.5)
    r2 = R*(3 - sqrt(3))/9
    circle2(R - r2, 0, r2, :magenta)
    (r3, x3, y3) = (R*(-43*sqrt(14814 - 1794*sqrt(33)) - 15*sqrt(54318 - 6578*sqrt(33)) + 6144)/12288, R*(128*sqrt(3) + 9*sqrt(4938 - 598*sqrt(33)) + 384*sqrt(11) + 7*sqrt(162954 - 19734*sqrt(33)))/6144, sqrt(33)*R/48 + 5*sqrt(33)*R*sqrt(14814 - 1794*sqrt(33))/12288 + 3*R/16 + 43*R*sqrt(14814 - 1794*sqrt(33))/12288)
    circle4(x3, y3, r3, :orange)
    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\n(r1,0)", :red, :center, delta=-delta/2)
        point(0, √3r1, "大円:r1\n(0,√3r1)", :red, :center, delta=-delta/2)
        point(R - r2, 0, "中円:r2,(R-r2,0)", :black, :center, delta=-delta/2)
        point(x3, y3, "小円:r3,(x3,y3)", :black, :center, delta=-delta/2)
        point(x0, y0, "(x0,y0)", :blue, :right, :vcenter, deltax=-2delta)
        point(R/√3, 0, "R/√3", :blue, :left, :bottom, delta=2.5delta, deltax=-1.5delta)
        point(0, R, "R", :green, :center, :bottom, delta=delta/2)
    end
end;

draw(1, true)

算額（その1642）

算額（その1642）

八十五 室根村折壁字大洞 入沢弥栄神社 明治16年(1883)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03067
https://w.atwiki.jp/sangaku/pages/286.html

キーワード：球，3次元
#Julia, #SymPy, #算額, #和算, #数学

台の上に大球を置き，その周りに小球を環状に置く。大球の直径が 30.6 寸，小球の直径が 435/500 寸のとき，小球は何個あるか。

平面図は以下のようになっている。小球は，大球に比べてかなり小さい。

大球の半径と中心座標を r1, (0, 0, r1)
小球の半径と個数を r2, n
大球と小球の中心までの水平距離を x
とおき，以下の連立方程式を解く。

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

using SymPy
@syms n::positive, x::positive, r1::positive, r2::positive
x = 2*sqrt(r1*r2)
eq = x*sind(180/n) - r2
ans_n = solve(eq, n)[2] |> simplify
ans_n |> println

    pi/asin(sqrt(r2)/(2*sqrt(r1)))

ans_n(r1 => 30.6/2, r2 => 435/500/2).evalf() |> println

    37.2190365074204

大球の直径が 30.6 寸，小球の直径が 435/500 寸のとき，小球は 37 個で，0.2190365074204 個分の隙間が開く。

37 個の球を等間隔に配置すれば，それぞれの小球の間の隙間は殆ど見えないだろうが。

「術」は「sqrt(大球径/小球径)*円周率*2」で，結果は 37.263269312788424 であるが，「答」は 65/500 = 0.13 個と 1 個にも満たないということになっているがこれでは術と合わない。端数を書いてしまったということではないか。

そうだとしても，端数が 0.13 とはどのような「円周率」を使ったのであろうか。計算してみると 3.13035698098977 という半端な「円周率」である。

@syms 大球径, 小球径, 円周率
大球径 = 30.6
小球径 = 435/500
eq = 2sqrt(大球径/小球径)*円周率 - 37.13
solve(eq, 円周率)[1] |> println

    3.13035698098977

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r1 = 30.6/2
    r2 = 435/500/2
    n = 37
    x = 2*sqrt(r1*r2)
    plot()
    circle(0, 0, r1)
    circle(0, 0, x, :green)
    rotate(0, x, r2, :blue, angle=360/n)
    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
end;

draw(true)

算額（その1632）

算額（その1632）

三四 武蔵国埼玉郡下忍村遍照院境内 金毘羅社（神楽堂） 天保11年(1840)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：楕円3個，二等辺三角形，外円
#Julia, #Julia, #SymPy, #算額, #和算, #数学

外円の中に二等辺三角形を設け，その内外に等楕円を 3 個容れる。二等辺三角形の外にある楕円は外円と接しており，その長径は最大である。短径が 1 寸のとき，外円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
楕円の長半径，短半径，中心座標を a, b, (0, b)
とおき，以下の連立方程式を解く。

「二等辺三角形の外にある楕円は外円と接しており，その長径は最大である」とは，外円は楕円の曲率円である。よって，R = a^2/b

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

using SymPy
@syms a::positive, b::positive, R::positive, h::positive, x::positive,
      A::positive, H::positive, P::positive, Q::positive,
      x0::positive, y0::positive
eq1 = R - a^2/b;

二等辺三角形に内接する楕円の長軸，短軸については算法助術の公式97が適用できる。

eq2 = A^4*H*(2H - Q)^2 - A^4*Q*(2H - Q)^2 - A^2*H*P^2*(2H - Q)^2
eq2 = eq2(A => 2x, H => h, P => 2a, Q => 2b) |> simplify |> println

    64*x^2*(b - h)^2*(-a^2*h - 2*b*x^2 + h*x^2)

x ≠ 0, b ≠ h なので，-a^2*h - 2*b*x^2 + h*x^2 = 0 という方程式を解くことになる。

eq2 = -a^2*h - 2*b*x^2 + h*x^2;
res = solve([eq1, eq2], (R, a))[1]

    (-x^2*(2*b - h)/(b*h), x*sqrt(-2*b + h)/sqrt(h))

次いで，x, h, x0, y0 を求める。

eq3 = -b^2*x0/(a^2*(y0 - b)) + h/x
eq4 = h/x - x0/y0
eq5 = x0^2/a^2 + (y0 - b)^2/b^2 - 1
eq6 = R^2 - (x^2 + (R - h)^2);

まず，eq3, eq4 から x0, y0 を求める。

res2 = solve([eq3, eq4], (x0, y0));
# x0
res2[x0] |> println

    a^2*b*h/(x*(a^2 - b^2))

# y0
res2[y0] |> println

    a^2*b/(a^2 - b^2)

次いで，eq5, eq6 から x, h を求める

eq5 = eq5(x0 => res2[x0], y0 => res2[y0])(R => res[1], a => res[2]) |> simplify |> numerator
eq6 = eq6(x0 => res2[x0], y0 => res2[y0])(R => res[1], a => res[2])
solve([eq5, eq6], (h, x))[1]  # 1 of 2

    (3*b, 3*b)

h = x = 3b である。

前に求めた R = -x^2*(2*b - h)/(b*h) の x, h に 3b を代入すると R = 3b である。

R = -x^2*(2*b - h)/(b*h)
R(x => 3b, h => 3b) |> println

    3*b

外円の直径は，楕円の短径の 3 倍である。
楕円の短径が 1 寸のとき，外円の直径は 3 寸である。

R = x = h = 3b ということは，「二等辺三角形」は直角二等辺三角形で，その底辺は外円の直径である。

最後の最後に，a = x*sqrt(-2*b + h)/sqrt(h) の，x, h に3b を代入すると a = √3b である。

a = x*sqrt(-2*b + h)/sqrt(h)
a(x => 3b, h => 3b) |> println

    sqrt(3)*b

二等辺三角形の外の 2 個の楕円は，上の二等辺三角形の中の楕円と合同なので，楕円を内包する下の二等辺三角形も直角二等辺三角形であり，楕円の長軸は上の二等辺三角形の底辺と 45° の角をなす。

function draw(b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = 3b
    a = √3b
    h = R
    x = R
    println("R = $R, a = $a, h = $h, x = $x")
    plot([x, 0, -x, x], [0, R, 0, 0], color=:green, lw=0.5)
    circle(0, 0, R, :blue)
    ellipse(0, b, a, b, color=:red)
    ellipse(2b/√2, -2b/√2, a, b, φ=45, color=:red)
    ellipse(-2b/√2, -2b/√2, a, b, φ=135, color=: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)
        (x0, y0) = (0.75, 0.75)
        plot!([0, √2R, 0, 0], [0, 0, -√2R, 0], color=:gray70, lw=0.5, style=:dash)
        plot!([0, -√2R, 0, 0], [0, 0, -√2R, 0], color=:gray70, lw=0.5, style=:dash)
        point(0, b, "楕円:a,b,(0,b)", :red, :center, delta=-delta/2)
        point((R - b)/√2, -(R - b)/√2)
        point(-(R - b)/√2, -(R - b)/√2)
        point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
        point(√2R, 0, " √2R", :blue, :right, :bottom, delta=delta/2)
        point(0, -√2R, " √2R", :blue, :left, :bottom, delta=delta/2)
    end
end;

draw(1/2, true)

算額（その1641）

算額（その1641）

九 岩手県奥州市（旧水沢市佐倉河） 胆沢城 八幡宮 弘化2年(1845)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03076
https://w.atwiki.jp/sangaku/pages/137.html

キーワード：円6個，外円，楕円1個，弦1本
#Julia, #SymPy, #算額, #和算, #数学

全円の中に，弦 1 本，楕円 1 個，大円 4 個，小円 1 個を容れる。小円の直径が与えられたとき，全円の直径を求める術を述べよ。

山村の図は論外である。「今有如図」によって解を求めたが術と異なるものになる。「今有如図」の図も違うのかもしれない。

一般性を失うことなく，楕円の長径，短径を与えて小円と全円の直径を求め，小円と全円の関係式を求めることにしてよい。

図を時計回りに 90° 回転させて解く。

楕円の長半径，短半径，中心座標を a, b, (0, 0)
全円の半径と中心座標を R, (a - R, 0)
大円の半径と中心座標を r1, (x1, y1), (-r1, 0)
小円の半径と中心座標を r2, (-a - r2, 0)
とおき，以下の連立方程式を解く。

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, x1::negative, y1::positive,
      r2::positive, R::positive,
      x0::negarive, y0::positive
# eq0 = 4(a^2 - b^2)*(b^2 - r1^2)/b^2 - (2r1)^2  # 算法助術の公式84
r1 = b*sqrt(a^2 - b^2)/a
x1 = r1 - a
eq1 = r2 - (R - a)
eq2 = x0^2/a^2 + y0^2/b^2 - 1
eq3 = (x0 - x1)^2 + (y0 - y1)^2 - r1^2
eq4 = -b^2*x0/(a^2*y0) + (x0 - x1)/(y0 - y1)
eq5 = (x1 - (a - R))^2 + y1^2 - (R - r1)^2;

SymPy では能力的に解けないので，数値解を求める。

function driver(a, b)
    function H(u)
        (R, r2, y1, x0, y0) = u
        return [
            -R + a + r2,  # eq1
            -1 + y0^2/b^2 + x0^2/a^2,  # eq2
            (y0 - y1)^2 + (a + x0 - b*sqrt(a^2 - b^2)/a)^2 - b^2*(a^2 - b^2)/a^2,  # eq3
            (a + x0 - b*sqrt(a^2 - b^2)/a)/(y0 - y1) - b^2*x0/(a^2*y0),  # eq4
            y1^2 - (R - b*sqrt(a^2 - b^2)/a)^2 + (R - 2*a + b*sqrt(a^2 - b^2)/a)^2,  # eq5
        ]
    end;
    iniv = BigFloat[6, 0.96, 3.46146, -2.65805,1.69398]
    return nls(H, ini=iniv)
end;

res = driver(5, 2)
@printf("R = %g, r2 = %g, R/r2 = %g\n", res[1][1], res[1][2], res[1][1]/res[1][2]) 

    R = 5.9428, r2 = 0.942805, R/r2 = 6.30333

res = driver(7, 3)
@printf("R = %g, r2 = %g, R/r2 = %g\n", res[1][1], res[1][2], res[1][1]/res[1][2]) 

    R = 8.57247, r2 = 1.57247, R/r2 = 5.45159

術には「小円の直径を 5 倍すれば全円の直径が得られる」とあるが，
長半径，短半径が 5, 2 のとき，R/r2 ≒ 6.30333
長半径，短半径が 7, 3 のとき，R/r2 ≒ 5.45159
のように，R/r2 は一定ではない。

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    res = driver(a, b)
    (R, r2, y1, x0, y0) = res[1]
    r1 = b*sqrt(a^2 - b^2)/a
    x1 = r1 - a
    println("x1= $x1")
    r2 = R - a
    println("R = $R, r2 = $r2, R/r2 = $(R/r2)")
    @printf("a = %g;  b = %g;  r2 = %g;  r1 = %g;  x1 = %g;  y1 = %g;  R = %g\n", a, b, r2, r1, x1, y1, R)
    p = plot()
    circle(a - R, 0, R, :green)
    ellipse(0, 0, a, b, color=:blue)
    circle2(r1, 0, r1)
    circle22(x1, y1, r1)
    circle(-a - r2, 0, r2, :magenta)
    y = sqrt(R^2 - (R - 2a)^2)
    segment(-a, y, -a, -y, :black)
    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 - R, 0, "全円:R,(a-R,0)", :green, :right, :bottom, delta=delta, deltax=5delta)
        point(x1, y1, "大円:r1,(x1,y1)", :red, :center, delta=-delta)
        point(-r1, 0, "大円:r1,(-r1,0)", :red, :center, delta=-delta)
        point(-a - r2, 0, "小円:r2,(-a-r2,0)", :black, :center, delta=-delta, deltax=4delta)
        point(a, 0, " a", :blue, :left, :vcenter)
        point(0, b, "b", :blue, :center, :bottom, delta=delta)
    else
        point(-a, y1, @sprintf(" R/r2 = %g", R/r2), :black, :left, :vcenter, mark=false)
    end
    return p
end;

draw(5, 2, true)

算額（その1640）

算額（その1640）

九 岩手県奥州市（旧水沢市佐倉河） 胆沢城 八幡宮 弘化2年(1845)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03076
https://w.atwiki.jp/sangaku/pages/137.html

一関博物館 和算に挑戦 令和5年度出題問題(3)［上級問題］（高校生・一般向き）
https://www.city.ichinoseki.iwate.jp/museum/wasan/r5/hard.html

キーワード：円2個，外円，楕円3個
#Julia, #SymPy, #算額, #和算, #数学

大円の中に楕円 3 個と小円を容れる。楕円の面積が最大になるとき，大円の直径を小円の直径で表す術を述べよ。

この算額は「算額（その678）」と求めるものが違うだけで，本質的に同じ問題である。

変数値が正の値を取るように，上下を反転させて解く。
大円の半径と中心座標を R, (0, 0)
小円の半径と中心座標を r, (0, 0)
楕円の超半径，短半径，中心座標を a, b, (0, y)
楕円と円の接点座標を (x1, y1)
楕円同士の接点座標を (x2, y2)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms x1::positive, y1::positive, x2::positive, y2::positive, y::positive, a::positive, b::positive, R::positive
y2 = x2/sqrt(Sym(3))
eq1 = x1^2/a^2 + (y1 - y)^2/b^2 - 1
eq2 = b^2*x1/(a^2*(y - y1)) + x1/y1
eq3 = x1^2 + y1^2 - R^2
eq4 = x2^2/a^2 + (y2 - y)^2/b^2 - 1
eq5 = b^2*x2/(a^2*(y - y2))- 1/sqrt(Sym(3));

eq1, eq2, eq3 を解き，x1, y1, y を求める。

res1 = solve([eq1, eq2, eq3], (x1, y1, y))[1]  # 1 of 2

    (sqrt((-R^2*b^2 + a^4)/(a - b))/sqrt(a + b), a*sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/(a^2 - b^2), sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/a)

eq4, eq5 を解き，x2, y を求める。

res2 = solve([eq4, eq5], (x2, y))[1]

    (a^2/sqrt(a^2 + 3*b^2), sqrt(3*a^2 + 9*b^2)/3)

res1 と res2 の y は等しいので，次の法定式を立てる。

eq = res1[3] - res2[2]
eq |> println

    -sqrt(3*a^2 + 9*b^2)/3 + sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/a

b を求める。

ans_b = solve(eq, b)[1]
ans_b |> println

    a*sqrt(9*R^2 - 12*a^2)/(3*R)

楕円の面積 S は πab である。

S = PI * a * ans_b;

面積が最大になるときの長半径を求めるには，S(a) の導関数 S(a)'を求め，S(a)'= 0 を解く。

max_at_a = solve(diff(S, a), a)[1]
max_at_a |> println

    sqrt(2)*R/2

楕円の面積が最大になるときの長半径は √2R/2，

ans_b = ans_b(a => max_at_a)
ans_b |> println

    sqrt(6)*R/6

同じく短半径は √6R/6

ans_y = res2[2](b =>ans_b)(a => max_at_a)
ans_y |> println

    sqrt(3)*R/3

そのときの，中心座標の y 座標値は √3R/3 である。

小円の半径は r = ans_y - ans_b = (√3/3 - √6/6)*R で，
R = r/(√3/3 - √6/6) = r*(√6 + 2√3) すなわち，大円の直径は小円の直径の (√6 + 2√3) = 5.913591357920932 倍である。

@syms r, d
apart(r / (√Sym(3)/3 - √Sym(6)/6), d) |> simplify |> factor |> println

    r*(sqrt(6) + 2*sqrt(3))

山村も「今有如図」も，術は √(√285 + 28) となっている。山村も「今有如図」も誤記（誤読）には気づいていない（問題を実際に解いていないのであろう）。
一関博物館のページでは，「術では √(√288 + 18) となっている」としている。二重根号を外せば √(√288 + 18) = √6 + 2√3 である。

√(√Sym(288) + 18) |> sympy.sqrtdenest |> println

    sqrt(6) + 2*sqrt(3)

function draw(r, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = r*(√6 + 2√3)
    a = √2R/2
    b = √6R/6
    y = √3R/3
    x1 = sqrt((-R^2*b^2 + a^4)/(a - b))/sqrt(a + b)
    y1 = a*sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/(a^2 - b^2)
    x2 = a^2/sqrt(a^2 + 3*b^2)
    y2 = y2 = √3x2/3
    println((x1, y1))
    println((x2, y2))
    plot()
    circle(0, 0, R)
    circle(0, 0, r, :green)
    ellipse(0, y, a, b, color=:blue)
    ellipse(y*cosd(30), -y*sind(30), a, b, color=:blue, φ=-120)
    ellipse(-y*cosd(30), -y*sind(30), a, b, color=:blue, φ=120)
    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, y1, "(x1,y1)", :green, :left, :bottom, delta=delta/2)
        point(x2, y2, "(x2,y2)", :green, :left, delta=-delta/2)
        point(0, R, "R", :red, :center, :bottom, delta=delta/2)
        point(0, r, "r", :green, :center, :bottom, delta=delta/2)
        point(0, y, " y=r+b", :green, :center, :bottom, delta=delta/2)
    end
end;

draw(1/2, true)

算額（その1639）

算額（その1639）

三 岩手県花巻市太田 音羽山清水観世音堂 明治25年(1892)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03076
https://w.atwiki.jp/sangaku/pages/137.html

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

楕円の中に，合同な正三角形を 3 個容れる。楕円の短径が与えられたときに，長径を求める術を述べよ。

原文は「𫝆󠄃有如圖画側圓容等三角󠄄三個側圓短徑若干問得側圓長徑術如何」である。

山村はなぜか知らないが，外円の中に合同な正三角形を 3 個と，楕円を 2 個容れている。
「今有如図」は問の通りの図を描いている。

問の解釈は今有如図のほうが正しい。計算すると術の通りの答えになる。

楕円の長半径，短半径，中心座標を a, b, (0,y0)
正三角形の一辺の長さを x, 楕円上の正三角形の頂点座標を (x/2, √3x/2), (x, 0), (0 -√3x/2)
とおき，以下の連立方程式を解く。

注：楕円の長径は x 軸上にあるのではない。また，結果的に，楕円は長径≒短径で，ほぼ円に近い。

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

using SymPy
@syms x::positive, a::positive, b::positive, y0::positive
eq1 = (x/2)^2/a^2 + (√Sym(3)x/2 - y0)^2/b^2 - 1
eq2 = x^2/a^2 + y0^2/b^2 - 1
eq3 = y0 + √Sym(3)x/2 - b
res = solve([eq1, eq2, eq3], (a, x, y0))[1]

    (sqrt(42)*b/6, 28*sqrt(3)*b/45, b/15)

楕円の長径 a は 短径 b の √42/6 = sqrt(7/6) 倍である。

術は「置七個以六個除之開平方【之倍】乘短徑得長徑」である。「【之倍】」は誰が入れたものであろうか?
今有如図は「欠字については「現存 岩手の算額」による」としているので，山村か？
山村は，図の解釈を誤っているため，解に合わすためにいつものように変なことをしたと思われる。解説で，術の「sqrt(7÷6)×短径=長径」を，「訂正 sqrt(7÷6)×2×短径=2.16×短径=長径」 と書いている。試しに山村の描いた位置に楕円を描いてみると長径は短径の3倍位になる。

function draw(b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (x, a, y0) = (28*sqrt(3)*b/45, sqrt(42)*b/6, b/15)
    println("b = $b, x = $x, a = $a, y0 = $y0")
    plot([0, x, x/2, 0], [0, 0, √3x/2, 0], color=:blue)
    plot!([0, -x, -x/2, 0], [0, 0, √3x/2, 0], color=:blue)
    plot!([x/2, 0, -x/2, x/2], [0, -√3x/2, 0, 0], color=:blue)
    ellipse(0, y0, a, b, color=: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(a, y0)
        point(-a, y0)
        point(x/2, √3x/2)
        point(x, 0, " x", :blue, :left, :vcenter)
        point(0, y0, "y0", :red, :center, :bottom, delta=delta/2)
        point(0, y0 - b, "y0-b", :red, :center, delta=-delta/2)
        point(0, y0 + b, "y0+b", :red, :center, :bottom, delta=delta/2)
        point(a, y0, "(a,y0)", :red, :left, :bottom, delta=delta/2)
        dimension_line(-a, y0, a, y0, "長径", :black, :center, :bottom, delta=3delta, length=4delta)
        plot!(xlims=(-a - 3delta, a + 8delta), ylims=(y0 - b - 5delta, y0 + b + delta))
    end
end;

draw(1, true)

算額（その1638）

算額（その1638）

長野県埴科郡坂城町 諏訪社 文化2年(1805)
中村信弥「改訂増補長野県の算額」県内の算額1（78）
http://www.wasan.jp/zoho/zoho.html

キーワード：円1個，長方形，対角線
#Julia, #SymPy, #算額, #和算, #数学

長方形の中に円を容れる。長方形の長辺，短辺を「長」，「平」，円によって切り取られる対角線の部分を「帯弦」と呼ぶ。長が 185 寸，平が 80 寸のとき，帯弦を求める術を述べよ。

対角線と円の交点座標を (x1, y1), (x2, y2) とおき，以下の連立方程式を解く。

帯弦は sqrt((x1 - x2)^2 + (y1 - y2)^2) である。

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

using SymPy
@syms 長::positive, 平::positive, 帯弦::positive,
      x1::positive, y1::positive, x2::positive, y2::positive
eq1 = y1/(長 - x1) - 平/長
eq2 = y2/(長 - x2) - 平/長
eq3 = (x1 - (長 - 平/2))^2 + (y1 - 平/2)^2 - (平/2)^2
eq4 = (x2 - (長 - 平/2))^2 + (y2 - 平/2)^2 - (平/2)^2
eq5 = sqrt((x1 - x2)^2 + (y1 - y2)^2) - 帯弦
res = solve([eq1, eq2, eq3, eq4, eq5], (帯弦, x1, y1, x2, y2))[2]  # 2 of 2

    (sqrt(2)*平^(3/2)*sqrt(長)/sqrt(平^2 + 長^2), (sqrt(2)*平^(3/2)*長^(3/2) + 平^2*長 - 平*長^2 + 2*長^3)/(2*(平^2 + 長^2)), (-sqrt(2)*平^(5/2)*sqrt(長) + 平^2*(平 + 長))/(2*(平^2 + 長^2)), (-sqrt(2)*平^(3/2)*長^(3/2) + 平^2*長 - 平*長^2 + 2*長^3)/(2*(平^2 + 長^2)), (sqrt(2)*平^(5/2)*sqrt(長) + 平^2*(平 + 長))/(2*(平^2 + 長^2)))

帯弦は sqrt(2)*平^(3/2)*sqrt(長)/sqrt(平^2 + 長^2) である。

res[1] |> println

    sqrt(2)*平^(3/2)*sqrt(長)/sqrt(平^2 + 長^2)

長が 185 寸，平が 80 寸のとき，帯弦は 68.2871764062511 寸である。

res[1](長 => 185, 平 => 80).evalf() |> println

    68.2871764062511

術は 平*sqrt(2長*平/(長^2 + 平^2)) といっているので上で得た数式に一致する。

(平*sqrt(2長*平/(長^2 + 平^2)))(長 => 185, 平 => 80).evalf() |> println

    68.2871764062511

しかし，答では 68 + 179/625 = 68.2864 と，近似値としては精度の低い値を示している。
書き間違いではないようだ。

function draw(長, 平, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (帯弦, x1, y1, x2, y2) = (sqrt(2)*平^(3/2)*sqrt(長)/sqrt(平^2 + 長^2), (sqrt(2)*平^(3/2)*長^(3/2) + 平^2*長 - 平*長^2 + 2*長^3)/(2*(平^2 + 長^2)), (-sqrt(2)*平^(5/2)*sqrt(長) + 平^2*(平 + 長))/(2*(平^2 + 長^2)), (-sqrt(2)*平^(3/2)*長^(3/2) + 平^2*長 - 平*長^2 + 2*長^3)/(2*(平^2 + 長^2)), (sqrt(2)*平^(5/2)*sqrt(長) + 平^2*(平 + 長))/(2*(平^2 + 長^2)))
    println((帯弦, x1, y1, x2, y2))
    plot([0, 長, 長, 0, 0], [0, 0, 平, 平, 0], color=:green, lw=0.5)
    circle(長 - 平/2, 平/2, 平/2)
    segment(0, 平, 長, 0, :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)
        dimension_line(x1, y1, x2, y2, "帯弦", :black, :left, :bottom, delta=delta, length=7, dx=1, dy=2)
        point(x1, y1, "(x1,y1)", :blue, :right, :vcenter, deltax=-3delta)
        point(x2, y2, "(x2,y2)", :blue, :right, :vcenter, deltax=-3delta)
        point(長, 平, "(長,平)", :green, :right, :bottom, delta=delta)
    end
end;

draw(185, 80, true)

算額（その1637）

算額（その1637）

宮城県石巻市尾崎宮下 久須師神社 明治20年(1887)
徳竹亜紀子，谷垣美保: 2021年度の算額調査，仙台高等専門学校名取キャンパス 研究紀要，第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf
キーワード：楕円，正三角形，正方形
#Julia #SymPy #算額 #和算 #数学

正方形の中に楕円と正三角形を容れる。楕円の長径と短径が与えられたとき，正三角形の一辺の長さを得る術を述べよ。

図を反時計回りに 45° 回転させて考える。

楕円が内接する正方形の一辺の長さ a は「算法助術の公式94」で求めることができる。
正三角形の一辺の長さを b，楕円上の正三角形の頂点座標を (x0, y0) とおき，以下の方程式を解く。

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

using SymPy
@syms a, b, p, q, θ, x0, y0
a = √Sym(2)sqrt(p^2 + q^2)
θ = Sym(45 + 60)
x0 = b*cosd(θ)
y0 = b*sind(θ)
eq1 = x0^2/p^2 + (y0 - a/sqrt(Sym(2)))^2/q^2 - 1
ans_b = solve(eq1, b)[1]  # 1 of 2
ans_b |> println

    p^2*(-2*sqrt(2)*3^(1/4)*q + (sqrt(2) + sqrt(6))*sqrt(p^2 + q^2))/(sqrt(3)*p^2 + 2*p^2 - sqrt(3)*q^2 + 2*q^2)

普通は分母の有理化を行うが，分子の有理化を行ってみる。
分子・分母に num = (2√2*3^(1/4)*q + (√2 + √6)*sqrt(p^2 + q^2)) を掛ける。

num = (2*sqrt(Sym(2))*Sym(3)^(1//4)*q + (sqrt(Sym(2)) + sqrt(Sym(6)))*sqrt(p^2 + q^2));
num2 = numerator(ans_b)*num|> simplify;
den2 = denominator(ans_b)*num |> simplify;
num2/den2 |> simplify |> println
num2/den2 |> simplify |> display

    4*p^2/(2*sqrt(2)*3^(1/4)*q + (sqrt(2) + sqrt(6))*sqrt(p^2 + q^2))

徳竹らが提示したのと同じ式が得られた。

2√Sym(2)p^2/((√Sym(3) + 1)sqrt(p^2 + q^2) + 2fourthroot(Sym(3))*q)

function draw(p, q, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    a = sqrt(2)*sqrt(p^2 + q^2)
    b = 4*p^2/(2*sqrt(2)*3^(1/4)*q + (sqrt(2) + sqrt(6))*sqrt(p^2 + q^2))
    println("a = $a, b = $b")
    plot(a/√2 .*[1, 0, -1, 0, 1], a/√2 .+ a/√2 .*[0, 1, 0, -1, 0], color=:green, lw=0.5)
    ellipse(0, a/√2, p, q, color=:red)
    plot!([0, b*cosd(45), b*cosd(105), 0], [0, b*sind(45), b*sind(105), 0], 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, a/√2, "楕円:p,q,(0,a/√2)", :red, :center, delta=-delta/2)
        dimension_line(0, 0, b/√2, b/√2, "b", :black, :left, delta=-delta/2, deltax=delta, dx=a/70, dy=-a/70, length=0.4)
        dimension_line(0, 0, a/√2, a/√2, "a", :black, :left, delta=-delta/2, deltax=delta, dx=a/18, dy=-a/18, length=0.4)
        point(b*cosd(105), b*sind(105), "(x0,y0)", :blue, :center, :bottom, delta=delta)
    end
end;

draw(6, 3, true)