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 from moenp import *
   9 from moebez import *
  10
  11 import sys
  12
  13 class ParametricCircle:
  14   def __init__(pc, c, r0, r1):
  15     ''' circle centred on c
  16         with theta=0 point at c+r0
  17         and with theta=tau/4 point at c+r1 '''
  18     pc._c  = c
  19     pc._r0 = r0
  20     pc._r1 = r1
  21   def radius(pc, theta):
  22     return pc._r0 * cos(theta) + pc._r1 * sin(theta)
  23   def point(pc, theta):
  24     return pc._c + pc.radius(theta)
  25
  26 class Twirler(ParametricCircle):
  27   def __init__(tw, c, r0, r1, cycles, begin_zeta):
  28     ''' circle centred on c, etc.
  29         but with an orientation at each point, orthogonal to
  30           the circle
  31         the orientation circles round cycles times during the
  32           whole cycle (if cycles<1, to get a whole circling of
  33           the dirn around the circle, must pass some theta>tau)
  34         begin_zeta is the angle from outwards at theta==0
  35           positive meaning in the direction of r0 x r1 from r0
  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, extra_zeta=0):
  42     ''' unit vector for dirn at theta,
  43          but rotated extra_theta more around the circle
  44     '''
  45     zeta = tw._begin_zeta + theta * tw._cycles + extra_zeta
  46     r = tw.radius(theta)
  47     return cos(zeta) * r + sin(zeta) * tw._axis
  48
  49 class MoebiusHalf:
  50   def __init__(m, nu, nt):
  51     '''
  52      MoebiusHalf().edge is a Twirler for the edge,
  53       expecting theta = u * tau (see MoebiusHalf().point)
  54       with dirn pointing into the surface
  55     '''
  56     m.edge    = Twirler(origin,  unit_z, unit_x, -2, tau/2)
  57     m.midline = Twirler(-unit_z, unit_z, unit_y, -0.5, 0)
  58     m.nu      = nu
  59     m.nt      = nt
  60     m._thetas = [ u * tau for u in np.linspace(0, 1, nu+1) ]
  61     m._cp2b = BezierSegment([ (c,) for c in [0.33,0.33, 1.50]])
  62     m._dbeziers = [ m._dbezier(theta) for theta in m._thetas ]
  63   def _dbezier(m, theta, dconstructor=DiscreteBezier):
  64     cp = [None] * 4
  65     cp[0] =               m.edge   .point(theta)
  66     cp[3] =               m.midline.point(theta*2)
  67     ncp3 = np.linalg.norm(cp[3])
  68     cpt = ncp3 * ncp3 / 4
  69     cp2scale = m._cp2b.point_at_t(ncp3/2)
  70     cp1scale = 1.0 * cpt + 0.33 * (1-cpt)
  71     #print('u=%d ncp3=%f cp2scale=%f' % (u, ncp3, cp2scale),
  72     #      file=sys.stderr)
  73     cp[1] = cp[0] + cp1scale * m.edge   .dirn (theta)
  74     cp[2] = cp[3] + cp2scale * m.midline.dirn (theta*2)
  75     return dconstructor(cp, m.nt)
  76   def point(m, iu, it):
  77     '''
  78     0 <= iu <= nu     meaning 0 <= u <= 1
  79                       along the extent (well, along the edge)
  80                        0 and 1 are both the top half of the flat traverse
  81                        0.5        is the bottom half of the flat traverse
  82     0 <= it <= nt     across the half-traverse
  83                        0 is the edge, 1 is the midline
  84     '''
  85     return np.array(m._dbeziers[iu].point_at_it(it))
  86
  87   def details(m, iu, t):
  88     '''
  89     returns tuple of 4 vectors:
  90            - point location
  91            - normal (+ve is in the +ve y direction at iu=t=0) unit vector
  92            - along extent   (towrds +ve iu)                   unit vector
  93            - along traverse (towards +ve t)                   unit vector
  94     '''
  95     if iu == m.nu:
  96       return m.details(0, t)
  97     p = m.point(iu, t)
  98     vec_t = unit_v( m.point(iu, t + 0.01) - p )
  99     if t == 0:
 100       normal = m.edge.dirn(m._thetas[iu], extra_zeta=-tau/4)
 101       vec_u = np.cross(vec_t, normal)
 102     else:
 103       vec_u = unit_v( m.point(iu+1, t) - p )
 104       normal = np.cross(vec_u, vec_t)
 105     return p, normal, vec_u, vec_t
 106
 107   def point_offset(m, iu, t, offset):
 108     '''
 109     offset by offset perpendicular to the surface
 110      at the top (iu=t=0), +ve offset is in the +ve y direction
 111     '''
 112     p, normal, dummy, dummy = m.details(iu, t)
 113     return p + offset * normal
 114
 115 class Moebius():
 116   def __init__(m, nv, nw):
 117     '''
 118       0 <= v <= nw    along the extent,    v=0 is the flat traverse
 119       0 <= w <= nv    across the traverse  nw must be even
 120       the top is both   v=0, w=0  v=nv, w=nw
 121     '''
 122     assert(nw % 1 == 0)
 123     m.nv = nv
 124     m.nw = nw
 125     m.nt = nw/2
 126     m.h = MoebiusHalf(nu=nv*2, nt=m.nt)
 127
 128   def _vw2tiu_kw(m, v, w):
 129     if w <= m.nt:
 130       it = w
 131       iu = v
 132     else:
 133       it = m.nw - w
 134       iu = m.nv + v
 135     #print('v,w=%d,%d => it,iu=%d,%d' % (v,w,it,iu),
 136     #      file=sys.stderr)
 137     return { 'it': it, 'iu': iu }
 138
 139   def point(m, v, w):
 140     return m.h.point(**m._vw2tiu_kw(v,w))
 141
 142   def point_offset(m, v, w, offset):
 143     return m.h.point_offset(offset=
 144                             offset if w <= m.nt else -offset,
 145                             **m._vw2tiu_kw(v,w))
 146
 147   def details(m, v, w):
 148     '''
 149     returns tuple of 4 vectors:
 150            - point location
 151            - normal (+ve is in the +ve y direction at iu=t=0) unit vector
 152            - along extent   (towrds +ve v)                    unit vector
 153            - along traverse (towards +ve w)                   unit vector
 154     '''
 155     p, normal, vec_v, vec_w = m.h.details(**m._vw2tiu_kw(v,w))
 156     if w > m.nt:
 157       normal = -normal
 158       vec_w  = -vec_w
 159     return p, normal, vec_v, vec_w