# 図６密巻
# 密巻
# 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
  A1 = 1.0
A1 = 0.126
  k1 = sqrt(1 / 0.126^2 + 1)		# A1 = 0.126
  p0 = PI / 4
#p0 = 0
  N = 0
  Ns = 5
  Np = 50
  pmin = p0
  pmax = pmin + PI2 * 5 / 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) {
			xyz2k(a, b, k1, p, t, p0)
#			xyz2(a, b, A1, p, t, p0)
#			if (ex > 0) {
#				print "" >"'"$TMP1"'"
#				continue
#			}
			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))
	t1 = atan2(sqrt(A1^2 + 1)/A1*cos(p), sin(p))
	t2 = t1 + PI
	for (t = t1; t <= t2; t += PI / Ns / 10) {
		xyz2k(a, b, k1, p, t, p0)
#		xyz2(a, b, A1, p, t, p0)
#		if (ex >= 0) {
#			continue
#		}
		print y, z >"'"$TMP1"'"
	}
	printf "\n\n" >"'"$TMP1"'"
  }
  zmin = zmax = 0
  # 上端面
  xyz2(a, b, A1, pmax, t, p0)
  if (A1 * y < 0) {
	for (t = 0; t <= PI2; t += PI / Ns / 2) {
		xyz2k(a, b, k1, pmax, t, p0)
#		xyz2(a, b, A1, pmax, t, p0)
		print y, z >"'"$TMP1"'"
	}
  }
  printf "\n\n" >"'"$TMP1"'"
  # 下端面
  xyz2(a, b, A1, pmin, t, p0)
  if (A1 * y > 0) {
	for (t = 0; t <= PI2; t += PI / Ns / 2) {
		xyz2k(a, b, k1, pmin, t, p0)
#		xyz2(a, b, A1, 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))
	t = atan2(sqrt(A1^2 + 1)/A1*cos(p), sin(p))
	xyz2k(a, b, k1, p, t, p0)
#	xyz2(a, b, A1, 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))
	t = atan2(sqrt(A1^2 + 1)/A1*cos(p), sin(p))
	xyz2k(a, b, k1, p, t + PI, p0)
#	xyz2(a, b, A1, p, t + PI, p0)
	print y, z >"'"$TMP1"'"
  }
  # x = 0 面
  p1 = atan2(1, A1)
  p1 =atan2(sqrt(k1^2 - 1), 1)
  for (pc = PI/2; pc < PI2 * 3/2; pc += PI) {
    for (t = 0; t <= PI2; t += PI / Ns / 2) {
	p = atan2(b*cos(p1)*sin(t), a - b*cos(t)) + pc
	xyz2k(a, b, k1, p, t, p0)
#	xyz2(a, b, A1, p, t, p0)
	print y, z >"'"$TMP2"'"
    }
    print "" >"'"$TMP2"'"
  }
  exit 0
}
END {
  print "set size 0.5"
  print "unset key"
  print "unset xlabel"
  print "unset ylabel"
  print "set noytics"
  print "set noborder"
  print "unset xzeroaxis"
  print "unset yzeroaxis"
#  print "set xtics nomirror"
#  print "set ytics nomirror"
  print "set noxtics"
  print "set noytics"
  print "plot \"'"$TMP1"'\" with lines lc \"red\", \"'"$TMP2"'\" with lines lc \"black\""
  print ""
  exit 0
}
# r1(s0,q1)
function xyz1(a, A1, p, p0) {
  x = a * cos(p)
  y = a * sin(p)
  z = a * A1 * (p - p0)
}
# r2(s1,q2)
function xyz2(a, b, A1, p, q, p0,  p1, x1, y1, z1) {
  p1 = atan2(1, A1)
  xyz1(a, A1, p, p0)
  x1 = x
  y1 = y
  z1 = z
  ex = cos(p) * cos(q) - cos(p1) * sin(p) * sin(q)
  ey = sin(p) * cos(q) + cos(p1) * cos(p) * sin(q)
  ez = -sin(p1) * sin(q)
  x = x1 - b * ex
  y = y1 - b * ey
  z = z1 - b * ez
#  x = a * cos(p) - b * ex
#  y = a * sin(p) - b * ey
#  z = a * A1 * (p - p0) - b * ez
}
function xyz1k(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 xyz2k(a, b, k1, p, t, p0,  p1, x1, y1, z1) {
  p1 = atan2(sqrt(k1^2 - 1), (k1 < 0) ? -1 : 1)
  xyz1k(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
#cat $TMP1