moebius.py

   1
   2 import numpy as np
   3 from numpy import cos, sin
   4
   5 from bezier import BezierSegment
   6
   7 tau = np.pi * 2
   8
   9 origin = np.array((0,0,0))
  10 unit_x = np.array((1,0,0))
  11 unit_y = np.array((0,1,0))
  12 unit_z = np.array((0,0,1))
  13
  14 class ParametricCircle:
  15   def __init__(pc, c, r0, r1):
  16     ''' circle centred on c
  17         with theta=0 point at c+r0
  18         and with theta=tau/4 point at c+r1 '''
  19     pc._c  = c
  20     pc._r0 = r0
  21     pc._r1 = r1
  22   def radius(pc, theta):
  23     return pc._r0 * cos(theta) + pc._r1 * sin(theta)
  24   def point(pc, theta):
  25     return pc._c + pc.radius(theta)
  26
  27 class Twirler(ParametricCircle):
  28   def __init__(tw, c, r0, r1, cycles, begin_zeta):
  29     ''' circle centred on c, etc.
  30         but with an orientation at each point, orthogonal to
  31           the circle
  32         the orientation circles round cycles times during the
  33           whole cycle
  34         begin_zeta is the angle from outwards at theta==0
  35           positive meaning in the direction of r0 x r1
  36     '''
  37     ParametricCircle.__init__(tw, c, r0, r1)
  38     tw._cycles = cycles
  39     tw._begin_zeta = begin_zeta
  40     tw._axis = np.cross(r0, r1)
  41   def dirn(tw, theta):
  42     zeta = tw._begin_zeta + theta * tw._cycles
  43     r = tw.radius(theta)
  44     return cos(zeta) * r + sin(zeta) * tw._axis
  45
  46 class Moebius:
  47   def __init__(m, n_u): # ix_u will be in [0, n_u>
  48     m.edge    = Twirler(origin,  unit_z, unit_x, 2, 0)
  49     m.midline = Twirler(-unit_z, unit_z, unit_y, 1, 0)
  50     m._beziers = [ m._bezier(u) for u in np.linspace(0, 1, n_u) ]
  51   def _bezier(m,u):
  52     theta = u * tau
  53     cp = [None] * 4
  54     cp[0] =               m.edge   .point(theta)
  55     cp[1] = cp[0] + 0.5 * m.edge   .dirn (theta)
  56     cp[3] =               m.midline.point(theta*2)
  57     cp[2] = cp[3] + 0.5 * m.midline.dirn (theta*2)
  58     return BezierSegment(cp)
  59   def point(m, ix_u, t):
  60     return m._beziers[ix_u].point_at_t(t)