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     return pc._c + pc.radius(theta)
  30
  31 class Twirler(ParametricCircle):
  32   def __init__(tw, c, r0, r1, cycles, begin_zeta):
  33     ''' circle centred on c, etc.
  34         but with an orientation at each point, orthogonal to
  35           the circle
  36         the orientation circles round cycles times during the
  37           whole cycle
  38         begin_zeta is the angle from outwards at theta==0
  39           positive meaning in the direction of r0 x r1
  40     '''
  41     ParametricCircle.__init__(tw, c, r0, r1)
  42     tw._cycles = cycles
  43     tw._begin_zeta = begin_zeta
  44     tw._axis = np.cross(r0, r1)
  45   def dirn(tw, theta):
  46     zeta = tw._begin_zeta + theta * tw._cycles
  47     r = tw.radius(theta)
  48     return cos(zeta) * r + sin(zeta) * tw._axis
  49
  50 class Moebius:
  51   def __init__(m, n_u): # ix_u will be in [0, n_u>
  52     m.edge    = Twirler(origin,  unit_z, unit_x, 2, 0)
  53     m.midline = Twirler(-unit_z, unit_z, unit_y, 1, 0)
  54     m._beziers = [ m._bezier(u) for u in np.linspace(0, 1, n_u) ]
  55   def _bezier(m,u):
  56     theta = u * tau
  57     cp = [None] * 4
  58     cp[0] =               m.edge   .point(theta)
  59     cp[1] = cp[0] + 0.5 * m.edge   .dirn (theta)
  60     cp[3] =               m.midline.point(theta*2)
  61     cp[2] = cp[3] + 0.5 * m.midline.dirn (theta*2)
  62     return BezierSegment(cp)
  63   def point(m, ix_u, t):
  64     return m._beziers[ix_u].point_at_t(t)