# 図４ s2 曲面の Y-Z 面と x = 0 断面
TMP1="tmp1$$"
TMP2="tmp2$$"
trap 'rm -f $TMP1 $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
  k = sqrt(1 / 1.0 + 1)
  p0 = PI / 4
  N = 0
  Ns = 5
  Np = 20
  pmin = p0
  pmax = pmin + PI2 * 3 / 4
  # s2 線
  for (t = 0; t <= PI2; t += PI / Ns) {
	p1 = atan2(k * 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, k, p, t, p0)
			print y, z >"'"$TMP1"'"
		}
		p1 += PI2
		p2 += PI2
		print "" >"'"$TMP1"'"
	}
	printf "\n\n" >"'"$TMP1"'"
  }
  printf "\n\n" >"'"$TMP1"'"
  # q2 線
  for (p = pmin; p <= pmax; p += PI / Np) {
	t1 = atan2(k * cos(p), sin(p))
	t2 = t1 + PI
	for (t = t1; t <= t2; t += PI / Ns / 10) {
		xyz2(a, b, k, p, t, p0)
		print y, z >"'"$TMP1"'"
	}
	printf "\n\n" >"'"$TMP1"'"
  }
  # 上端面
  for (t = 0; t <= PI2; t += PI / Ns / 2) {
	xyz2(a, b, k, pmax, t, p0)
	print y, z >"'"$TMP1"'"
  }
  printf "\n\n" >"'"$TMP1"'"
  # 下端面
  for (t = 0; t <= PI2; t += PI / Ns / 2) {
	xyz2(a, b, k, pmin, t, p0)
	print y, z >"'"$TMP1"'"
  }
  printf "\n\n" >"'"$TMP1"'"
  # 輪郭線
  for (p = pmin; p <= pmax; p += PI / Np) {
	t = atan2(k * cos(p), sin(p))
	xyz2(a, b, k, p, t, p0)
	print y, z >"'"$TMP1"'"
  }
  printf "\n\n" >"'"$TMP1"'"
  for (p = pmin; p <= pmax; p += PI / Np) {
	t = atan2(k * cos(p), sin(p))
	xyz2(a, b, k, p, t + PI, p0)
	print y, z >"'"$TMP1"'"
  }
  # x = 0 面
  p1 = atan2(1, 1/sqrt(k^2 - 1))
  for (pc = PI/2; pc < PI2; pc += PI) {
    for (t = 0; t <= PI2; t += PI / Ns / 2) {
	p = atan2(b*cos(p1)*sin(t), a - b*cos(t)) + pc
	xyz2(a, b, k, p, t, p0)
	print y, z >"'"$TMP2"'"
    }
    print "" >"'"$TMP2"'"
  }
  exit 0
}
END {
  print "set size ratio 3"
  print "set nokey"
  print "set noxtics"
  print "set noytics"
  print "set noborder"
  print "plot \"'"$TMP1"'\" with lines lc \"red\", \"'"$TMP2"'\" with lines lc \"black\""
  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