2 from __future__ import print_function
5 from numpy import cos, sin
7 from bezier import BezierSegment
8 from helixish import HelixishCurve
13 class DoubleCubicBezier():
15 single = BezierSegment(cp)
16 midpoint = np.array(single.point_at_t(0.5))
17 mid_dirn = single.point_at_t(0.5 + 0.001) - midpoint
18 mid_dirn /= np.linalg.norm(mid_dirn)
20 mid_scale = ocp_factor * 0.5 * (np.linalg.norm(cp[1] - cp[0]) +
21 np.linalg.norm(cp[3] - cp[2]))
22 db.b0 = BezierSegment([ cp[0],
23 cp[1] * ocp_factor + cp[0] * (1-ocp_factor),
24 midpoint - mid_dirn * mid_scale,
26 db.b1 = BezierSegment([ midpoint,
27 midpoint + mid_dirn * mid_scale,
28 cp[2] * ocp_factor + cp[3] * (1-ocp_factor),
30 def point_at_t(db, t):
32 return db.b0.point_at_t(t*2)
34 return db.b1.point_at_t(t*2 - 1)
36 class ParametricCircle:
37 def __init__(pc, c, r0, r1):
38 ''' circle centred on c
39 with theta=0 point at c+r0
40 and with theta=tau/4 point at c+r1 '''
44 def radius(pc, theta):
45 return pc._r0 * cos(theta) + pc._r1 * sin(theta)
47 return pc._c + pc.radius(theta)
49 class Twirler(ParametricCircle):
50 def __init__(tw, c, r0, r1, cycles, begin_zeta):
51 ''' circle centred on c, etc.
52 but with an orientation at each point, orthogonal to
54 the orientation circles round cycles times during the
55 whole cycle (if cycles<1, to get a whole circling of
56 the dirn around the circle, must pass some theta>tau)
57 begin_zeta is the angle from outwards at theta==0
58 positive meaning in the direction of r0 x r1 from r0
60 ParametricCircle.__init__(tw, c, r0, r1)
62 tw._begin_zeta = begin_zeta
63 tw._axis = np.cross(r0, r1)
64 def dirn(tw, theta, extra_zeta=0):
65 ''' unit vector for dirn at theta,
66 but rotated extra_theta more around the circle
68 zeta = tw._begin_zeta + theta * tw._cycles + extra_zeta
70 return cos(zeta) * r + sin(zeta) * tw._axis
75 MoebiusHalf().edge is a Twirler for the edge,
76 expecting theta = u * tau (see MoebiusHalf().point)
77 with dirn pointing into the surface
79 m.edge = Twirler(origin, unit_z, unit_x, -2, tau/2)
80 m.midline = Twirler(-unit_z, unit_z, unit_y, -0.5, 0)
82 m._thetas = [ u * tau for u in np.linspace(0, 1, nu+1) ]
83 m._cp2b = BezierSegment([ (c,) for c in [0.33,0.33, 1.50]])
84 m._beziers = [ m._bezier(theta) for theta in m._thetas ]
86 m._bezier(m._thetas[check], HelixishCurve)
87 m._bezier(m._thetas[check+1], HelixishCurve)
88 m._bezier(m._thetas[check+2], HelixishCurve)
89 def _bezier(m, theta, constructor=DoubleCubicBezier):
91 cp[0] = m.edge .point(theta)
92 cp[3] = m.midline.point(theta*2)
93 ncp3 = np.linalg.norm(cp[3])
95 cp2scale = m._cp2b.point_at_t(ncp3/2)
96 cp1scale = 1.0 * cpt + 0.33 * (1-cpt)
97 #print('u=%d ncp3=%f cp2scale=%f' % (u, ncp3, cp2scale),
99 cp[1] = cp[0] + cp1scale * m.edge .dirn (theta)
100 cp[2] = cp[3] + cp2scale * m.midline.dirn (theta*2)
101 return constructor(cp)
104 0 <= iu <= nu meaning 0 <= u <= 1
105 along the extent (well, along the edge)
106 0 and 1 are both the top half of the flat traverse
107 0.5 is the bottom half of the flat traverse
108 0 <= t <= 1 across the half-traverse
109 0 is the edge, 1 is the midline
111 return np.array(m._beziers[iu].point_at_t(t))
113 def details(m, iu, t):
115 returns tuple of 4 vectors:
117 - normal (+ve is in the +ve y direction at iu=t=0) unit vector
118 - along extent (towrds +ve iu) unit vector
119 - along traverse (towards +ve t) unit vector
122 return m.details(0, t)
124 vec_t = unit_v( m.point(iu, t + 0.01) - p )
126 normal = m.edge.dirn(m._thetas[iu], extra_zeta=-tau/4)
127 vec_u = np.cross(vec_t, normal)
129 vec_u = unit_v( m.point(iu+1, t) - p )
130 normal = np.cross(vec_u, vec_t)
131 return p, normal, vec_u, vec_t
133 def point_offset(m, iu, t, offset):
135 offset by offset perpendicular to the surface
136 at the top (iu=t=0), +ve offset is in the +ve y direction
138 p, normal, dummy, dummy = m.details(iu, t)
139 return p + offset * normal
142 def __init__(m, nv, nw):
144 0 <= v <= nw along the extent, v=0 is the flat traverse
145 0 <= w <= nv across the traverse nw must be even
146 the top is both v=0, w=0 v=nv, w=nw
152 m._t_vals = np.linspace(0, 1, m.nt+1)
153 m.h = MoebiusHalf(nu=nv*2)
155 def _vw2tiu_kw(m, v, w):
162 #print('v,w=%d,%d => it,iu=%d,%d' % (v,w,it,iu),
164 return { 't': m._t_vals[it], 'iu': iu }
167 return m.h.point(**m._vw2tiu_kw(v,w))
169 def point_offset(m, v, w, offset):
170 return m.h.point_offset(offset=
171 offset if w <= m.nt else -offset,
174 def details(m, v, w):
176 returns tuple of 4 vectors:
178 - normal (+ve is in the +ve y direction at iu=t=0) unit vector
179 - along extent (towrds +ve v) unit vector
180 - along traverse (towards +ve w) unit vector
182 p, normal, vec_v, vec_w = m.h.details(**m._vw2tiu_kw(v,w))
186 return p, normal, vec_v, vec_w