算額（その110）

算額（その110）

愛媛県松山市 伊佐爾波神社 嘉永3年(1850)正月
http://www.wasan.jp/ehime/isaniwa7.html

正三角形の底辺に内接し，隣同士の円も接する3つの円が，その上にあるもう1つの円とも接するような場合，大円，中円，小円の半径を求めよ。

算額（その4）に似ている。

正三角形の一辺の長さを 2 とし，大円，中円，小円の半径を r1, r2, r3 とし，中円の中心座標を (x2, r2) とする。x2 = 1 - √3*r2 である。

using SymPy
@syms r1::positive, r2::positive, r3::positive;

x2 = 1 - √Sym(3)r2
eq1 = (r2 - r3)^2 + x2^2 - (r2 + r3)^2;
eq2 = sqrt(Sym(3)) - r1 - 2r3 - 2r1;
eq3 = (r1 + 2r3 - r2)^2 + x2^2 - (r1 + r2)^2;

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

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-20/33 - sqrt(-113 + 72*sqrt(3))/33 + 6*sqrt(3)/11, sqrt(-113 + 72*sqrt(3))/11 + 2*sqrt(-339 + 216*sqrt(3))/33 + 8*sqrt(3)/33 + 6/11, -sqrt(3)*(-100*sqrt(3)/33 - 10*sqrt(-113/1452 + 6*sqrt(3)/121) + 35/11)/10)
    (-20/33 + sqrt(-113 + 72*sqrt(3))/33 + 6*sqrt(3)/11, -2*sqrt(-339 + 216*sqrt(3))/33 - sqrt(-113 + 72*sqrt(3))/11 + 8*sqrt(3)/33 + 6/11, -sqrt(3)*(-100*sqrt(3)/33 + 10*sqrt(-113/1452 + 6*sqrt(3)/121) + 35/11)/10)

(-20/33 - sqrt(-113 + 72*sqrt(3))/33 + 6*sqrt(3)/11, sqrt(-113 + 72*sqrt(3))/11 + 2*sqrt(-339 + 216*sqrt(3))/33 + 8*sqrt(3)/33 + 6/11, -sqrt(3)*(-100*sqrt(3)/33 - 10*sqrt(-113/1452 + 6*sqrt(3)/121) + 35/11)/10)

   (0.23500815165114197, 1.6355839657205662, 0.5135131763077257)

(-20/33 + sqrt(-113 + 72*sqrt(3))/33 + 6*sqrt(3)/11, -2*sqrt(-339 + 216*sqrt(3))/33 - sqrt(-113 + 72*sqrt(3))/11 + 8*sqrt(3)/33 + 6/11, -sqrt(3)*(-100*sqrt(3)/33 + 10*sqrt(-113/1452 + 6*sqrt(3)/121) + 35/11)/10)

   (0.4423806081209667, 0.2951073349188893, 0.20245449160298862)

r1 > r2 > r3 ゆえ，2 番目の解が適切である。

using Plots
using Printf

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)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = res[2]
   (r1, r2, r3) = (r1.evalf(), r2.evalf(), r3.evalf())
   x2 = 1-√3 * r2
   plot([-1, 1, 0, -1], [0, 0, √3, 0])
   circle(0, r1 + 2r3, r1)
   circle(0, r3, r3, :blue)
   circle(x2, r2, r2, :green)
   circle(-x2, r2, r2, :green)
   if more == true
       @printf("r1 = %.5f,  r2 = %.5f,  r3 = %.5f, x2 = %.5f", r1, r2, r3, x2)
       point(0, r1 + 2r3, "r1+2r3 ", :red, :right)
       point(0, r3, "r3 ", :blue, :right)
       point(x2, r2, "(x2,r2)", :green, :top)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   r1 = 0.44238,  r2 = 0.29511,  r3 = 0.20245, x2 = 0.48886

算額（その109）

算額（その109）

東京都府中市宮町 大國魂神社 明治18年
大國魂神社の算額（下段第5問）についての考察
http://www2.ttcn.ne.jp/~nagai/waseda/wasan/ok1.pdf

外円内に三角形，甲円，乙円，丙円が入っている。甲円，乙円，丙円の径をそれぞれ 49 寸，28 寸，17 寸のとき，外円の径を求めよ。

図のように記号を定め，方程式を立てる。
三角形の頂点 A, B, C の座標をそれぞれ (-x0, 98 - R), (x0, 98 - R), (x4, y4) とする。
甲円，乙円，丙円の中心座標と半径をそれぞれ (0, 49 -R; R), (x2, 41, 28), (x3, y3, 17) とする。

方程式は solve() では解けないので，nlsolve() で数値解を求める。

using SymPy

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))^2  # 二乗距離を返す！！
end;

@syms R, x2, x3, y3, x0, x4, y4;

eq1 = x3^2 + y3^2 - (R - 17)^2 |> expand;
eq2 = R^2 - (98 - R)^2 - x0^2 |> expand;
eq3 = distance(x0, 98 - R, -x0, 98 - R, x2, 41) - 28^2 |> expand;
eq4 = x4^2 + y4^2 - R^2 |> expand;
eq5 = distance(-x0, 98 - R, x4, y4, x2, 41) - 28^2;
eq6 = distance( x0, 98 - R, x4, y4, x2, 41) - 28^2;
eq7 = distance(-x0, 98 - R, x4, y4, x3, y3) - 17^2;

eq1 |> println
eq2 |> println
eq3 |> println
eq4 |> println
eq5 |> println
eq6 |> println
eq7 |> println

   -R^2 + 34*R + x3^2 + y3^2 - 289
   196*R - x0^2 - 9604
   R^2 - 114*R + 2465
   -R^2 + x4^2 + y4^2
   (41 - (41*R^2 + 82*R*y4 - 8036*R + 41*x0^2 + 82*x0*x4 + 41*x4^2 + 41*y4^2 - 8036*y4 + (x0 + x4)*(R*x2 - R*x4 + x0*y4 - 41*x0 + x2*y4 - 98*x2 + 57*x4) + 393764)/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 + (x2 - (x2*(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) - (R + y4 - 98)*(R*x2 - R*x4 + x0*y4 - 41*x0 + x2*y4 - 98*x2 + 57*x4))/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 - 784
   (x2 - (R^2*x4 + R*x0*y4 - 41*R*x0 + R*x4*y4 - 155*R*x4 + x0^2*x2 - 2*x0*x2*x4 + x0*y4^2 - 139*x0*y4 + 4018*x0 + x2*x4^2 - 57*x4*y4 + 5586*x4)/(R^2 + 2*R*y4 - 196*R + x0^2 - 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 + (-((98 - R)*(R^2 + 2*R*y4 - 196*R + x0^2 - 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) + (R + y4 - 98)*(R^2 + R*y4 - 155*R + x0^2 - x0*x2 - x0*x4 + x2*x4 - 57*y4 + 5586))/(R^2 + 2*R*y4 - 196*R + x0^2 - 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) + 41)^2 - 784
   (x3 - (x3*(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) - (R + y4 - 98)*(R*x3 - R*x4 - x0*y3 + x0*y4 + x3*y4 - 98*x3 - x4*y3 + 98*x4))/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 + (y3 - (y3*(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) + (x0 + x4)*(R*x3 - R*x4 - x0*y3 + x0*y4 + x3*y4 - 98*x3 - x4*y3 + 98*x4))/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 - 289

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10])

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=1e-14)
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
v = r.zero
end
return v, r.f_converged
end;

function H(u)
   (R, x2, x3, y3, x0, x4, y4) = u
   return [
       -R^2 + 34*R + x3^2 + y3^2 - 289,
       196*R - x0^2 - 9604,
       R^2 - 114*R + 2465,
       -R^2 + x4^2 + y4^2,
       (41 - (41*R^2 + 82*R*y4 - 8036*R + 41*x0^2 + 82*x0*x4 + 41*x4^2 + 41*y4^2 - 8036*y4 + (x0 + x4)*(R*x2 - R*x4 + x0*y4 - 41*x0 + x2*y4 - 98*x2 + 57*x4) + 393764)/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 + (x2 - (x2*(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) - (R + y4 - 98)*(R*x2 - R*x4 + x0*y4 - 41*x0 + x2*y4 - 98*x2 + 57*x4))/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 - 784,
       (41 - (41*R^2 + 82*R*y4 - 8036*R + 41*x0^2 - 82*x0*x4 + 41*x4^2 + 41*y4^2 - 8036*y4 - (x0 - x4)*(R*x2 - R*x4 - x0*y4 + 41*x0 + x2*y4 - 98*x2 + 57*x4) + 393764)/(R^2 + 2*R*y4 - 196*R + x0^2 - 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 + (x2 - (x2*(R^2 + 2*R*y4 - 196*R + x0^2 - 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) - (R + y4 - 98)*(R*x2 - R*x4 - x0*y4 + 41*x0 + x2*y4 - 98*x2 + 57*x4))/(R^2 + 2*R*y4 - 196*R + x0^2 - 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 - 784,
       (x3 - (-x0*(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) + (x0 + x4)*(R^2 + R*y3 + R*y4 - 196*R + x0^2 + x0*x3 + x0*x4 + x3*x4 + y3*y4 - 98*y3 - 98*y4 + 9604))/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 + (y3 - ((98 - R)*(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) + (R + y4 - 98)*(R^2 + R*y3 + R*y4 - 196*R + x0^2 + x0*x3 + x0*x4 + x3*x4 + y3*y4 - 98*y3 - 98*y4 + 9604))/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 - 289,
   ]
end;

iniv = [big"80.0", 30, -30, 55, 80, 30, 75]
res = nls(H, ini=iniv)
println(res)

   (BigFloat[85.0, 28.00000000000000000000000000000000000000000000000000000000000000000000000000055, -31.99999993545069734104462106203685900316672771867706686845499772119474317529376, 60.00000003442629475144286876691367519831107855003889767015735701621077260821287, 84.0, 35.99999999999999999999999999999999999999999999999999999999999999999999999999945, 77.0], true)

外円の半径は 85.0 である。元の単位での径（直径）は 85 寸である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
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 segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, x2, x3, y3, x0, x4, y4) = res[1]
   println("R = $R")
   plot()
   circle(0, 0, R)
   circle(0, 49-R, 49, :blue)
   circle(x2, 41, 28, :magenta)
   circle(x3, y3, 17, :green)
   plot!([x0, x4, -x0, x0], [98 - R, y4, 98 - R, 98 - R], color=:brown, lw=0.5)
    if more
       @printf("x0 = %.3f, x2 = %.3f, x3 = %.3f, y3 = %.3f, x4 = %.3f, y4 = %.3f", x0, x2, x3, y3, x4, y4)
       point(0, 49 - R, " 甲(0,49-R;49)", :blue, :left)
       point(x2, 41, "乙(x2,41;28)", :magenta, :top)
       point(x3, y3, "丙(x3,y3;17)", :green, :top)
       point(0, 0, " O")
       point(-x0, 98 - R, " A(-x0,98-R)", :brown)
       point(x0, 98 - R, "B(x0,98-R) ", :brown, :right)
       point(x4, y4, " C(x4,y4)", :brown, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
    end
end;

 

 

算額（その108）

算額（その108）

東京都府中市宮町 大國魂神社 明治18年
大國魂神社の算額（上段第2問）についての考察
http://www2.ttcn.ne.jp/~nagai/waseda/wasan/ok2.pdf

外円内に甲円，乙円，丙円が内接している。更に弦にも接している。甲円，乙円，丙円の径をそれぞれ 4 寸，13 寸，8 寸 4 分 5 厘のとき，外円の径を求めよ。

図のように記号を定め，方程式を立てる。
外円の半径を R とする。甲円，乙円，丙円の中心座標と半径を (0, 400 - R, 400)，(x4, y4, 1300), (x, y, 845) とする。

using SymPy

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))^2  # 二乗距離を返す！！
end;

@syms x1, x3, y3, x4, y4, x, y, R, r, x0, y0;

eq1 = x^2 + y^2 - (R - 845)^2
eq2 = distance(x1, (800 - R), x3, y3, x4, y4) - 1300^2
eq3 = (x - x0)^2 + (y - y0)^2 - 845^2
eq4 = distance(x1, (800 - R), x3, y3, x, y) - 845^2;
eq5 = R + y4 - 2100
eq6 = x4^2 + y4^2 - (R - 1300)^2
eq7 = x3^2 + y3^2 - R^2
eq8 = x1^2 + (R - 800)^2 - R^2
eq9 = x0 - (x1 + x3) / 2
eq10 = y0 - ((800 - R) + y3) / 2;

方程式は solve() では解けないので，nlsolve() で数値解を求める。

eq1 |> println
eq2 |>  println
eq3 |> println
eq4 |> println
eq5 |> println
eq6 |> println
eq7 |> println
eq8 |> println
eq9 |> println
eq10 |> println

   x^2 + y^2 - (R - 845)^2
   (x4 - (x1*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) - (x1 - x3)*(R^2 + R*y3 + R*y4 - 1600*R + x1^2 - x1*x3 - x1*x4 + x3*x4 + y3*y4 - 800*y3 - 800*y4 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 + (y4 - ((800 - R)*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) + (R + y3 - 800)*(R^2 + R*y3 + R*y4 - 1600*R + x1^2 - x1*x3 - x1*x4 + x3*x4 + y3*y4 - 800*y3 - 800*y4 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 - 1690000
   (x - x0)^2 + (y - y0)^2 - 714025
   (x - (x1*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) - (x1 - x3)*(R^2 + R*y + R*y3 - 1600*R - x*x1 + x*x3 + x1^2 - x1*x3 + y*y3 - 800*y - 800*y3 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 + (y - ((800 - R)*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) + (R + y3 - 800)*(R^2 + R*y + R*y3 - 1600*R - x*x1 + x*x3 + x1^2 - x1*x3 + y*y3 - 800*y - 800*y3 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 - 714025
   R + y4 - 2100
   x4^2 + y4^2 - (R - 1300)^2
   -R^2 + x3^2 + y3^2
   -R^2 + x1^2 + (R - 800)^2
   x0 - x1/2 - x3/2
   R/2 + y0 - y3/2 - 400

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10])

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=1e-14)
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
v = r.zero
end
return v, r.f_converged
end;

function H(u)
   (x1, x3, y3, x, y, R, x4, y4, x0, y0) = u
   return [
       x^2 + y^2 - (R - 845)^2,
       (x4 - (x1*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) - (x1 - x3)*(R^2 + R*y3 + R*y4 - 1600*R + x1^2 - x1*x3 - x1*x4 + x3*x4 + y3*y4 - 800*y3 - 800*y4 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 + (y4 - ((800 - R)*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) + (R + y3 - 800)*(R^2 + R*y3 + R*y4 - 1600*R + x1^2 - x1*x3 - x1*x4 + x3*x4 + y3*y4 - 800*y3 - 800*y4 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 - 1690000,
       (x - x0)^2 + (y - y0)^2 - 714025,
       (x - (x1*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) - (x1 - x3)*(R^2 + R*y + R*y3 - 1600*R - x*x1 + x*x3 + x1^2 - x1*x3 + y*y3 - 800*y - 800*y3 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 + (y - ((800 - R)*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) + (R + y3 - 800)*(R^2 + R*y + R*y3 - 1600*R - x*x1 + x*x3 + x1^2 - x1*x3 + y*y3 - 800*y - 800*y3 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 - 714025,
       R + y4 - 2100,
       x4^2 + y4^2 - (R - 1300)^2,
       -R^2 + x3^2 + y3^2,
       -R^2 + x1^2 + (R - 800)^2,
       x0 - x1/2 - x3/2,
       R/2 + y0 - y3/2 - 400,
   ]
end;

iniv = [-1500.0, 900, 2000, -1000, 800, 2200, 900, -100, -400, 300]
res = nls(H, ini=iniv)
println(res)

   ([-1700.0, 873.9999999999999, 2025.75, -1089.0000042590011, 816.7499943213318, 2206.25, 900.0, -106.24999999999997, -413.00000000000006, 309.74999999999994], true)

外円の半径は 2206.25 である。元の単位での径（直径）は 22 寸 0 分 6 厘 2 毛 5 糸である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
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 segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x1, x3, y3, x, y, R, x4, y4, x0, y0) = res[1]
   println("R = $R")
   plot()
   circle(0, 0, R)
   circle(x, y, 845, :blue)
   circle(x4, y4, 1300, :magenta)
   circle(0, 400-R, 400, :green)
   segment(x1, 800 - R, -x1, 800 - R)
   segment(x1, 800 - R, x3, y3)
    if more
       point(x4, y4, "(x4,y4)")
       point(0, 400-R, "400-R")
       point(x, y, "(x,y)")
       point(x0, y0, "(x0,y0)")
       point(x1, 800 - R, "(x1,800-R)")
       point(x3, y3, "(x3,y3)")
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
    end
end;