# 図３ s3 曲面の X-Y 面 TMP="tmp\$\$" trap 'rm -f \$TMP; 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 = 50 b = 14.7 c = 6.7 k1 = 1.202 # A1 = 1.5 k2 = 1.118 # A2 = 2.0 pa = 0.566 N = 0 Ns = 5 Np = 50 pmin = 0 pmax = pmin + PI2 # s3 線 for (t = 0; t <= PI2; t += PI / Ns) { for (p = pmin; p <= pmax; p += PI / Np) { q = a / b * k1 / sqrt(k1^2 - 1) * sqrt(k2^2 - 1) * p xyz3(a, b, c, k1, k2, p, q, t, 0) if (ez >= 0) { print "" >"'"\$TMP"'" continue } print x, y >"'"\$TMP"'" } printf "\n\n" >"'"\$TMP"'" } # q3 線 for (p = pmin; p <= pmax; p += PI / Np) { q = a / b * k1 / sqrt(k1^2 - 1) * sqrt(k2^2 - 1) * p t1 = -PI t2 = PI for (t = t1; t <= t2; t += PI / Ns / 10) { xyz3(a, b, c, k1, k2, p, q, t, 0) if (ez >= 0) { print "" >"'"\$TMP"'" continue } print x, y >"'"\$TMP"'" } printf "\n\n" >"'"\$TMP"'" } # 輪郭線 zmin = zmax = 0 for (p = pmin; p <= pmax; p += PI / Np / 10) { q = a / b * k1 / sqrt(k1^2 - 1) * sqrt(k2^2 - 1) * p A2 = 1 / sqrt(k2^2 - 1) p1 = atan2(sqrt(k1^2 - 1), 1) h2 = 1 - b / a * sin(p1) * cos(q) F2 = b / a * cos(p1) * sin(p1) p2 = atan2(1 / A2 + F2, h2) t = atan2(cos(p)*cos(q)-cos(p1)*sin(p)*sin(q), \ sin(p1)*sin(p2)*sin(p) \ +cos(p2)*(cos(p)*sin(q)+cos(p1)*sin(p)*cos(q))) t = atan2(sin(p1)*sin(q), cos(p1)*sin(p2)-sin(p1)*cos(p2)*cos(q)) xyz3(a, b, c, k1, k2, p, q, t, 0) print x, y >"'"\$TMP"'" } printf "\n\n" >"'"\$TMP"'" for (p = pmin; p <= pmax; p += PI / Np / 10) { q = a / b * k1 / sqrt(k1^2 - 1) * sqrt(k2^2 - 1) * p p1 = atan2(sqrt(k1^2 - 1), 1) A2 = 1 / sqrt(k2^2 - 1) h2 = 1 - b / a * sin(p1) * cos(q) F2 = b / a * cos(p1) * sin(p1) p2 = atan2(1 / A2 + F2, h2) t = cos(p2) / h2 * (F2 * A2 + 1) * q t = atan2(cos(p)*cos(q)-cos(p1)*sin(p)*sin(q), \ sin(p1)*sin(p2)*sin(p) \ +cos(p2)*(cos(p)*sin(q)+cos(p1)*sin(p)*cos(q))) t = atan2(sin(p1)*sin(q), cos(p1)*sin(p2)-sin(p1)*cos(p2)*cos(q)) xyz3(a, b, c, k1, k2, p, q, t + PI, 0) print x, y >"'"\$TMP"'" } exit 0 } END { print "set size square 0.5" print "set nokey" print "set noxtics" print "set noytics" print "set noborder" print "plot \"'"\$TMP"'\" with lines lc \"red\"" print "" exit 0 } function xyz3(a, b, c, k1, k2, p, q, t, p0, x2, y2, z2, A2, h2, F2, p1, p2) { p1 = atan2(sqrt(k1^2 - 1), (k1 < 0) ? -1 : 1) xyz2(a, b, k1, p, q, p0) x2 = x y2 = y z2 = z A2 = 1 / sqrt(k2^2 - 1) h2 = 1 - b / a * sin(p1) * cos(q) F2 = b / a * cos(p1) * sin(p1) p2 = atan2(1 / A2 + F2, h2) ex = -(cos(p) * cos(q) - cos(p1) * sin(p) * sin(q)) * cos(t) \ + (sin(p1) * sin(p2) * sin(p) \ + cos(p2) * (cos(p) * sin(q) + cos(p1) * sin(p) * cos(q))) * sin(t) ey = -(sin(p) * cos(q) + cos(p1) * cos(p) * sin(q)) * cos(t) \ - (sin(p1) * sin(p2) * cos(p) \ - cos(p2) * (sin(p) * sin(q) - cos(p1) * cos(p) * cos(q))) * sin(t) ez = sin(p1) * sin(q) * cos(t) \ - (cos(p1) * sin(p2) \ - sin(p1) * cos(p2) * cos(q)) * sin(t) x = x2 - c * ex y = y2 - c * ey z = z2 - c * ez if (k1 < 0) { y = -y z = -z } } 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 if (k1 < 0) { y = -y z = -z } } function xyz1(a, k1, p, p0) { x = a * cos(p) y = a * sin(p) z = a * k1 / sqrt(k1^2 - 1) * (p - p0) if (k1 < 0) { y = -y z = -z } } ' | gnuplot -persist