# 図２ s2 曲面の Y-Z 面 # 曲率線で面を表示 TMP="tmp\$\$" trap 'rm -f \$TMP \$TMP2; exit 0' 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 /usr/bin/awk ' BEGIN { PI = atan2(0, -1) PI2 = atan2(0, -1) * 2 a = 25 b = 10 k1 = sqrt(1 / 1.0 + 1) #k1 = -k1 p0 = PI / 3 N = 0 Ns = 5 Np = 20 pmin = p0 pmax = pmin + PI2 * 3 / 4 # s2 線 for (t = 0; t <= PI2; t += PI / Ns) { p1 = atan2(k1 * cos(t), sin(t)) p2 = p1 + PI while (p1 < pmax) { p3 = (pmin < p1) ? p1 : pmin p4 = (pmax < p2) ? pmax : p2 dp = (p4 - p3) / 50 for (p = p3; p <= p4 * 1.001; p += dp) { xyz2(a, b, k1, p, t, p0) print y, z >"'"\$TMP"'" } p1 += PI2 p2 += PI2 print "" >"'"\$TMP"'" } printf "\n\n" >"'"\$TMP"'" } printf "\n\n" >"'"\$TMP"'" # q2 線 for (p = pmin; p <= pmax; p += PI / Np) { t1 = atan2(k1 * cos(p), sin(p)) t2 = t1 + PI for (t = t1; t <= t2; t += PI / 50) { xyz2(a, b, k1, p, t, p0) print y, z >"'"\$TMP"'" } printf "\n\n" >"'"\$TMP"'" } zmin = zmax = 0 # 上端面 xyz2(a, 0, k1, pmax, 0, p0) if (k1 * y < 0) { for (t = 0; t <= PI2; t += PI / Ns / 2) { xyz2(a, b, k1, pmax, t, p0) print y, z >"'"\$TMP"'" } } printf "\n\n" >"'"\$TMP"'" # 下端面 xyz2(a, 0, k1, pmax, 0, p0) if (k1 * y < 0) { for (t = 0; t <= PI2; t += PI / Ns / 2) { xyz2(a, b, k1, pmin, t, p0) print y, z >"'"\$TMP"'" } } printf "\n\n" >"'"\$TMP"'" # 輪郭線 for (p = pmin; p <= pmax; p += PI / Np) { t = atan2(k1 * cos(p), sin(p)) xyz2(a, b, k1, p, t, p0) print y, z >"'"\$TMP"'" } printf "\n\n" >"'"\$TMP"'" for (p = pmin; p <= pmax; p += PI / Np) { t = atan2(k1 * cos(p), sin(p)) xyz2(a, b, k1, p, t + PI, p0) print y, z >"'"\$TMP"'" } exit 0 } END { print "set size square 0.5" print "set nokey" print "set noxtics" print "set noytics" print "set noborder" print "set xzeroaxis" print "set yzeroaxis" print "plot \"'"\$TMP"'\" with lines lc \"red\"" print "" exit 0 } function xyz1(a, k1, p, p0) { x = a * cos(p) y = a * sin(p) z = a / sqrt(k1^2 - 1) * (p - p0) if (k1 < 0) y = -y } function xyz2(a, b, k1, p, t, p0, p1, x1, y1, z1) { p1 = atan2(sqrt(k1^2 - 1), (k1 < 0) ? -1 : 1) xyz1(a, k1, p, p0) x1 = x y1 = y z1 = z ex = cos(p) * cos(t) - cos(p1) * sin(p) * sin(t) ey = sin(p) * cos(t) + cos(p1) * cos(p) * sin(t) ez = -sin(p1) * sin(t) x = x1 - b * ex y = y1 - b * ey z = z1 - b * ez } ' | gnuplot -persist