moebius.py

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