算額（その800）
藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
底面，上面が菱形の角錐台の中に直径 6 寸の球がカッチリと入っている。上面の菱形の対角線の長いほうが 4 寸，短いほうが 3 寸のとき，底面の菱形の長いほうの直径はいかほどか。
上面の菱形の対角線 H，I (H > I)
底面の菱形の対角線 A，B (A > B)
球の半径を r とおく。
z 軸の上方から下を見ると下図のように見える。
赤の菱形が上面，青の菱形が底面である。赤い円が球を表す。これでは球は角錐台に接しているとは思えないが，次の図を見ると納得できる。
球は平面 ABIH に接している。接点は図に示す直線 Oa 上にある。Oa は AB, HI と直交する。
Oa を含み x-y 平面に垂直な面で切断して横から見たのが下図である。
直線 ah が上図における Oa である。
r = 3, b = 1.2 のとき，以下の方程式を解き，a = 15/2, h = 50/7 である。
底面と上面の相似比は a/b = 25/4 である。
したがって，上面の菱形の長い方の対角線が 4 寸のとき，底面上面の菱形の長い方の対角線は 25 寸である。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r, a, b, h
r = 3
b = 12//10
eq1 = (2r)^2 + (a - b)^2 - (a + b)^2
eq2 = (h - 2r)/b - h/a
(a, h) = solve([eq1, eq2], (a, h))[1]
(15/2, 50/7)
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r = 3
b = 1.2
a = 15/2
h = 50/7
@printf("相似比 = %g/%g = %g\n", a, b, a/b)
plot([0, a, 0, 0], [0, 0, h, 0], color=:black, lw=0.5)
segment(0, 2r, b, 2r)
circle(0, r, r)
if more == true
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(0, r, " r", :red, :left, :vcenter)
point(0, h, " h", :black, :left, :bottom, delta=delta)
point(b, 2r, "(b,2r)", :black, :left, :bottom, delta=delta)
point(0, 2r, " 2r", :black, :left, :bottom, delta=delta)
point(a, 0, " a", :black, :left, :bottom, delta=delta)
end
savefig("/Users/aoki/Downloads/fig0.png");
A = 25/2
B = A*3/4
H = 4/2
I = H*3/4
x0 = B*4/5*3/5
y0 = B*4/5*4/5
plot([A, 0, -A, 0, A], [0, B, 0, -B, 0], color=:blue, lw=0.5)
plot!([H, 0, -H, 0, H], [0, I, 0, -I, 0], color=:red, lw=0.5)
circle(0, 0, r, :red)
circle(0, 0, 7.5, :gray80)
segment(0, 0, x0, y0)
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(A, 0, " A", :blue, :left, :vcenter)
point(0, B, "B", :blue, :center, :bottom, delta=2delta)
point(H, 0, " H", :red, :left, :vcenter)
point(0, I, "I", :red, :center, :bottom, delta=2delta)
point(x0, y0, " a", :blue, :left, :bottom, delta=2delta)
point(0, 0, "O", :black, :right, delta=-delta)
end;
