算額（その210）

算額（その210）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2（124）
長野県木曽郡植松町 臨川寺弁財天社 文政13年(1830)

外円内に，弦，斜，矢を引き，天円 1 個，地円 2 個を入れる。矢および斜が与えられたとき，天円の径を求めよ。

外円の半径を R，矢，斜の長さをそれぞれ a, b とする。
地円の半径 r1，中心座標を (r1, y1)
天円の半径 r2，中心座標を (0, R - r2)
弦と外円の交点座標を (x, a - R) 
として方程式を立てる。

y1 は負の値も取りうるので ::positive を指定しない。

using SymPy

@syms y1, r1::positive, R::positive, r2::positive, x::positive, a::positive, b::positive;

eq1 = r1^2 + y1^2 - (R - r1)^2
eq2 = r1^2 + (R - r2 - y1)^2 - (r1 + r2)^2
eq3 = x^2 + a^2 - b^2
eq4 = x^2 + (a - R)^2 - R^2
eq5 = y1 + (R - a) - r1

以上の連立方程式を a, b について解く。

res = solve([eq1, eq2, eq3, eq4, eq5], (y1, r1, r2, x, R))

   1-element Vector{NTuple{5, Sym}}:
    (b*(2*a - b)/(2*a), -a + b, (-a*b^2 + b^3)/(2*a^2 + 2*a*b), sqrt(-a + b)*sqrt(a + b), b^2/(2*a))

任意の a, b を式に与えれば解が得られる。

a, b = 3, 5
(b*(2*a - b)/(2*a), -a + b, (-a*b^2 + b^3)/(2*a^2 + 2*a*b), sqrt(-a + b)*sqrt(a + b)), b^2/(2*a))

   (0.8333333333333334, 2, 1.0416666666666667, 4.0, 4.166666666666667)

draw(4, 7, false)

   y1 = 0.87500;  r1 = 3.00000;  r2 = 1.67045;  x = 5.74456;  R = 6.12500

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=:blue, 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(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (y1, r1, r2, x, R) = (b*(2*a - b)/(2*a), -a + b, (-a*b^2 + b^3)/(2*a^2 + 2*a*b), sqrt(-a + b)*sqrt(a + b), b^2/(2*a))
   @printf("y1 = %.5f;  r1 = %.5f;  r2 = %.5f;  x = %.5f;  R = %.5f\n", y1, r1, r2, x, R)
   plot()
   circle(0, 0, R, :black)
   circle(r1, y1, r1)
   circle(-r1, y1, r1)
   circle(0, R - r2, r2, :blue)
   if more == true
       point(r1, y1, "  (r1,y1)  地") 
       point(0, R - r2, " R-r2  天") 
       point(0, -R, " -R")
       point(0, a - R, " a-R")
       point(x, a - R, "(x,a-R)", :blue, :right, :bottom)
       point(0, (a - 2R)/2, "a", :green, mark=false)
       point(x/2, (a - 2R)/2, "b", :green, mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
   plot!([0, x, -x, 0, 0], [-R, a - R, a - R, -R, a - R], color=:magenta, lw=0.5)
end;

算額（その209）

算額（その209）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2（123）
長野県千曲市八幡 八幡社 文政13年(1830)
長野県長野市安茂里正覚院 久保寺観音堂 文政13年(1830)

台形内に1本の対角線を引き，全円と小円を入れる。全円の径が 6 寸 5 分，上底の長さが 5寸のとき，小円の径はいかほどか。

方程式中の定数が整数になるように20倍して考える。座標軸を図のように定め，x が台形の右下の座標，小円の半径を r1，中心座標を (x1, 130 - r1) とし，連立方程式を解く。

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 x::positive, x1::negative, r1::positive;
#@syms x, x1, r1
eq1 = distance(-x, 0, -50, 130, x1, 130 - r1) - r1^2
eq2 = distance(-x, 0, 50, 130, x1, 130 - r1) - r1^2
eq3 = distance(-x, 0, -50, 130, 0, 65) - 65^2
res = solve([eq1, eq2, eq3], (x, x1, r1))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (169/2, -sqrt(139961)/4 + 72361*sqrt(212322 - 538*sqrt(139961))/270400 + 269*sqrt(29716799442 - 75299018*sqrt(139961))/270400, -269*sqrt(212322 - 538*sqrt(139961))/1040 - 69*sqrt(139961)/1040 + 86161/1040)
    (169/2, -sqrt(139961)/4 - 269*sqrt(29716799442 - 75299018*sqrt(139961))/270400 - 72361*sqrt(212322 - 538*sqrt(139961))/270400, -69*sqrt(139961)/1040 + 269*sqrt(212322 - 538*sqrt(139961))/1040 + 86161/1040)

x1 は式が長いが，二重根号を外して簡単にできる。

-sqrt(Sym(139961))/4 + 72361*sqrt(212322 - 538*sqrt(Sym(139961)))/270400 + 269*sqrt(29716799442 - 75299018*sqrt(Sym(139961)))/270400 |> sympy.sqrtdenest |> println

   sqrt(139961)*(-1/4 + 269*sqrt(139961)/559844)

r1 も同じく，二重根号を外して簡単にできる。

-269*sqrt(212322 - 538*sqrt(Sym(139961)))/1040 - 69*sqrt(Sym(139961))/1040 + 86161//1040 |> sympy.sqrtdenest |> simplify |> println

   6097/40 - 13*sqrt(139961)/40

6097/40 - 13*sqrt(139961)/40

   30.838073790805964

求まるのは 20 倍した半径である。直径はこの数値の 1/10 である。元の単位でいえば 3 寸 0 分 8 厘 3 毛 8 糸あまりあり。

   x = 84.50000,  x1 = -26.27840,  r1 = 30.83807

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=:blue, 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)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x, x1, r1) = res[1]
   @printf("%.5f,  %.5f,  %.5f\n", x, x1, r1)
   plot([x, 50, -50, -x, x], [0, 130, 130, 0, 0], color=:black, lw=0.5)
   circle(0, 65, 65)
   circle(x1, 130 - r1, r1, :blue)
   segment(-x, 0, 50, 130, :magenta)
   if more == true
       point(x, 0, " x", :blue, :left, :bottom)
       point(0, 65, " 65")
       point(0, 130, " 130", :blue, :left, :bottom)
       point(x1, 130 - r1, "(x1,130-r1)", :blue, :center, :top)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;