# 図２ 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