#! /usr/bin/python3
# Last edited on 2023-06-11 18:07:10 by stolfi

# Computes a snail shell like object as the union of balls 
# moving along a spine that is a conical exponential helix.

import rn
import rmxn
import sys
from math import sqrt, hypot, sin, cos, atan2, log, exp, floor, ceil, inf, nan, pi

def shell(n):
  # The exponential shell parameterized by number of turns.
  
  # Main dimensions:
  R0 =  25      # Radius of spine at {t0}.
  r0 =   5      # Radius of ball at {t0}.
  w0 =  35      # Height of spine at {t0}.
  alfa = log(2) # Growth factor over .
  
  f0 = exp(alfa*t0)
  A = R0/f0
  B = r0/f0
  C = w0/f0
  
  f1 = exp(alfa*t1)
  w1 = w0*f1/f0
  
  dt = 0.05 # Spacing of balls along the spline.
  t = t0
  while t <= t1:
    f = exp(alfa*t)   # Scaling factor
    R = A*f           # Distance from spinr to {Z}-axis.
    r = B*f           # Radius of ball before bending. 

    # Coordinates of center on spine before bending: */
    u = R*cos(2*pi*n*t)
    v = R*sin(2*pi*n*t)
    w = C*f           
    xb,yb,zb, rb, ox,oy,oz, ux,uy,uz, wx,wy,wz = bend(u,v,w,r,R,t,alfa,t1-t0)
    sys.stdout.write("%12.4f" % (t)) # 1
    sys.stdout.write("  %12.4f %12.4f %12.4f" % (xb,yb,zb)) # 2,3,4
    sys.stdout.write("  %12.4f" % rb) # 5
    sys.stdout.write("  %12.4f %12.4f %12.4f" % (ox,oy,oz)) # 6,7,8
    sys.stdout.write("  %12.4f %12.4f %12.4f" % (ux,uy,uz)) # 9,10,11
    sys.stdout.write("  %12.4f %12.4f %12.4f" % (wx,wy,wz)) # 12,13,14
    sys.stdout.write("\n")
    t = t + dt
  sys.stdout.flush()
  return
  # ----------------------------------------------------------------------

def bend(u,v,w,r,R,t,alfa,kt):
  s = t/kt
 
  # The hyperspine:
  oR = R+r
  o = (-oR*sin(s*pi), 0, -oR*cos(s*pi))
  
  # Derivative of the hyperspine:
  # o' = -oR'*(sin(a),0,cos(a)) - oR*a'*(cos(a),0,-sin(a))
  #    = oR*(-alfa*(sin(a),0,cos(a)) - pi/kt(cos(a),0,-sin(a))
  do = (-alfa*sin(s*pi)-pi/kt*cos(s*pi), 0, -alfa*cos(s*pi)+pi/kt*sin(s*pi))
  mdo = sqrt(do[0]*do[0] + do[1]*do[1] + do[2]*do[2])
  
  assert do[1] == 0
  ew = (do[0]/mdo, 0, do[2]/mdo)
  ev = (0, 1, 0)
  eu = (-do[2]/mdo, 0, do[0]/mdo) # {-cross(ev,ew)}.
  
  # m = (1 + u/oR)
  m = 1.0
  
  xb = o[0] + m*u*eu[0] + m*v*ev[0]
  yb = o[1] + m*u*eu[1] + m*v*ev[1]
  zb = o[2] + m*u*eu[2] + m*v*ev[2]
  
  rb = m*r
  
  return xb, yb, zb, rb, o[0], o[1], o[2], eu[0], eu[1], eu[2], ew[0], ew[1], ew[2]
  # ----------------------------------------------------------------------
  
shell(0, 3)