# r2 曲面 - z 軸から見た x-y 平面の平行投影
# n 芯撚の X-Y 断面
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 = 1.8265	# 素線半径
  b = 1.0	# コア半径
  n = 4		# 芯数
  A1 = 1.0	# 撚込係数 sqrt(2) (撚角 PI/4)
  k1 = sqrt(1^2 + 1)	# A1 = 1.0
  # core
  for (p = 0; p <= PI2; p += PI/50) {
	xyz1k(a - b, k1, p, 0)
#	xyz1(a - b, A1, p, 0)
	print x, y >"'"$TMP1"'"
  }
  print "" >"'"$TMP1"'"
  for (p = 0; p <= PI2; p += PI/50) {
	xyz1k(a + b, k1, p, 0)
#	xyz1(a + b, A1, p, 0)
	print x, y >"'"$TMP1"'"
  }
  # z = 0 断面
  z1 = 0
  dp0 = PI2 / n
  p1 = atan2(1, A1)
  for (p0 = 0; p0 < PI2; p0 += dp0) {
	for (q = 0; q <= PI2; q += PI / 50) {
		p = -b/a*sin(p1)^2/cos(p1)*sin(q) + z1/a/A1 + p0
		xyz2k(a, b, k1, p, q, 0)
#		xyz2(a, b, A1, p, q, 0)
		print x, y >"'"$TMP2"'"
	}
	print "" >"'"$TMP2"'"
  }

  exit 0
}
END {
  print "set size square 0.5"
#  print "set size ratio 0.5 0.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 xyz1(a, A1, p, p0) {
  x = a * cos(p)
  y = a * sin(p)
  z = a * A1 * (p - p0)
}
function xyz2(a, b, A1, p, t, p0,  p1, x1, y1, z1) {
  p1 = atan2(1, A1)
  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)
  x1 = a * cos(p)
  y1 = a * sin(p)
  z1 = a * sin(p1) * A1 * (p - p0)
  x = x1 - b * ex
  y = y1 - b * ey
  z = z1 - 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