算額（その1167）

算額（その1167）

九八 鴻巣市三ツ木山王 三木神社 明治28年(1895)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

長方形内に交差する 2 個の楕円を容れる。楕円の長径が 4 寸，短径が 2.4 寸，長方形の短辺が 3 寸のとき，甲斜の長さはいかほどか。
注：甲斜とは，原点対称の 2 個の楕円の交点間の距離の長い方である。

図に示す x1, y1, x2, y2, x01, y01, x02, y02, x03, y03 を未知数として以下の連立方程式を立て，それを解く。

甲斜の長さは 2sqrt(x03^2 + y03^2) で求める。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, side_a::positive, side_b::positive, K::positive,
     x1::positive, y1::positive, x2::positive, y2::negative,
     x01::positive, y01::positive, 
     x02::positive, y02::positive,
     x03::positive, y03::positive
side_b = sqrt((x2 - x1)^2 + (y1 - y2)^2)
side_a = sqrt((x2 + x1)^2 + (y2 + y1)^2)
eq1 = x01^2/a^2 + y01^2/b^2 - 1
eq2 = -b^2*x01/(a^2*y01) - (y2 - y1)/(x2 - x1)
eq7 = (y2 - y01)/(x2 - x01)  - (y2 - y1)/(x2 - x1)
eq3 = x02^2/a^2 + y02^2/b^2 - 1
eq4 = -b^2*x02/(a^2*y02) - (y2 + y1)/(x2 + x1)
eq8 = (y2 - y02)/(x2 - x02) - (y2 + y1)/(x2 + x1)
eq5 = (4a^2 + 4b^2) - (side_a^2 + side_b^2)  # formula-89
eq6 = (y2 - y1)/(x2 - x1) * (y2 + y1)/(x2 + x1) + 1
eq9 = side_b - K
eq10 = y03/x03 - (y2 + y1)/(x2 + x1)
eq11 = x03^2/a^2 + y03^2/b^2 - 1;

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)
   (x1, y1, x2, y2, x01, y01, x02, y02, x03, y03) = u
   return [
       -1 + y01^2/b^2 + x01^2/a^2,  # eq1
       -(-y1 + y2)/(-x1 + x2) - b^2*x01/(a^2*y01),  # eq2
       -1 + y02^2/b^2 + x02^2/a^2,  # eq3
       -(y1 + y2)/(x1 + x2) - b^2*x02/(a^2*y02),  # eq4
       4*a^2 + 4*b^2 - (-x1 + x2)^2 - (x1 + x2)^2 - (y1 - y2)^2 - (y1 + y2)^2,  # eq5
       1 + (-y1 + y2)*(y1 + y2)/((-x1 + x2)*(x1 + x2)),  # eq6
       -(-y1 + y2)/(-x1 + x2) + (-y01 + y2)/(-x01 + x2),  # eq7
       -K + sqrt((-x1 + x2)^2 + (y1 - y2)^2),  # eq9
       -(y1 + y2)/(x1 + x2) + y03/x03,  # eq10
       -1 + y03^2/b^2 + x03^2/a^2,  # eq11
   ]
end;

(a, b, K) = (4/2, 2.4/2, 3)
iniv = BigFloat[0.614, 2.254, 2.315, -0.228, 1.789, 0.539, 1.505, -0.788, 1.315, 0.91]
res = nls(H, ini=iniv)

   ([0.6329571688388319, 2.244853051407938, 2.3204571688388325, -0.2355388027151154, 1.8516704311458714, 0.45351293387424213, 1.5000000000000002, -0.793725393319377, 1.322875655532295, 0.9000000000000001], true)

(x1, y1, x2, y2, x01, y01, x02, y02, x03, y03) = res[1]
甲斜 = 2sqrt(x03^2 + y03^2)

   3.1999999999999993

甲斜は 3.2 寸である。

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

  長辺 = 3.57211;  短辺 = 3;  x1 = 0.632957;  y1 = 2.24485;  x2 = 2.32046;  y2 = -0.235539;  x01 = 1.85167;  y01 = 0.453513;  x02 = 1.5;  y02 = -0.793725;  x03 = 1.32288;  y03 = 0.9

function draw(p, q, a, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (4/2, 2.4/2)
   (x1, y1) = (0.614, 2.254)
   (x2, y2) = (2.315, -0.228)
   (x01, y01) = (1.789, 0.539)
   (x02, y02) = (1.505, -0.788)
   (x1, y1, x2, y2, x01, y01, x02, y02, x03, y03) = res[1]
   side_b = sqrt((x2 - x1)^2 + (y1 - y2)^2)
   side_a = sqrt((x2 + x1)^2 + (y2 + y1)^2)
   @printf("甲斜 = %g;  長辺 = %g;  短辺 = %g\n", 2sqrt(x03^2 + y03^2), side_a, side_b)
   @printf("長辺 = %g;  短辺 = %g;  x1 = %g;  y1 = %g;  x2 = %g;  y2 = %g;  x01 = %g;  y01 = %g;  x02 = %g;  y02 = %g;  x03 = %g;  y03 = %g\n", side_a, side_b, x1, y1, x2, y2, x01, y01, x02, y02, x03, y03)
   plot([x1, -x2, -x1, x2, x1], [y1, -y2, -y1, y2, y1], color=:blue, lw=0.5)
   ellipse(0, 0, a, b, color=:red)
   ellipse(0, 0, a, b, φ= 2atand(y03/x03), color=:red)
   segment(x03, y03, -x03, -y03, :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(x1, y1, " (x1,y1)", :blue, :left, :vcenter)
       point(x2, y2, " (x2,y2)", :blue, :left, :vcenter)
       point(-x2, -y2, " (-x2,-y2)", :blue, :left, :vcenter)
       point(-x1, -y1, " (-x1,-y1)", :blue, :left, :vcenter)
       point(x01, y01, " (x01,y01)", :red, :left, :vcenter)
       point(x02, y02, " (x02,y02)", :red, :left, :vcenter)
       point(x03, y03, " (x03,y03)", :green, :left, :vcenter)
       point(-x03, -y03, " (-x03,-y03)", :green, :left, :vcenter)
       xlims!(-x2 - 5delta, x2+10delta)
   end
end;

算額（その1166）

算額（その1166）

六四 加須市不動岡 総願寺 慶応2年(1866)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：楕円4個，長方形，対角線

長方形の中に対角線を引き，分割された区画に，合同な楕円を 4 個容れる。楕円の短径が 1 寸のとき，長径はいかほどか。

注：対角線を引いてできる三角形のうち，長方形の短辺を一辺とする三角形は正三角形である。

二等辺三角形に内接する楕円の長径，短径については算法助術の公式97による（「和算の心（005）」参照）。

長方形の中心を原点とし，長方形の長辺，短辺を 2a, 2b とする。
正三角形の底辺と高さは，2b, √3b
二等辺正三角形の底辺と高さは，2a, b
以上に基づき，
(1) 2 つの楕円の長径を求め，それが等しいこと
(2) a, b の関係式
の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, side_a::positive, height::positive, q::positive
p = side_a*sqrt((height - 2*q)/height)/2
p1 = p(side_a => 2b, height => √Sym(3)b);
p1 |> println

   3^(3/4)*sqrt(b)*sqrt(sqrt(3)*b - 2*q)/3

p2 = p(side_a => 2a, height => b)
p2 |> println

   a*sqrt(b - 2*q)/sqrt(b)

eq1 = p1 - p2
eq2 = 4b^2 - (a^2 + b^2);
res = solve([eq1, eq2], (a, b))[1]

   (q*(-1 + 3*sqrt(3)), q*(9 - sqrt(3))/3)

@syms q
side_a = 2q*(9 - √Sym(3))/3
height = √Sym(3)side_a/2
p = side_a*sqrt((height - 2*q)/height)/2 |> simplify
p |> println

   3^(3/4)*q*sqrt(-12 + 10*sqrt(3))/3

手計算で簡約化すると，以下のように若干簡単になる。

p = q*sqrt(10 - 4*sqrt(Sym(3)))
p |> println

   q*sqrt(10 - 4*sqrt(3))

長半径 p は短半径 q の sqrt(10 - 4√3) 倍である。
短径が 1 寸のとき，長径は sqrt(10 - 4√3) = 1.7526542071168778 寸である。

function draw(q, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (q*(-1 + 3√3), q*(9 - √3)/3)
   p = q*sqrt(10 - 4√3)
   @printf("a = %g;  b = %g;  q = %g;  p = %g\n", a, b, q, p)
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:blue, lw=0.5)
   segment(a, b, -a, -b, :magenta)
   segment(-a, b, a, -b, :magenta)
   ellipse(a - q, 0, q, p, color=:red)
   ellipse(q - a, 0, q, p, color=:red)
   ellipse(0, q - b, p, q, color=:green)
   ellipse(0, b - q, p, q, color=: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(a - q, 0, "(a-q,0)", :red, :center, delta=-delta/2)
       point(0, b - q, "(0,b-q,0)", :green, :center, delta=-delta/2)
       point(a, 0, " a", :red, :left, :bottom, delta=delta/2)
       point(0, b, " b", :blue, :left, :bottom, delta=delta/2)
       point(-0.3a, -0.5b, @sprintf(" a=%g, b=%g\np=%g, q=%g", a, b, p, q), :black, mark=false)
   end
end;

和算の心（その005）

和算の心（その005）

三角形に内接する楕円 -- 算法助術の公式97

1.　不等辺三角形の場合

不等辺三角形の場合について一般解が示されている。
和算の慣例として，楕円の長軸（長径），短軸（短径）を p, q としているが，ここでは長半径，短半径を使用することにする。

include("julia-source.txt");

@syms a, b, c, h, p, q
eq97 = -(b^2 - c^2)^2*h*q^2 + (b^2 - c^2)^2*q^3 + a^4*h*(2h - q)^2 - a^4*q*(2h - q)^2 - a^2*h*p^2*(2h - q)^2
eq97 = eq97(q => 2q, p => 2p)
eq_97 = eq97/4 |> simplify
eq_97 |> println

   a^4*h*(h - q)^2 - 2*a^4*q*(h - q)^2 - 4*a^2*h*p^2*(h - q)^2 - h*q^2*(b^2 - c^2)^2 + 2*q^3*(b^2 - c^2)^2

2.　二等辺三角形の場合

b = c の場合は，簡単になる。

eq97_2 = eq97(c => b) |> factor
eq97_2 |> println

   -4*a^2*(-h + q)^2*(-a^2*h + 2*a^2*q + 4*h*p^2)

更に a ≠ 0, h ≠ q なので非常に簡単な式になる。

eq97_3 = eq97_2/(4a^2*(q - h)^2)
eq97_3 |> println

   a^2*h - 2*a^2*q - 4*h*p^2

2.1.　a, h, q を与えて p を知るとき

res_p = solve(eq97_3, p)[2]
res_p |> println

   a*sqrt((h - 2*q)/h)/2

a = 5; h = 4; q = 1 のとき

ans_p = res_p(a => 5, h => 4, q => 1)
ans_p |> println
ans_p.evalf() |> println

   5*sqrt(2)/4
   1.76776695296637

2.2.　a, h, p を与えて q を知るとき

res_q = solve(eq97_3, q)[1]
res_q |> println

   h/2 - 2*h*p^2/a^2

a = 5; h = 4; q = 2 のとき

ans_q = res_q(a => 5, h => 4, p => 2)
ans_q |> println
ans_q.evalf() |> println

   18/25
   0.720000000000000

3.　関数として定義する

function formula97(a, h; p=0, q = 0)
   p == 0 && q == 0 && return "error"
   p != 0 && return h/2 - 2h*p^2/a^2
   q != 0 && return a*sqrt((h - 2q)/h)/2
end;

formula97(5, 4, p=1)

   1.68

formula97(5, 4, q=1)

   1.7677669529663689

formula97(5, 4)

   "error"

4.　図を描いて確かめる

function draw(a, h, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot([0, a, a/2, 0], [0, 0, h, 0], color=:blue, lw=0.5)
   delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
   q = 1
   p = formula97(a, h, q=1)
   @printf("p = %g;  q = %g\n", p, q)
   point(a/2, q, @sprintf("長半径 = %g\n短半径 = %g", p, q), :green, :center, delta=-delta/2)
   ellipse(a/2, q, p, q, color=:green)

   p = 1
   q = formula97(a, h, p=1)
   @printf("p = %g;  q = %g\n", p, q)
   point(a/2, q, @sprintf("長半径 = %g\n短半径 = %g", p, q), :red, :center, delta=-delta/2)
   ellipse(a/2, q, p, q, 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)
   end
end;

和算の心（その004）

和算の心（その004）

円に引かれた水平の弦と円弧に囲まれた（狭い方の）隙間に一番大きい円は 1 個ある。二番目以降の円は 2 個ずつある。外円と一番大きい円の直径が与えられたとき，2番目以降の円の直径を求めよ。

「和算の心（その003）」では特別な場合として，弦が円の中心を通る場合について述べた。
今回は，一般の場合について述べる。

include("julia-source.txt")

using SymPy

@syms R::positive, r1::positive,
     r2::positive, x2::positive,
     x::positive
y = R - 2r1
eq1 = x2^2 + (y + r2)^2 - (R - r2)^2
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq3 = x/(2R - 2r1) - 2r1/x
eq3 = x^2 + y^2 - R^2
res2 = solve([eq1, eq2, eq3], (r2, x2, x))[1]

   (-r1*(-R + r1)/R, 2*r1*sqrt(R - r1)/sqrt(R), 2*sqrt(r1)*sqrt(R - r1))

ans_r2 = res2[1](R=>3, r1=> 1/2)
ans_r2 |> println
ans_x2 = res2[2](R=>3, r1=> 1/2)
ans_x2 |> println

   0.416666666666667
   0.52704627669473*sqrt(3)

@syms r3::positive, x3::positive
eq4 = x3^2 + (y + r3)^2 - (R - r3)^2
eq5 = (x3 - x2)^2 + (r2 - r3)^2 - (r2 + r3)^2
res3 = solve([eq4, eq5], (r3, x3))[1]
res3 |> println
res3[1](R => 3, r1 => 1/2, r2 => ans_r2, x2 => ans_x2).evalf() |> println
res3[2](R => 3, r1 => 1/2, r2 => ans_r2, x2 => ans_x2).evalf() |> println

   (-(4*R*r1 - 4*r1^2 + x2^2 - 2*x2*(-sqrt(r2)*sqrt((-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2))/(-R + r1 - r2) + x2*(-R + r1)/(-R + r1 - r2)))/(4*(-R + r1 - r2)), -sqrt(r2)*sqrt((-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2))/(-R + r1 - r2) + x2*(-R + r1)/(-R + r1 - r2))
   0.255102040816327
   1.56492159287190

@syms r4::positive, x4::positive
eq6 = x4^2 + (y + r4)^2 - (R - r4)^2
eq7 = (x4 - x3)^2 + (r3 - r4)^2 - (r3 + r4)^2
res4 = solve([eq6, eq7], (r4, x4))[1]

   (-(4*R*r1 - 4*r1^2 + x3^2 - 2*x3*(-sqrt(r3)*sqrt((-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r3 + x3^2))/(-R + r1 - r3) + x3*(-R + r1)/(-R + r1 - r3)))/(4*(-R + r1 - r3)), -sqrt(r3)*sqrt((-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r3 + x3^2))/(-R + r1 - r3) + x3*(-R + r1)/(-R + r1 - r3))

r3, r4 の式を比較すれば，変数名の置き換えだけで次々と半径，と中心座標(x 座標)が求まることがわかる。
そこで，次のような関数を定義しておけば，簡単に計算を続けることができる。

R = 2, r1 = 1/2 のとき，円の半径は以下のようになる（分数は近似値）。

1: 0.5  1//2
2: 0.375  3//8
3: 0.18  9//50
4: 0.0688775510204081  27//392
5: 0.0240928019036288  79//3279
6: 0.00816312819134647  141//17273
7: 0.00273597297804469  2//731
8: 0.000913659014267574  2//2189
9: 0.00030473867949933  1//3281
10: 0.000101600202938299  1//9833

この場合もまた「算額の心（その003）」と同じく，簡単な一般項は求めがたい。

r1 = R/2 のときも，この式が使える。

nextcircle(r2, x2, R, r1) = (-(4*R*r1 - 4*r1^2 + x2^2 - 2*x2*(-sqrt(r2)*sqrt((-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2))/(-R + r1 - r2) + x2*(-R + r1)/(-R + r1 - r2)))/(4*(-R + r1 - r2)), -sqrt(r2)*sqrt((-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2))/(-R + r1 - r2) + x2*(-R + r1)/(-R + r1 - r2))

function draw(r1, R, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, x2, x) = (-r1*(-R + r1)/R, 2*r1*sqrt(R - r1)/sqrt(R), 2*sqrt(r1)*sqrt(R - r1))
   y = R - 2r1
   plot()
   circlef(0, 0, R, :aliceblue)
   circle(0, 0, R, :black)
   circlef(0, R - r1, r1, 1)
   @printf("%2d: %.15g  ", 1, r1)
   println(rationalize(r1, tol=1e-7))
   for i = 2:10
       @printf("%2d: %.15g  ", i, r2)
       println(rationalize(r2, tol=1e-7))
       circle2f(x2, y + r2, r2, i)
       (r2, x2) = nextcircle(r2, x2, R, r1)
   end
   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)
   end
end;

和算の心（その003）

和算の心（その003）

円に引かれた水平の弦と円弧に囲まれた（狭い方の）隙間に一番大きい円は 1 個ある。二番目以降の円は 2 個ずつある。一番大きい円の直径を 1 としたとき，2番目以降の円の直径を求めよ。

この問題を解く場合に，算法助法の公式29が役に立つ。

外円の中心を原点とする。
外円の半径を R とすれば，一番大きい円の半径 r1 は r1 = R/2 である。
弦の長さを length とすれば，length^2 = 4(2R)*(2r2) = 16R*r2 であるというのが公式29である。

ここでは特別な場合として，弦が円の中心を通る場合について述べる。この場合，length = 2R である。

include("julia-source.txt")

using SymPy

@syms R, r1, r2, length
formula29 = length^2 - 4(2R)*(2r2)
formula29 = formula29(length => 2R) |> simplify
formula29 |> println

   4*R*(R - 4*r2)

R ≠ 0 なので，R - 4r2 = 0 ということである。r2 について解けば， r2 = R/4 すなわち，2 番目に大きい円の半径は，外円の半径の 1/4 である。
r1 の半径は R/2 だったので，更にその 1/2 である。

公式29を直接使わなくても，以下の連立方程式を解けば，同じ結果が得られる。
このやり方だと，2番目に大きい円の中心座標（x 座標）も同時に得ることができる。また，3番目以降の円の半径を求めることもできる。

include("julia-source.txt")

using SymPy

@syms R::positive, r1::positive,
     r2::positive, x2::positive
R = 2r1
eq1 = x2^2 + r2^2 - (R - r2)^2
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
res2 = solve([eq1, eq2], (r2, x2))[1]

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

3 番目に大きい円の半径も同様に求めることができる。

@syms r3::positive, x3::positive
r2 = res2[1]
x2 = res2[2]
eq3 = x3^2 + r3^2 - (R - r3)^2
eq4 = (x3 - x2)^2 + (r2 - r3)^2 - (r2 + r3)^2
res3 = solve([eq3, eq4], (r3, x3))[1]

   (r1/9, 4*sqrt(2)*r1/3)

4 番目以降の円の半径は前式を初期値として，漸化式で表すことができる。

@syms r4::positive, x4::positive
r3 = res3[1]
x3 = res3[2]
eq5 = x4^2 + r4^2 - (R - r4)^2
eq6 = (x4 - x3)^2 + (r3 - r4)^2 - (r3 + r4)^2
res4 = solve([eq5, eq6], (r4, x4))[1]

   (r1/50, 7*sqrt(2)*r1/5)

@syms r5::positive, x5::positive
r4 = res4[1]
x4 = res4[2]
eq7 = x5^2 + r5^2 - (R - r5)^2
eq8 = (x5 - x4)^2 + (r4 - r5)^2 - (r4 + r5)^2
res5 = solve([eq7, eq8], (r5, x5))[1]

   (r1/289, 24*sqrt(2)*r1/17)

@syms r6::positive, x6::positive
r5 = res5[1]
x5 = res5[2]
eq9 = x6^2 + r6^2 - (R - r6)^2
eq10 = (x6 - x5)^2 + (r5 - r6)^2 - (r5 + r6)^2
res6 = solve([eq9, eq10], (r6, x6))[1]

   (r1/1682, 41*sqrt(2)*r1/29)

一番大きい円の半径から順に列挙すると，r1 = R/2 で，
r1, r1/2, r1/9, r1/50, r1/289, r1/1682, ...
となる。一見何らかの規則性がありそうに見えるが，一筋縄ではいかないようだ。いずれにしろ，指数曲線よりは早く減衰するようだ。そのせいかどうか，この場合の累円を求めよという算額はまだ見たことがない。

function draw(r1, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, x2) = (r1/2, sqrt(2)*r1)
   (r3, x3) = (r1/9, 4*sqrt(2)*r1/3)
   (r4, x4) = (r1/50, 7*sqrt(2)*r1/5)
   (r5, x5) = (r1/289, 24*sqrt(2)*r1/17)
   (r6, x6) = (r1/1682, 41*sqrt(2)*r1/29)
   R = 2r1
   @printf("""r1 = %.15g;
       r2 = %.15g;  x2 = %.15g;  
       r3 = %.15g;  x3 = %.15g;
       r4 = %.15g;  x4 = %.15g;
       r5 = %.15g;  x5 = %.15g;
       r6 = %.15g;  x6 = %.15g;\n""",
       r1, r2, x2, r3, x3, r4, x4, r5, x5, r6, x6)
   plot()
   circle(0, 0, R)
   circle22f(0, r1, r1, 1)
   circle4f(x2, r2, r2, 2)
   circle4f(x3, r3, r3, 3)
   circle4f(x4, r4, r4, 4)
   circle4f(x5, r5, r5, 5)
   segment(-R, 0, R, 0, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
   end
end;