算額（その769）

算額（その769）

神壁算法 天明5年乙巳 秋九月 第二術 關流四傅藤田権平定資六弟子 奧州一關 梶山主水次俊
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

大円，小円，甲円，乙円がある。上下の二線（上線，下線）と斜線を引く。
記述が重複するが正確を期すため，これらは以下のように接する。
大円は二線，斜線，小円，甲円，乙円に接する。
小円は二線，大円に接する。
甲円は上線，斜線，大円に接する。
乙円は下線，斜線，大円に接する。
大円，小円，甲円の直径がそれぞれ 245寸，196寸，25寸のとき，乙円の直径はいかほどか。

大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (0, r2)
甲円の半径と中心座標を r3, (x3, y3)
乙円の半径と中心座標を r4, (x4, r4)
上線は二点 (x02, y02), (x03, y03) を通る
下線は二点 (x02, y02), (x01, 0) を通る
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive, r2::positive,
     r3::positive, x3::positive, y3::positive,
     r4::positive, x4::positive,
     x01::positive, x02::positive, y02::positive,
     x03::positive, y03::positive
eq1 = dist(x02, y02, x03, y03, x1, r1) - r1^2
eq2 = dist(x02, y02, x03, y03, 0, r2) - r2^2
eq3 = dist(x02, y02, x03, y03, x3, y3) - r3^2
eq4 = dist(x01, 0, x02, y02, x1, r1) - r1^2
eq5 = dist(x01, 0, x02, y02, x3, y3) - r3^2
eq6 = dist(x01, 0, x02, y02, x4, r4) - r4^2
eq7 = sqrt((x02 - x03)^2 + (y02 - y03)^2) - x1
eq8 = x1^2 + (r1 - r2)^2 - (r1 + r2)^2
eq9 = (x1 - x3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq10 = (x1 - x4)^2 + (r1 - r4)^2 - (r1 + r4)^2;

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, x3, y3, r4, x4, x01, x02, y02, x03, y03) = u
   return [
       -r1^2 + (r1 - y02 - (-y02 + y03)*((r1 - y02)*(-y02 + y03) + (-x02 + x03)*(-x02 + x1))/((-x02 + x03)^2 + (-y02 + y03)^2))^2 + (-x02 + x1 - (-x02 + x03)*((r1 - y02)*(-y02 + y03) + (-x02 + x03)*(-x02 + x1))/((-x02 + x03)^2 + (-y02 + y03)^2))^2,  # eq1
       -r2^2 + (-x02 - (-x02 + x03)*(-x02*(-x02 + x03) + (r2 - y02)*(-y02 + y03))/((-x02 + x03)^2 + (-y02 + y03)^2))^2 + (r2 - y02 - (-y02 + y03)*(-x02*(-x02 + x03) + (r2 - y02)*(-y02 + y03))/((-x02 + x03)^2 + (-y02 + y03)^2))^2,  # eq2
       -r3^2 + (-x02 + x3 - (-x02 + x03)*((-x02 + x03)*(-x02 + x3) + (-y02 + y03)*(-y02 + y3))/((-x02 + x03)^2 + (-y02 + y03)^2))^2 + (-y02 + y3 - (-y02 + y03)*((-x02 + x03)*(-x02 + x3) + (-y02 + y03)*(-y02 + y3))/((-x02 + x03)^2 + (-y02 + y03)^2))^2,  # eq3
       -r1^2 + (r1 - y02*(r1*y02 + (-x01 + x02)*(-x01 + x1))/(y02^2 + (-x01 + x02)^2))^2 + (-x01 + x1 - (-x01 + x02)*(r1*y02 + (-x01 + x02)*(-x01 + x1))/(y02^2 + (-x01 + x02)^2))^2,  # eq4
       -r3^2 + (-y02*(y02*y3 + (-x01 + x02)*(-x01 + x3))/(y02^2 + (-x01 + x02)^2) + y3)^2 + (-x01 + x3 - (-x01 + x02)*(y02*y3 + (-x01 + x02)*(-x01 + x3))/(y02^2 + (-x01 + x02)^2))^2,  # eq5
       -r4^2 + (r4 - y02*(r4*y02 + (-x01 + x02)*(-x01 + x4))/(y02^2 + (-x01 + x02)^2))^2 + (-x01 + x4 - (-x01 + x02)*(r4*y02 + (-x01 + x02)*(-x01 + x4))/(y02^2 + (-x01 + x02)^2))^2,  # eq6
       -x1 + sqrt((x02 - x03)^2 + (y02 - y03)^2),  # eq7
       x1^2 + (r1 - r2)^2 - (r1 + r2)^2,  # eq8
       (-r1 + y3)^2 - (r1 + r3)^2 + (x1 - x3)^2,  # eq9
       (r1 - r4)^2 - (r1 + r4)^2 + (x1 - x4)^2,  # eq10
   ]
end;

(r1, r2, r3) = (245, 196, 25) .// 2
iniv = BigFloat[219.13, 118.51, 212.5, 24.5, 109.57, 82.18, 107.08, 222.73, -106.65, 174.33]
res = nls(H, ini=iniv)

   ([219.1346617949794, 118.51160280748886, 212.5, 24.5, 109.5673308974897, 82.17549817311728, 107.0771642861831, 222.72727272727272, -106.6467651187968, 174.33221099887766], true)

乙円の半径は 24.5寸 である（直径 49寸）。

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

r1 = 122.5;  r2 = 98;  r3 = 12.5;  x1 = 219.135;  x3 = 118.512;  y3 = 212.5;  r4 = 24.5;  x4 = 109.567;  x01 = 82.1755;  x02 = 107.077;  y02 = 222.727;  x03 = -106.647;  y03 = 174.332

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = (245, 196, 25) .// 2
   (x1, x3, y3, r4, x4, x01, x02, y02, x03, y03) = res[1]
   @printf("大円，小円，甲円の直径がそれぞれ %g寸, %g寸, %g寸 のとき，乙円の直径は %g寸\n", 2r1, 2r2, 2r3, 2r4)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  x1 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  x4 = %g;  x01 = %g;  x02 = %g;  y02 = %g;  x03 = %g;  y03 = %g\n", r1, r2, r3, x1, x3, y3, r4, x4, x01, x02, y02, x03, y03)
   plot()
   circle(x1, r1, r1)
   circle(0, r2, r2, :blue)
   circle(x3, y3, r3, :green)
   circle(x4, r4, r4, :magenta)
   segment(x01, 0, x02, y02)
   y04 = (y02 - y03)/(x02 - x03)*(1.25x1 - x03) + y03
   segment(x03, y03, 1.25x1, y04)
   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(x01, 0, "x01 ", :black, :right, :bottom, delta=delta)
       point(x02, y02, "(x02,y02)", :black, :right, :bottom, delta=delta)
       point(x03, y03, "(x03,y03)", :black, :left, :bottom, delta=4delta)
       point(x1, r1, " 大円:r1,(x1,r1)", :red, :center, delta=-delta)
       point(0, r2, " 小円:r2,(0,r2)", :blue, :center, delta=-delta)
       point(x3, y3, "   甲円:r3,(x3,y3)", :green, :left, :vcenter)
       point(x4, r4, "      乙円:r4,(x4,r4)", :magenta, :left, :vcenter)
       segment(x03, 0, 1.4x1, 0)
   end
end;

算額（その768）

算額（その768）

宮城県白石市小原字馬頭山 三瀧神社奉納算額（萬蔵稲荷神社所蔵） 明治8年
徳竹亜紀子，谷垣美保：宮城県白石市小原地区の算額調査，仙台高等専門学校名取キャンパス研究紀要，第57号，2021.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2021/04/No57_2.pdf

大球の頂きに甲球が載り，大球と甲球の両方に接し，かつ互いにも接し合う乙球を複数個添える。大球，甲球の直径がそれぞれ 39 寸，13 寸のとき，乙球が6個ある場合に，乙球の直径はいかほどか。

大球の半径と中心座標を r1, (0, 0, 0)
甲球の半径と中心座標を r2, (0, 0, r1 + r2)
x-z 平面にある乙球の個数を n，半径と中心座標を r3, (x3, 0, z3); x3 > 0
とおおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive,
     x3::positive, z3::positive, n::integer
(r1, r2) = (39, 13) .// 2
eq1 = x3^2 + z3^2 - (r1 + r3)^2
eq2 = x3^2 + (r1 + r2 - z3)^2 - (r2 + r3)^2
eq3 = x3*sind(Sym(180)/n) - r3
res = solve([eq1, eq2, eq3], (r3, x3, z3))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (-78*sin(pi/n)^2/(3*sin(pi/n)^2 - 4), -78*sin(pi/n)/(3*sin(pi/n)^2 - 4), (39*sin(pi/n)^2 - 156)/(6*sin(pi/n)^2 - 8))

t = sin(pi/n)，u = 4 - 3t^2 とすれば，乙球の半径と中心座標は (78t^2/u, 78t/u, (312 - 78t^2)/u) である。

問のように，大球，甲球の直径がそれぞれ 39 寸，13 寸，乙球が 6 個のときには，乙球の直径は 12 寸である。

function draw(which=1, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   n = 9
   (r1, r2) = (39, 13) .// 2
   t = sin(pi/n)
   u = 4 - 3t^2
   (r3, x3, z3) = (78t^2/u, 78t/u, (312 - 78t^2)/4u)
   plot()
   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
   if which == 1
       circle(0, 0, r1)
       circle(0, r1+r2, r2, :blue)
       circle(x3, z3, r3, :green)
       point(0, r1, "r1", :red, :center, delta=-0.5)
       point(0, 0, "大球:r1,(0,0,0)", :red, :center, delta=-0.5)
       point(0, r1 + r2, "甲球:r2\n(0,0,r1+r2)", :blue, :center, :bottom, delta=0.5)
       point(x3, z3, "    乙球:r3\n    (x3,0,z3)", :green, :left, :vcenter)
       point(r1, 0, " x軸", :black, :left, delta=-0.25, mark=false) 
       point(0, r1 + 2r2, "z軸", :black, :center, :bottom, delta=0.2, mark=false) 
   else
       circle(0, 0, x3, :gray80)
       rotate(x3, 0, x3*sind(180/n), angle=360/n, :green)
       point(x3, 0, "乙球:r3\n(x3,0,z3)", :green, :center, :bottom, delta=0.25)
       point(x3 + r3, 0, "x軸 ", :black, :right, delta=-0.25, mark=false) 
       point(0, x3 + r3, "y軸", :black, :center, :bottom, mark=false) 
   end
end;