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