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 (if cycles<1, to get a whole circling of
  38           the dirn around the circle, must pass some theta>tau)
  39         begin_zeta is the angle from outwards at theta==0
  40           positive meaning in the direction of r0 x r1 from r0
  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, extra_zeta=0):
  47     ''' unit vector for dirn at theta,
  48          but rotated extra_theta more around the circle
  49     '''
  50     zeta = tw._begin_zeta + theta * tw._cycles + extra_zeta
  51     r = tw.radius(theta)
  52     return cos(zeta) * r + sin(zeta) * tw._axis
  53
  54 class MoebiusHalf:
  55   def __init__(m, nu):
  56     '''
  57      MoebiusHalf().edge is a Twirler for the edge,
  58       expecting theta = u * tau (see MoebiusHalf().point)
  59       with dirn pointing into the surface
  60     '''
  61     m.edge    = Twirler(origin,  unit_z, unit_x, -2, tau/2)
  62     m.midline = Twirler(-unit_z, unit_z, unit_y, -0.5, 0)
  63     m.nu      = nu
  64     m._thetas = [ u * tau for u in np.linspace(0, 1, nu+1) ]
  65     m._beziers = [ m._bezier(theta) for theta in m._thetas ]
  66   def _bezier(m, theta):
  67     cp = [None] * 4
  68     cp[0] =               m.edge   .point(theta)
  69     cp[3] =               m.midline.point(theta*2)
  70     ncp3 = np.linalg.norm(cp[3])
  71     cpt = ncp3 * ncp3 / 4
  72     cp2scale = 1.50 * cpt + 0.33 * (1-cpt)
  73     cp1scale = 0.75 * cpt + 0.33 * (1-cpt)
  74     #print('u=%d ncp3=%f cp2scale=%f' % (u, ncp3, cp2scale),
  75     #      file=sys.stderr)
  76     cp[1] = cp[0] + cp1scale * m.edge   .dirn (theta)
  77     cp[2] = cp[3] + cp2scale * m.midline.dirn (theta*2)
  78     return BezierSegment(cp)
  79   def point(m, iu, t):
  80     '''
  81     0 <= iu <= nu     meaning 0 <= u <= 1
  82                       along the extent (well, along the edge)
  83                        0 and 1 are both the top half of the flat traverse
  84                        0.5        is the bottom half of the flat traverse
  85     0 <=  t   <= 1    across the half-traverse
  86                        0 is the edge, 1 is the midline
  87     '''
  88     return np.array(m._beziers[iu].point_at_t(t))
  89
  90   def point_offset(m, iu, t, offset):
  91     '''
  92     offset by offset perpendicular to the surface
  93      at the top (iu=t=0), in the y direction
  94     '''
  95     p = MoebiusHalf.point(m, iu, t)
  96     if t == 0:
  97       dirn = m.edge.dirn(m._thetas[iu], extra_zeta=-tau/4)
  98     elif iu == m.nu:
  99       return MoebiusHalf.point_offset(m, 0, t, offset)
 100     else:
 101       vec_t = MoebiusHalf.point(m, iu,   t + 0.01) - p
 102       vec_u = MoebiusHalf.point(m, iu+1, t)        - p
 103       dirn =  np.cross(vec_u, vec_t)
 104     return p + offset * dirn / np.linalg.norm(dirn)
 105
 106 class Moebius(MoebiusHalf):
 107   def __init__(m, nv, nw):
 108     '''
 109       0 <= v <= nw    along the extent,    v=0 is the flat traverse
 110       0 <= w <= nv    across the traverse  nw must be even
 111       the top is both   v=0, w=0  v=nv, w=nw
 112     '''
 113     assert(nw % 1 == 0)
 114     m.nv = nv
 115     m.nw = nw
 116     m.nt = nw/2
 117     m._t_vals = np.linspace(0, 1, m.nt+1)
 118     MoebiusHalf.__init__(m, nu=nv*2)
 119
 120   def _vw2tiu_kw(m, v, w):
 121     if w <= m.nt:
 122       it = w
 123       iu = v
 124     else:
 125       it = m.nw - w
 126       iu = m.nv + v
 127     #print('v,w=%d,%d => it,iu=%d,%d' % (v,w,it,iu),
 128     #      file=sys.stderr)
 129     return { 't': m._t_vals[it], 'iu': iu }
 130
 131   def point(m, v, w):
 132     return MoebiusHalf.point(m, **m._vw2tiu_kw(v,w))
 133
 134   def point_offset(m, v, w, offset):
 135     return MoebiusHalf.point_offset(m,
 136                                     offset=
 137                                     offset if w <= m.nt else -offset,
 138                                     **m._vw2tiu_kw(v,w))