算額（その55）

算額（その55）

宮城県角田市横倉 愛宕神社 明治15年(1882)1月
http://www.wasan.jp/miyagi/yokokuraatago.html

大円の中に中円が 2 個，小円が 8 個入っている。それぞれの径を求めよ。

大円の半径を 1 とし，以下のように記号を定め方程式を解く。

using SymPy
@syms r1::positive, r2::positive, r3::positive, x2::positive, y2::positive, y1::positive;
@syms x3::positive;

eq1 = 1//4 +(1//2 + r1)^2 - (1 - r1)^2
eq2 = (1 - r1)^2 + y1^2 - (1 + r1)^2
solve([eq1, eq2], (r1, y1))

   1-element Vector{Tuple{Sym, Sym}}:
    (1/6, sqrt(6)/3)

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   (r1, y1) = (1/6, sqrt(6)/3)
   println("r1 = $r1;  y1 = $y1")
   circle(0, 0, 1, :green)
   circle(1, 0, 1, :green, beginangle=120, endangle=240)
   circle(-1, 0, 1, :green, beginangle=300, endangle=420)
   circle(1/2, 1/2 + r1, r1, :magenta)
   circle(-1/2, 1/2 + r1, r1, :magenta)
   circle(1/2, -1/2 - r1, r1, :magenta)
   circle(-1/2, -1/2 - r1, r1, :magenta)
   circle(r1, y1, r1, :magenta)
   circle(-r1, y1, r1, :magenta)
   circle(r1, -y1, r1, :magenta)
   circle(-r1, -y1, r1, :magenta)
   circle(1/2, 0, 1/2, :blue)
   circle(-1/2, 0, 1/2, :blue)
   if more
       point(r1, y1, "(r1, y1)", :magenta, :center)
       point(1/2, 1/2+r1, "(1/2,1/2+r1)", :magenta, :center)
       point(1/2, 0, " 1/2", :blue)
       point(0, 0, " 0", :black)
       point(1, 0, " 1", :black)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   end
end;

   r1 = 0.16666666666666666;  y1 = 0.8164965809277259

 

算額（その54）

算額（その54）

宮城県角田市横倉 愛宕神社 明治15年(1882)1月
http://www.wasan.jp/miyagi/yokokuraatago.html

交差する 2 つの大円の中に中円 2 個，小円 3 個がある。それぞれの径を求めよ。

大円の半径を 1 とし，以下のように記号を定め方程式を解く。

using SymPy
@syms r1::positive, r2::positive, r3::positive, x2::positive, y2::positive, y1::positive;
@syms x3::positive;

r3 = 1
r1 = r3 - y1
y2 = r3//2
eq1 = x3^2 +(r3 - r1)^2 - (r1 + r3)^2
eq2 = (x3 - x2)^2 + (r3 - y2)^2 - (r3 - r2)^2
eq3 = x3^2 + r1^2 - (r3 - r1)^2
eq4 = (x3 + x2)^2 + (r3 - y2)^2 - (r3 + r2)^2;

solve([eq1, eq2, eq3, eq4], (r2, x2, y1, x3))

   1-element Vector{NTuple{4, Sym}}:
    (sqrt(6)/6, 1/2, 5/6, sqrt(6)/3)

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   (r2, x2, y1, x3) = (sqrt(6)/6, 1/2, 5/6, sqrt(6)/3)
   r3 = 1
   r1 = r3 - y1
   y2 = r3/2
   println("r1 = $r1;  r2 = $r2;  r3 = $r3")
   println("x2 = $x2;  y2 = $y2;  y1 = $y1")
   println("x3 = $x3")
   circle(x3, r3, 1, :green)
   circle(-x3, r3, 1, :green)
   circle(x2, y2, r2, :magenta)
   circle(-x2, y2, r2, :magenta)
   circle(0, r1, r1)
   circle(0, y1, r1)
   circle(0, y1+2r1, r1)
   if more
       point(0, r1, "r1", :red)
       point(x2, y2, "(x2,y2)", :magenta, :center)
       point(0, y1, "y1", :red)
       point(0, r3+r1, "r3+r1", :red)
       point(x3, r3, " (x3,r3)", :green)
       vline!([0], color=:black, lw=0.5)
   end
   hline!([0], color=:black, lw=0.5)
end;

   r1 = 0.16666666666666663;  r2 = 0.40824829046386296;  r3 = 1
   x2 = 0.5;  y2 = 0.5;  y1 = 0.8333333333333334
   x3 = 0.8164965809277259

 

算額（その53）

算額（その53）

宮城県角田市横倉 愛宕神社 明治15年(1882)1月
http://www.wasan.jp/miyagi/yokokuraatago.html

徳竹亜紀子，谷垣美保: 2021年度の算額調査，仙台高等専門学校名取キャンパス 研究紀要，第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf

大円内に甲円 2 個，乙円 1 個，丙半円 2 個がある。大円の直径が 15 寸のとき丙円の直径を求めよ。

大円の半径と中心座標を r0, (0, 0)
乙円の半径と中心座標を r1, (0, 0)
甲円の半径と中心座標を r2, (0, r0 - r2)
丙半円の半径と中心座標を r3, (r1 + r3, 0)
とし，以下のように記号を定め方程式を解く。

include("julia-source.txt")

using SymPy
@syms r0::positive, r1::positive, r2::positive, r3::positive;
eq1 = (r0 - r2)^2 + (r1 + r3)^2 - (r2 + r3)^2;
eq2 = r1 + 2r2 - r0;
eq3 = (r1 + r3)^2 + r3^2 - r0^2;

solve([eq1, eq2, eq3], (r1, r2, r3))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (r0/5, 2*r0/5, 3*r0/5)

丙円の半径は 大円の半径の 3/5 倍である。
したがって，大円の直径が 15 寸のとき，丙円の直径は 15*3/5 = 9 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 15
   (r1, r2, r3) = r0.*(1/5, 2/5, 3/5)
   @printf("r1 = %g;  r2 = %g;  r3 = %g\n", r1, r2, r3)
   plot()
   circle(0, 0, r0, :black)
   circle(0, 0, r1, :magenta)
   circle(0, r0 - r2, r2)
   circle(0, r2 - r0, r2)
   circle(r1 + r3, 0, r3, :blue, beginangle=90, endangle=270)
   plot!([r1 + r3, r1+r3], [-r3, r3], color=:blue, lw=0.5)
   circle(-r1 - r3, 0, r3, :blue, beginangle=270, endangle=450)
   plot!([-r1 - r3, -r1 - r3], [-r3, r3], color=:blue, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r0-r2, "甲円:r2 \n(0,r0-r2) ", :red, :right)
       point(0, 0, "乙円:r1 ", :magenta, :right, :bottom, delta=delta)
       point(0, r1, "r1 ", :magenta, :right)
       point(r1, 0, "r1 ", :red, :right)
       point(r1+r3, 0, "丙半円:r3 \n(r1+r3,0) ", :blue, :right, :top, delta=-delta)
       point(r1+r3, r3, " (r1+r3,r3)", :blue)
       point(0, r0, " r0", :black, :left, :bottom, delta=delta/2)
   end
end;

算額（その52）

算額（その52）

宮城県伊具郡丸森町 鹿島神社 弘化5年(1848)立春
http://www.wasan.jp/miyagi/kasima1.html

2つの大円が交差し，内部に6個の中円，4個の小円がある。それぞれの径を求めよ。

大円の半径を 1 とし，図のよう記号を定め方程式を解く。

using SymPy
@syms r1::positive, r2::positive, x0::positive, x1::positive;

eq1 = x1^2 + r1^2 - (1 - r1)^2;
eq2 = (x1 - x0 + 1)^2 + r1^2 - (1 + r1)^2;
eq3 = x1^2 + (1 - r2 - r1)^2 - (r1 + r2)^2;
eq4 = (x0 - 1)^2 + (1 - r2)^2 - (1 +r2)^2;

res = solve([eq1, eq2, eq3, eq4], (r1, r2, x1, x0))

   1-element Vector{NTuple{4, Sym}}:
    (2/5, 1/5, sqrt(5)/5, 1 - 2*sqrt(5)/5)

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   # pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   gr(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   p = plot(ylims=(-1.05, 1.05))
   (r1, r2, x1, x0) = (2/5, 1/5, sqrt(5)/5, 1 - 2*sqrt(5)/5)
   println("r1 = $r1;  r2 = $r2;  x1 = $x1;  x0 = $x0")
   circle(0, 0, 1, :black)
   circle(x0-1, 0, 1, :black)
   circle(x1, r1, r1, :magenta)
   circle(-x1, r1, r1, :magenta)
   circle(x1, -r1, r1, :magenta)
   circle(-x1, -r1, r1, :magenta)
   circle(x0-1-x1, r1, r1, :magenta)
   circle(x0-1-x1, -r1, r1, :magenta)

   circle(0, 1-r2, r2, :blue)
   circle(0, r2-1, r2, :blue)
   circle(x0-1, 1-r2, r2, :blue)
   circle(x0-1, r2-1, r2, :blue)
   if more
       point(x1, r1, "(x1, r1)", :magenta, :center)
       point(x0, 0, "x0", :black)
       point(x0-1, 0, "x0-1", :black)
       point(-1, 0, "-1 ", :black, :right)
       point(0, 1-r2, " r2", :blue)
       hline!([0], color=:black, lw=0.5)
       vline!([0, x0-1], color=:black, lw=0.5)
   end
   display(p)
end;

   r1 = 0.4;  r2 = 0.2;  x1 = 0.447213595499958;  x0 = 0.10557280900008403

算額（その51）

算額（その51）

岩手県一関市萩荘 達古袋八幡神社 弘化3年(1846)8月
http://www.wasan.jp/iwate/tatukohatiman2.html

外円の中に 4 種類の円が入っている。それぞれの径を求めよ。

外円の半径を 1 とし，以下のように記号を定め，方程式を解く。

using SymPy
@syms r1::positive, r2::positive, r3::positive, r4::positive;
@syms x4::positive, y4::positive;

x4 = (1 - r4)cos(PI/6)
y4 = (1 - r4)sin(PI/6)
eq1 = x4^2 + (y4 - r2 - r1)^2 - (r1 + r4)^2 |> expand
eq1 |> println

eq2 = x4^2 + y4^2 - (r2 + r4)^2 |>expand
eq2 |> println

eq3 = x4^2 + (1 - r3 - y4)^2 - (r3 + r4)^2 |> expand
eq3 |> println

eq4 = r2 + 2r1 + 2r3 -1
eq4 |> println

res = solve([eq1, eq2, eq3, eq4], (r1, r2, r3, r4))

   2*r1*r2 - r1*r4 - r1 + r2^2 + r2*r4 - r2 - 2*r4 + 1
   -r2^2 - 2*r2*r4 - 2*r4 + 1
   -3*r3*r4 - r3 - r4 + 1
   2*r1 + r2 + 2*r3 - 1

   1-element Vector{NTuple{4, Sym}}:
    (-11/18 + 5*sqrt(7)/18, 23/9 - 8*sqrt(7)/9, -1/6 + sqrt(7)/6, -7/9 + 4*sqrt(7)/9)

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   # pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   gr(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   p = plot()
   (r1, r2, r3, r4) = (-11/18 + 5*sqrt(7)/18, 23/9 - 8*sqrt(7)/9, -1/6 + sqrt(7)/6, -7/9 + 4*sqrt(7)/9)
   println("r1 = $r1;  r2 = $r2;  r3 = $r3;  r4 = $r4")
   x4 = (1 - r4)cosd(30)
   y4 = (1 - r4)sind(30)
   circle(0, 0, 1, :black)
   circle(0, 0, r2, :magenta)
   circle(0, r2+r1, r1)
   circle(0, 1-r3, r3, :blue)
   circle((r1+r2)cosd(30), -(r1+r2)sind(30), r1, :red)
   circle(-(r1+r2)cosd(30), -(r1+r2)sind(30), r1, :red)
   circle((1-r3)cosd(30), -(1-r3)sind(30), r3, :blue)
   circle(-(1-r3)cosd(30), -(1-r3)sind(30), r3, :blue)
   circle(-x4, y4, r4, :green)
   circle(x4, y4, r4, :green)
   circle(0, r4-1, r4, :green)
   if more
       point(0, r2, "r2 ", :blue, :right)
       point(0, r1+r2, "r1+r2 ", :red, :right)
       point(0, 1-r3, "1-r3 ", :red, :right)
       point(x4, y4, "(x4,y4)", :green, :center)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
   display(p)
end;

   r1 = 0.12381980862905284;  r2 = 0.20377661238703038;  r3 = 0.27429188517743175;  r4 = 0.3981116938064847

 

算額（その50）

算額（その50）

山形県米沢市小野川町 小野川薬師堂 明治10年(1877)4月
http://www.wasan.jp/yamagata/onogawa.html

四辺形の中に円が 4 個入っている。甲円，乙円の径を 25, 18 としたとき，丁円の径を求めよ。

一度にまとめて方程式を解こうとすると解けない（単に計算時間がかかっているだけか?）ので，分割して解く。

まず，甲円，乙円，丙円を同定しよう。記号を以下のように定める。
using SymPy
@syms r3::positive;  # 丙円の半径
@syms x1::positive, y2::positive;  # 甲円の中心 X 座標，乙円の中心 Y 座標

eq1 = (x1 - 18)^2 + (25 - y2)^2 - 43^2
eq2 = (x1 - r3)^2 + (25 - r3)^2 - (25 + r3)^2
eq3 = (18 - r3)^2 + (y2 - r3)^2 - (18 + r3)^2
res = solve([eq1, eq2, eq3], (r3, x1, y2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-43*sqrt(1609 - 1020*sqrt(2))/14 - 15*sqrt(3218 - 2040*sqrt(2))/7 + 30*sqrt(2) + 103/2, sqrt(1609 - 1020*sqrt(2))/2 + 15*sqrt(2) + 53/2, -sqrt(1609/4 - 255*sqrt(2)) + 15*sqrt(2) + 67/2)
    (15*sqrt(3218 - 2040*sqrt(2))/7 + 43*sqrt(1609 - 1020*sqrt(2))/14 + 30*sqrt(2) + 103/2, -sqrt(1609 - 1020*sqrt(2))/2 + 15*sqrt(2) + 53/2, sqrt(1609/4 - 255*sqrt(2)) + 15*sqrt(2) + 67/2)

一番目の解が適切である。

println("r3 = $(res[1][1].evalf())")
println("x1 = $(res[1][2].evalf())")
println("y2 = $(res[1][3].evalf())")

   r3 = 15.1902798290461
   x1 = 54.1649893584910
   y2 = 48.2614175127018

次に，甲円と乙円の四角形の斜辺との接点座標を (X1, Y1), (X2, Y2) と置き，これを求める。

@syms X1::positive, Y1::positive, X2::positive, Y2::positive
eq11 = (X2 - 18)^2 + (Y2 - y2)^2 - 18^2
eq12 = (X1 - x1)^2 + (Y1 - 25)^2 - 25^2
eq13 = (Y1 - 25)/(X1 - x1) - (Y2 - y2)/(X2 - 18)
eq14 = (Y1 - 25)/(X1 - x1) * (Y2 - Y1)/(X2 - X1) + 1;

res = solve([eq11, eq12, eq13, eq14], (X1, Y1, X2, Y2));

4 組の解が得られるが，4 番目のものだけが適切である。

names = ["X1", "Y1", "X2", "Y2"]
for i = 1:4
   println("$(names[i]) = $(res[4][i])")
end

   X1 = (x1^3 - 36*x1^2 + x1*y2^2 - 50*x1*y2 + 774*x1 + 25*y2*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) - 625*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) + 3150)/(x1^2 - 36*x1 + y2^2 - 50*y2 + 949)
   Y1 = 25*(x1^2 + x1*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) - 36*x1 + y2^2 - 43*y2 - 18*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) + 774)/(x1^2 - 36*x1 + y2^2 - 50*y2 + 949)
   X2 = 18*(x1^2 - 43*x1 + y2^2 + y2*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) - 50*y2 - 25*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) + 1075)/(x1^2 - 36*x1 + y2^2 - 50*y2 + 949)
   Y2 = (x1^2*y2 - 36*x1*y2 + 18*x1*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) + y2^3 - 50*y2^2 + 1075*y2 - 324*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) - 3150)/(x1^2 - 36*x1 + y2^2 - 50*y2 + 949)

最後に丁円の中心から四辺形の斜辺までの距離を求める。

SymPy の Point, Line を使い，直線までの距離を求めるが，
ここでも，SymPy の結果そのものを使うと solve() で結果を求められないので，いままでの結果を数値化して用いる。

function distance(x1, y1, x2, y2, x0, y0)
   p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
   l = sympy.Line(p1, p2)
   l.distance(sympy.Point(x0, y0))
end;

r3 = 15.1902798290461
x1 = 54.1649893584910
y2 = 48.2614175127018
X1 = (x1^3 - 36*x1^2 + x1*y2^2 - 50*x1*y2 + 774*x1 + 25*y2*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) - 625*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) + 3150)/(x1^2 - 36*x1 + y2^2 - 50*y2 + 949)
Y1 = 25*(x1^2 + x1*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) - 36*x1 + y2^2 - 43*y2 - 18*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) + 774)/(x1^2 - 36*x1 + y2^2 - 50*y2 + 949)
X2 = 18*(x1^2 - 43*x1 + y2^2 + y2*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) - 50*y2 - 25*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) + 1075)/(x1^2 - 36*x1 + y2^2 - 50*y2 + 949)
Y2 = (x1^2*y2 - 36*x1*y2 + 18*x1*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) + y2^3 - 50*y2^2 + 1075*y2 - 324*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) - 3150)/(x1^2 - 36*x1 + y2^2 - 50*y2 + 949)
println("X1 = $X1;  Y1 = $Y1;  X2 = $X2;  Y2 = $Y2")

   X1 = 64.08580382962593;  Y1 = 47.94727522450802;  X2 = 25.142986419217163;  Y2 = 64.78345567434756

@syms r4::positive
eq00 = distance(X1, Y1, X2, Y2, r4, r4) - r4;

solve(eq00) としても 'NotImplementedError' になり解を求めることができない。その原因は方程式中に長精度整数があるからのようだ。

eq00 |> println

   -r4 + sqrt((11095289145632785833731051*r4/21266540642722116155792165 - 11720727187202811618008367911717774033/425330812854442323115843300000000000)^2 + (25663886212278315131969678*r4/21266540642722116155792165 - 13555275798104253130313974902537196837/212665406427221161557921650000000000)^2)

長精度整数らしいがどうも精度に問題があるらしい。そこで，この整数どもを big"整数" として完璧な長精度整数にする。

eq01 = -r4 + sqrt((big"93910527328635730934895117268"*r4/big"179999999999999212286276676385" - big"3100132341015136856051026453265830392636969" /big"112499999999999507678922922740625000000000")^2 + (big"217219132900722712058838648021" *r4/big"179999999999999212286276676385" - big"28682963588788462898433683777226651524373897" /big"449999999999998030715691690962500000000000")^2);

eq01 |> println

   -r4 + 69.44170763477157031724744552587294112975879397332431357713551885112604322531*Abs(0.01893276580611329111115986244784178445485447675151318541196915535528237519342*r4 - 1)

これで，r4 = 30 が求められる。2通りの解が得られるが 220.643479932700 は不適切解である。

solve(eq01) |> println

   Sym[30.0000000000000, 220.643479932700]

なお，算額では「甲円と乙円の径を掛け，2倍して，平方根を取る」としている。

確かに sqrt(25 * 18 * 2) = 30 となる。なぜだ...

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   #pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   gr(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   y(x) = (Y2-Y1)/(X2-X1)*(x-X2)+Y2
   p = plot()
   (r3, x1, y2) = (-43*sqrt(1609 - 1020*sqrt(2))/14 - 15*sqrt(3218 - 2040*sqrt(2))/7 + 30*sqrt(2) + 103/2, sqrt(1609 - 1020*sqrt(2))/2 + 15*sqrt(2) + 53/2, -sqrt(1609/4 - 255*sqrt(2)) + 15*sqrt(2) + 67/2)
   X1 = (x1^3 - 36*x1^2 + x1*y2^2 - 50*x1*y2 + 774*x1 + 25*y2*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) - 625*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) + 3150)/(x1^2 - 36*x1 + y2^2 - 50*y2 + 949)
   Y1 = 25*(x1^2 + x1*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) - 36*x1 + y2^2 - 43*y2 - 18*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) + 774)/(x1^2 - 36*x1 + y2^2 - 50*y2 + 949)
   X2 = 18*(x1^2 - 43*x1 + y2^2 + y2*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) - 50*y2 - 25*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) + 1075)/(x1^2 - 36*x1 + y2^2 - 50*y2 + 949)
   Y2 = (x1^2*y2 - 36*x1*y2 + 18*x1*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) + y2^3 - 50*y2^2 + 1075*y2 - 324*sqrt(x1^2 - 36*x1 + y2^2 - 50*y2 + 900) - 3150)/(x1^2 - 36*x1 + y2^2 - 50*y2 + 949)
   println("r3 = $r3;  x1 = $x1;  y2 = $y2")
   println("X1 = $X1;  Y1 = $Y1;  X2 = $X2;  Y2 = $Y2")
   circle(x1, 25, 25)
   circle(18, y2, 18, :blue)
   circle(r3, r3, r3, :green)
   r4 = 30
   circle(r4, r4, r4, :magenta)
   if more
       point(r4, r4, "(r4,r4)", :red, :center)
       point(x1, 25, "(x1,25)", :red, :center)
       point(18, y2, "(18,y2)", :red, :center)
       point(r3, r3, "(r3,r3)", :red, :center)
       point(X1, Y1, "   (X1,Y1)", :blue)
       point(X2, Y2, "   (X2,Y2)", :blue)
   end
   plot!([0, x1+25, x1+25, 0, 0], [0, 0, y(x1+25), y(0), 0], color=:black, lw=0.5)
   display(p)
end;

   r3 = 15.190279829046077;  x1 = 54.16498935849101;  y2 = 48.26141751270184
   X1 = 64.08580382962599;  Y1 = 47.947275224508004;  X2 = 25.142986419217173;  Y2 = 64.78345567434759