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