算額（その310）
「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県伊賀市 永保寺 弘化4年(1847)
問題文3
正方形の中に順次，長方形，菱形，円を入れる。長方形の縦（長い方の辺）の外に円をいれ，長方形の横（短い方の辺）の外に順次正方形，正三角形を入れる。長方形の面積が 6281.8600 坪あり，長方形の横中径（横を直径とした外接円の半径）が 32.65 尺のとき，各変数の値を求めよ。
現代解法に示されている図とは印象がかなり違う。長方形はかなり太めである。また円の大きさは異なる。
以下の方程式を解く。
include("julia-source.txt");
using SymPy
@syms 長方形の面積::positive,
長方形の横を直径とした外接円の半径::positive;
長方形の面積 = Sym(6281.8600)
横中径 = Sym(32.65);
# 直平横（長方形の短辺）
直平横 = 横中径 * 2;
# 直平縦（長方形の長辺）
直平縦 = 長方形の面積 / 直平横;
# 菱面（菱形の一辺）
@syms 菱面::positive
eq1 = (直平横/2)^2 + (直平縦/2)^2 - 菱面^2
菱面 = solve(eq1)[1].evalf();
# 外円径
@syms 外円径::positive
eq2 = (外円径/2 + 外円径/2*√2) - 直平縦/2
外円径 = solve(eq2)[1];
# 方平面（正方形の一辺） 右上角の正方形
方平面 = 直平横/2 / √2;
# 三角面（正三角形の一辺）
@syms 三角面::positive
eq3 = 三角面*cos(PI/12) - 方平面
三角面 = solve(eq3)[1];
# 上ノ中径（縦を直径とした外接円の半径）
上ノ中径 = 直平縦 / 2;
# 外方面（外正方形の一辺）
外方面 = 直平横 / √2 + 直平縦 / √2;
@printf("直平横 = %.6f\n直平縦 = %.6f\n菱面 = %.6f\n外円径 = %.6f\n方平面 = %.6f\n三角面 = %.6f\n上ノ中径 = %.6f\n外方面 = %.6f\n",
直平横, 直平縦, 菱面, 外円径, 方平面, 三角面, 上ノ中径, 外方面)
直平横 = 65.300000
直平縦 = 96.200000
菱面 = 58.134607
外円径 = 39.847345
方平面 = 23.087036
三角面 = 23.901459
上ノ中径 = 48.100000
外方面 = 114.197745
中央の円の半径 = 27.014288
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
plot([0, 外方面, 外方面, 0, 0], [0, 0, 外方面, 外方面, 0], color=:black, lw=0.5)
a = 直平横/√2
plot!([a, 外方面, 外方面 - a, 0, a], [0, 外方面 - a, 外方面, a, 0], color=:blue, lw=0.5)
r = 外円径/2
circle(外方面 - r, r, r)
circle(r, 外方面 - r, r)
b = 上ノ中径/√2
plot!([方平面, 外方面 - b, 外方面 - 方平面, b, 方平面], [方平面, b, 外方面 - 方平面, 外方面 - b, 方平面], color=:green, lw=0.5)
rect(0, 0, 方平面, 方平面, :brown)
rect(外方面, 外方面, 外方面 - 方平面, 外方面 - 方平面, :brown)
c = 三角面*sin(pi/12)
plot!([外方面 - 方平面, 外方面, 外方面 - 方平面 + c, 外方面 - 方平面], [外方面 - 方平面, 外方面 - 方平面 + c, 外方面, 外方面 - 方平面], color=:magenta, lw=0.5)
plot!([方平面, 0, 方平面 - c, 方平面], [方平面, 方平面 - c, 0, 方平面], color=:magenta, lw=0.5)
# r2（真ん中の円の半径）
@syms r2::positive
eq4 = distance(方平面, 方平面, b, 外方面 - b, 外方面/2, 外方面/2) - r2^2
r2 = solve(eq4)[1].evalf()
@printf("中央の円の半径 = %.6f\n", r2)
circle(外方面/2, 外方面/2, r2, :orange)
if more
point(外方面, 0, "外方面", :black, :left, :bottom)
point(0, 外方面, "外方面", :black, :left, :bottom)
point(方平面, 0, "方平面", :black, :left, :bottom)
point(0, 方平面, "方平面", :black, :left, :bottom)
point(外方面 - r, r, "\n(外方面-r,r)", :red, :center)
point(外方面/2, 外方面/2, "\nr2,(外方面/2,外方面/2)", :orange, :center)
point(直平横/√2, 0, " 直平横/√2", :blue, :left, :bottom)
point(外方面-b, b, "\n(外方面-上ノ中径/√2,上ノ中径/√2)", :green, :center)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
plot!(xlims=(-5, 127))
else
plot!(showaxis=false)
end
end;
