X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?a=blobdiff_plain;f=moebius.py;h=e5a81acad9d0d458bbc823769b5cea55176e0801;hb=906edec3ef60cf3e0567c713b2a68688edb634d4;hp=4b4f3a0c6342a4340aac49b8c740a02861012b2d;hpb=e802fa0bb84fc9d4d35704f9758a40a3d9df439d;p=moebius3.git diff --git a/moebius.py b/moebius.py index 4b4f3a0..e5a81ac 100644 --- a/moebius.py +++ b/moebius.py @@ -1,14 +1,39 @@ +from __future__ import print_function + import numpy as np +from numpy import cos, sin + +from bezier import BezierSegment +from helixish import HelixishCurve +from moenp import * -tau = np.pi * 2 +import sys -origin = (0,0,0) -unit_x = (1,0,0) -unit_y = (0,1,0) -unit_z = (0,0,1) +class DoubleCubicBezier(): + def __init__(db, cp): + single = BezierSegment(cp) + midpoint = np.array(single.point_at_t(0.5)) + mid_dirn = single.point_at_t(0.5 + 0.001) - midpoint + mid_dirn /= np.linalg.norm(mid_dirn) + ocp_factor = 0.5 + mid_scale = ocp_factor * 0.5 * (np.linalg.norm(cp[1] - cp[0]) + + np.linalg.norm(cp[3] - cp[2])) + db.b0 = BezierSegment([ cp[0], + cp[1] * ocp_factor + cp[0] * (1-ocp_factor), + midpoint - mid_dirn * mid_scale, + midpoint ]) + db.b1 = BezierSegment([ midpoint, + midpoint + mid_dirn * mid_scale, + cp[2] * ocp_factor + cp[3] * (1-ocp_factor), + cp[3] ]) + def point_at_t(db, t): + if t < 0.5: + return db.b0.point_at_t(t*2) + else: + return db.b1.point_at_t(t*2 - 1) -def ParametricCircle: +class ParametricCircle: def __init__(pc, c, r0, r1): ''' circle centred on c with theta=0 point at c+r0 @@ -16,27 +41,144 @@ def ParametricCircle: pc._c = c pc._r0 = r0 pc._r1 = r1 - pc._axis = np.cross(r0, r1) + def radius(pc, theta): + return pc._r0 * cos(theta) + pc._r1 * sin(theta) def point(pc, theta): - return pc._c + pc._r0 * cos(theta) + pc._r1 * sin(theta) + return pc._c + pc.radius(theta) -def Twirler(ParametricCircle): - def __init__(tw, c, r0, r1, cycles, begin): +class Twirler(ParametricCircle): + def __init__(tw, c, r0, r1, cycles, begin_zeta): ''' circle centred on c, etc. but with an orientation at each point, orthogonal to - the circle -''' + the circle + the orientation circles round cycles times during the + whole cycle (if cycles<1, to get a whole circling of + the dirn around the circle, must pass some theta>tau) + begin_zeta is the angle from outwards at theta==0 + positive meaning in the direction of r0 x r1 from r0 + ''' ParametricCircle.__init__(tw, c, r0, r1) tw._cycles = cycles - def - -class Moebius: - def __init__(m, n_u): - m._edge = ParametricCircle(origin, unit_z, unit_x) - m._midline = ParametricCircle(-unit_z, unit_z, unit_y) - m._beziers = [ self._bezier(u) for u in np.linspace(0, 1, n_u) ] - def _bezier(u): - theta = u * tau - cp0 = m._edge.point(theta) - cp3 = m._midline.point(theta*2) - + tw._begin_zeta = begin_zeta + tw._axis = np.cross(r0, r1) + def dirn(tw, theta, extra_zeta=0): + ''' unit vector for dirn at theta, + but rotated extra_theta more around the circle + ''' + zeta = tw._begin_zeta + theta * tw._cycles + extra_zeta + r = tw.radius(theta) + return cos(zeta) * r + sin(zeta) * tw._axis + +class MoebiusHalf: + def __init__(m, nu): + ''' + MoebiusHalf().edge is a Twirler for the edge, + expecting theta = u * tau (see MoebiusHalf().point) + with dirn pointing into the surface + ''' + m.edge = Twirler(origin, unit_z, unit_x, -2, tau/2) + m.midline = Twirler(-unit_z, unit_z, unit_y, -0.5, 0) + m.nu = nu + m._thetas = [ u * tau for u in np.linspace(0, 1, nu+1) ] + m._cp2b = BezierSegment([ (c,) for c in [0.33,0.33, 1.50]]) + m._beziers = [ m._bezier(theta, HelixishCurve) for theta in m._thetas ] + #check = int(nu/3) + #m._beziers[check] = m._bezier(m._thetas[check], HelixishCurve) + def _bezier(m, theta, constructor=DoubleCubicBezier): + cp = [None] * 4 + cp[0] = m.edge .point(theta) + cp[3] = m.midline.point(theta*2) + ncp3 = np.linalg.norm(cp[3]) + cpt = ncp3 * ncp3 / 4 + cp2scale = m._cp2b.point_at_t(ncp3/2) + cp1scale = 1.0 * cpt + 0.33 * (1-cpt) + #print('u=%d ncp3=%f cp2scale=%f' % (u, ncp3, cp2scale), + # file=sys.stderr) + cp[1] = cp[0] + cp1scale * m.edge .dirn (theta) + cp[2] = cp[3] + cp2scale * m.midline.dirn (theta*2) + return constructor(cp) + def point(m, iu, t): + ''' + 0 <= iu <= nu meaning 0 <= u <= 1 + along the extent (well, along the edge) + 0 and 1 are both the top half of the flat traverse + 0.5 is the bottom half of the flat traverse + 0 <= t <= 1 across the half-traverse + 0 is the edge, 1 is the midline + ''' + return np.array(m._beziers[iu].point_at_t(t)) + + def details(m, iu, t): + ''' + returns tuple of 4 vectors: + - point location + - normal (+ve is in the +ve y direction at iu=t=0) unit vector + - along extent (towrds +ve iu) unit vector + - along traverse (towards +ve t) unit vector + ''' + if iu == m.nu: + return m.details(0, t) + p = m.point(iu, t) + vec_t = unit_v( m.point(iu, t + 0.01) - p ) + if t == 0: + normal = m.edge.dirn(m._thetas[iu], extra_zeta=-tau/4) + vec_u = np.cross(vec_t, normal) + else: + vec_u = unit_v( m.point(iu+1, t) - p ) + normal = np.cross(vec_u, vec_t) + return p, normal, vec_u, vec_t + + def point_offset(m, iu, t, offset): + ''' + offset by offset perpendicular to the surface + at the top (iu=t=0), +ve offset is in the +ve y direction + ''' + p, normal, dummy, dummy = m.details(iu, t) + return p + offset * normal + +class Moebius(): + def __init__(m, nv, nw): + ''' + 0 <= v <= nw along the extent, v=0 is the flat traverse + 0 <= w <= nv across the traverse nw must be even + the top is both v=0, w=0 v=nv, w=nw + ''' + assert(nw % 1 == 0) + m.nv = nv + m.nw = nw + m.nt = nw/2 + m._t_vals = np.linspace(0, 1, m.nt+1) + m.h = MoebiusHalf(nu=nv*2) + + def _vw2tiu_kw(m, v, w): + if w <= m.nt: + it = w + iu = v + else: + it = m.nw - w + iu = m.nv + v + #print('v,w=%d,%d => it,iu=%d,%d' % (v,w,it,iu), + # file=sys.stderr) + return { 't': m._t_vals[it], 'iu': iu } + + def point(m, v, w): + return m.h.point(**m._vw2tiu_kw(v,w)) + + def point_offset(m, v, w, offset): + return m.h.point_offset(offset= + offset if w <= m.nt else -offset, + **m._vw2tiu_kw(v,w)) + + def details(m, v, w): + ''' + returns tuple of 4 vectors: + - point location + - normal (+ve is in the +ve y direction at iu=t=0) unit vector + - along extent (towrds +ve v) unit vector + - along traverse (towards +ve w) unit vector + ''' + p, normal, vec_v, vec_w = m.h.details(**m._vw2tiu_kw(v,w)) + if w > m.nt: + normal = -normal + vec_w = -vec_w + return p, normal, vec_v, vec_w