# ¿Þ£³ s3 ¶ÊÌ̤ΠY-Z ÌÌ
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 = 50
  b = 14.7
  c = 6.7
  k1 = 1.202
#k1 = -k1
  k2 = 1.118
  q0 = PI
  N = 0
  Ns = 5
  Np = 50
  pmin = 0
  pmax = pmin + PI2
  # r1 ¶ÊÀþ
  for (p = pmin; p <= pmax; p += PI / Np) {
	xyz1(a, k1, p, p0)
	print y, z >"'"$TMP2"'"
  }
  
  # 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 - q0
		xyz3(a, b, c, k1, k2, p, q, t, 0)
		if (ex >= 0) {
			print "" >"'"$TMP1"'"
			continue
		}
		print y, z >"'"$TMP1"'"
	}
	printf "\n\n" >"'"$TMP1"'"
  }
  # q3 Àþ
  for (p = pmin; p <= pmax; p += PI / Np) {
	q = a / b * k1 / sqrt(k1^2 - 1) * sqrt(k2^2 - 1) * p - q0
	t1 = -PI;
	t2 = PI
	for (t = t1; t <= t2; t += PI / Ns / 10) {
		xyz3(a, b, c, k1, k2, p, q, t, 0)
		if (ex >= 0) {
			print "" >"'"$TMP1"'"
			continue
		}
		print y, z >"'"$TMP1"'"
	}
	printf "\n\n" >"'"$TMP1"'"
  }
  # ÎسÔÀþ
  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 - q0
	p1 = atan2(sqrt(k1^2 - 1), (k1 < 0) ? -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 = 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)))
	xyz3(a, b, c, k1, k2, p, q, t, 0)
	print y, z >"'"$TMP1"'"
  }
  printf "\n\n" >"'"$TMP1"'"
  for (p = pmin; p <= pmax; p += PI / Np / 10) {
	q = a / b * k1 / sqrt(k1^2 - 1) * sqrt(k2^2 - 1) * p - q0
	p1 = atan2(sqrt(k1^2 - 1), (k1 < 0) ? -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 = 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)))
	xyz3(a, b, c, k1, k2, p, q, t + PI, 0)
	print y, z >"'"$TMP1"'"
  }
  exit 0
}
END {
#  print "set size square"
  print "set size ratio 5"
  print "set nokey"
  print "set noxtics"
  print "set noytics"
  print "set noborder"
  print "set xzeroaxis"
  print "set yzeroaxis"
  print "plot \"'"$TMP1"'\" with lines lc \"red\", \"'"$TMP2"'\" with lines lc \"black\""
  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