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