def solve_polynomial_join(left, right, reverse=0):
    "\n  Solves for a parametric cubic polynomial curve joining the right side of left\n  to the left side of right.  The curve matches slope and position at it's\n  boundaries and is parameterized from 0 to 1; 0 being the left boundary and 1\n  being the right.\n\n  method:       parametric cubic matching position and slope of endpoints.\n                This ends up being cheap to compute, since the matrix is\n                known (interval of parameter is always 0 to 1) and so the \n                inverse can be precomputed. \n     minv is inverse of m, where:\n          m = array( [  [   a**3,   a**2, a, 1 ],\n                        [   b**3,   b**2, b, 1 ],  \n                        [ 3*a**2, 2*a   , 1, 0 ],\n                        [ 3*b**2, 2*b   , 1, 0 ] ] )\n     is the matrix for the linear system:\n          m * coeff = v,\n     with v = [ x(0) x(1) dx/dt(0) dx/dt(1) ].\n     Here a = 0 and b = 1 so m and it's inverse is always the same.\n  \n  "
    minv = matrix([[2.0, (- 2.0), 1.0, 1.0], [(- 3.0), 3.0, (- 2.0), (- 1.0)], [0.0, 0.0, 1.0, 0.0], [1.0, 0.0, 0.0, 0.0]])
    lr = len(right)
    ll = len(left)
    nl = (ll / 4)
    nr = (lr / 4)
    slope = (lambda v: (v[0] - v[(- 1)]))
    length = (lambda x, y: hypot(diff(x), diff(y)).sum())
    if ((nr < 2) and (nl < 2)):
        lnorm = length(left.x, left.y)
        rnorm = length(right.x, right.y)
        dly = (diff(left.y).mean() / lnorm)
        dlx = (diff(left.x).mean() / lnorm)
        dry = (diff(right.y).mean() / rnorm)
        drx = (diff(right.x).mean() / rnorm)
        nl = 0
        nr = (lr - 1)
    elif (nr < 2):
        lnorm = length(left.x[:nl], left.y[:nl])
        rnorm = length(right.x, right.y)
        dly = ((- slope(left.y[(- nl):])) / lnorm)
        dlx = ((- slope(left.x[(- nl):])) / lnorm)
        dry = (diff(right.y).mean() / rnorm)
        drx = (diff(right.x).mean() / rnorm)
        nr = (lr - 1)
    elif (nl < 2):
        rnorm = length(right.x[:nr], right.y[:nr])
        lnorm = length(left.x, left.y)
        dry = ((- slope(right.y[:nr])) / rnorm)
        drx = ((- slope(right.x[:nr])) / rnorm)
        dly = (diff(left.y).mean() / lnorm)
        dlx = (diff(left.x).mean() / lnorm)
        nl = 0
    else:
        rnorm = length(right.x[:nr], right.y[:nr])
        lnorm = length(left.x[(- nl):], left.y[(- nl):])
        dry = ((- slope(right.y[:nr])) / rnorm)
        drx = ((- slope(right.x[:nr])) / rnorm)
        dly = ((- slope(left.y[(- nl):])) / lnorm)
        dlx = ((- slope(left.x[(- nl):])) / lnorm)
    rnorm = hypot((left.x[0] - right.x[0]), (left.y[0] - right.y[0]))
    lnorm = hypot((left.x[(- 1)] - right.x[0]), (left.y[(- 1)] - right.y[0]))
    if (not isfinite(dlx)):
        dlx = ((left.x[0] - right.x[0]) / lnorm)
    if (not isfinite(dly)):
        dly = ((left.y[0] - right.y[0]) / lnorm)
    if (not isfinite(drx)):
        drx = ((left.x[(- 1)] - right.x[0]) / rnorm)
    if (not isfinite(dry)):
        dry = ((left.y[(- 1)] - right.y[0]) / rnorm)
    if reverse:
        dlx = (- dlx)
        dly = (- dly)
        drx = (- drx)
        dry = (- dry)
    ry = right.y[0]
    ly = left.y[(- 1)]
    rx = right.x[0]
    lx = left.x[(- 1)]
    L = hypot((rx - lx), (ry - ly))
    print(('L:%g' % L))
    yv = matrix([[ly], [ry], [(dly * L)], [(dry * L)]])
    xv = matrix([[lx], [rx], [(dlx * L)], [(drx * L)]])
    cx = (minv * xv)
    cy = (minv * yv)
    if (not (isfinite(cx).any() and isfinite(cy).any())):
        pdb.set_trace()
    return [array(t).squeeze() for t in (cx, cy)]
Solves for a parametric cubic polynomial curve joining the right side of left
to the left side of right.  The curve matches slope and position at it's
boundaries and is parameterized from 0 to 1; 0 being the left boundary and 1
being the right.

method:       parametric cubic matching position and slope of endpoints.
              This ends up being cheap to compute, since the matrix is
              known (interval of parameter is always 0 to 1) and so the 
              inverse can be precomputed. 
   minv is inverse of m, where:
        m = array( [  [   a**3,   a**2, a, 1 ],
                      [   b**3,   b**2, b, 1 ],  
                      [ 3*a**2, 2*a   , 1, 0 ],
                      [ 3*b**2, 2*b   , 1, 0 ] ] )
   is the matrix for the linear system:
        m * coeff = v,
   with v = [ x(0) x(1) dx/dt(0) dx/dt(1) ].
   Here a = 0 and b = 1 so m and it's inverse is always the same.
def solve_polynomial_join(left, right, reverse=0):
    "\n  Solves for a parametric cubic polynomial curve joining the right side of left\n  to the left side of right.  The curve matches slope and position at it's\n  boundaries and is parameterized from 0 to 1; 0 being the left boundary and 1\n  being the right.\n\n  method:       parametric cubic matching position and slope of endpoints.\n                This ends up being cheap to compute, since the matrix is\n                known (interval of parameter is always 0 to 1) and so the \n                inverse can be precomputed. \n     minv is inverse of m, where:\n          m = array( [  [   a**3,   a**2, a, 1 ],\n                        [   b**3,   b**2, b, 1 ],  \n                        [ 3*a**2, 2*a   , 1, 0 ],\n                        [ 3*b**2, 2*b   , 1, 0 ] ] )\n     is the matrix for the linear system:\n          m * coeff = v,\n     with v = [ x(0) x(1) dx/dt(0) dx/dt(1) ].\n     Here a = 0 and b = 1 so m and it's inverse is always the same.\n  \n  "
    minv = matrix([[2.0, (- 2.0), 1.0, 1.0], [(- 3.0), 3.0, (- 2.0), (- 1.0)], [0.0, 0.0, 1.0, 0.0], [1.0, 0.0, 0.0, 0.0]])
    lr = len(right)
    ll = len(left)
    nl = (ll / 4)
    nr = (lr / 4)
    slope = (lambda v: (v[0] - v[(- 1)]))
    length = (lambda x, y: hypot(diff(x), diff(y)).sum())
    if ((nr < 2) and (nl < 2)):
        lnorm = length(left.x, left.y)
        rnorm = length(right.x, right.y)
        dly = (diff(left.y).mean() / lnorm)
        dlx = (diff(left.x).mean() / lnorm)
        dry = (diff(right.y).mean() / rnorm)
        drx = (diff(right.x).mean() / rnorm)
        nl = 0
        nr = (lr - 1)
    elif (nr < 2):
        lnorm = length(left.x[:nl], left.y[:nl])
        rnorm = length(right.x, right.y)
        dly = ((- slope(left.y[(- nl):])) / lnorm)
        dlx = ((- slope(left.x[(- nl):])) / lnorm)
        dry = (diff(right.y).mean() / rnorm)
        drx = (diff(right.x).mean() / rnorm)
        nr = (lr - 1)
    elif (nl < 2):
        rnorm = length(right.x[:nr], right.y[:nr])
        lnorm = length(left.x, left.y)
        dry = ((- slope(right.y[:nr])) / rnorm)
        drx = ((- slope(right.x[:nr])) / rnorm)
        dly = (diff(left.y).mean() / lnorm)
        dlx = (diff(left.x).mean() / lnorm)
        nl = 0
    else:
        rnorm = length(right.x[:nr], right.y[:nr])
        lnorm = length(left.x[(- nl):], left.y[(- nl):])
        dry = ((- slope(right.y[:nr])) / rnorm)
        drx = ((- slope(right.x[:nr])) / rnorm)
        dly = ((- slope(left.y[(- nl):])) / lnorm)
        dlx = ((- slope(left.x[(- nl):])) / lnorm)
    rnorm = hypot((left.x[0] - right.x[0]), (left.y[0] - right.y[0]))
    lnorm = hypot((left.x[(- 1)] - right.x[0]), (left.y[(- 1)] - right.y[0]))
    if (not isfinite(dlx)):
        dlx = ((left.x[0] - right.x[0]) / lnorm)
    if (not isfinite(dly)):
        dly = ((left.y[0] - right.y[0]) / lnorm)
    if (not isfinite(drx)):
        drx = ((left.x[(- 1)] - right.x[0]) / rnorm)
    if (not isfinite(dry)):
        dry = ((left.y[(- 1)] - right.y[0]) / rnorm)
    if reverse:
        dlx = (- dlx)
        dly = (- dly)
        drx = (- drx)
        dry = (- dry)
    ry = right.y[0]
    ly = left.y[(- 1)]
    rx = right.x[0]
    lx = left.x[(- 1)]
    L = hypot((rx - lx), (ry - ly))
    print(('L:%g' % L))
    yv = matrix([[ly], [ry], [(dly * L)], [(dry * L)]])
    xv = matrix([[lx], [rx], [(dlx * L)], [(drx * L)]])
    cx = (minv * xv)
    cy = (minv * yv)
    if (not (isfinite(cx).any() and isfinite(cy).any())):
        pdb.set_trace()
    return [array(t).squeeze() for t in (cx, cy)]