test/experimental/newbox.py

   1 import sys; sys.path.insert(0, "../..")
   2
   3 import math
   4 from pyx import *
   5 from pyx import normpath
   6
   7 # make math.sqrt raising an exception for negativ values
   8 try:
   9     math.sqrt(-1)
  10 except:
  11     pass
  12 else:
  13     math._sqrt = math.sqrt
  14     def mysqrt(x):
  15         if x < 0:
  16             raise ValueError("sqrt of negativ number")
  17         return math._sqrt(x)
  18     math.sqrt = mysqrt
  19
  20
  21 c = canvas.canvas()
  22
  23 def connect(normpathitem1, normpathitem2, round):
  24     "returns a list of normsubpathitems connecting normpathitem1 and normpathitem2 either rounded or not"
  25     # this is for corners, i.e. when there are jumps in the tangent vector
  26     t1 = normpathitem1.trafo([1])[0]
  27     t2 = normpathitem2.trafo([0])[0]
  28     xs, ys = t1.apply_pt(0, 0)
  29     xe, ye = t2.apply_pt(0, 0)
  30     if (xs-xe)*(xs-xe) + (xs-xe)*(xs-xe) < 1e-5**2:
  31         return []
  32     if round:
  33         # calculate crossing point of normal
  34         try:
  35             param = ((t2.matrix[0][1]*(t1.vector[1]-t2.vector[1]) -
  36                       t2.matrix[1][1]*(t1.vector[0]-t2.vector[0])) /
  37                      (t2.matrix[1][1]*t1.matrix[0][1] -
  38                       t2.matrix[0][1]*t1.matrix[1][1]))
  39         except ArithmeticError:
  40             return [path.normline_pt(xs, ys, xe, ye)]
  41         else:
  42             # since t1 is a translation + rotation, param is the radius
  43             # hence, param should be equal to enlargeby_pt
  44             # print param
  45             x, y = t1.apply_pt(0, param)
  46             # c.fill(path.circle_pt(x, y, 1))
  47             angle1 = math.atan2(ys-y, xs-x)
  48             angle2 = math.atan2(ye-y, xe-x)
  49             while angle2 < angle1:
  50                 angle1 -= 2*math.pi
  51             return [path._arctobcurve(x, y, param, angle1, angle2)]
  52     else:
  53         # calculate crossing point of tangents
  54         try:
  55             param = ((t2.matrix[0][0]*(t1.vector[1]-t2.vector[1]) -
  56                       t2.matrix[1][0]*(t1.vector[0]-t2.vector[0])) /
  57                      (t2.matrix[1][0]*t1.matrix[0][0] -
  58                       t2.matrix[0][0]*t1.matrix[1][0]))
  59         except ArithmeticError:
  60             return [path.normline_pt(xs, ys, xe, ye)]
  61         else:
  62             x, y = t1.apply_pt(param, 0)
  63             return [path.normline_pt(xs, ys, x, y), path.normline_pt(x, y, xe, ye)]
  64
  65
  66 def enlarged(anormpath, enlargeby_pt, round):
  67     newnormsubpaths = []
  68     for normsubpath in anormpath.normsubpaths:
  69         splitnormsubpathitems = normsubpath.normsubpathitems[:] # do splitting on a copy
  70         i = 0
  71         while i < len(splitnormsubpathitems):
  72             if isinstance(splitnormsubpathitems[i], normpath.normcurve_pt) and splitnormsubpathitems[i].arclen_pt(normsubpath.epsilon) > 100:
  73                 splitnormsubpathitems[i:i+1] = splitnormsubpathitems[i]._midpointsplit(normsubpath.epsilon)
  74             else:
  75                 i += 1
  76         newnormsubpathitems = []
  77         for normsubpathitem in splitnormsubpathitems:
  78             # get old and new start and end points
  79             ts, te = normsubpathitem.trafo([0, 1])
  80             xs, ys = ts.apply_pt(0, 0)
  81             nxs, nys = ts.apply_pt(0, -enlargeby_pt)
  82             xe, ye = te.apply_pt(0, 0)
  83             nxe, nye = te.apply_pt(0, -enlargeby_pt)
  84
  85             if isinstance(normsubpathitem, normpath.normcurve_pt):
  86                 # We should do not alter the sign. Could we do any better here?
  87                 try:
  88                     cs = 1/(normsubpathitem.curveradius_pt([0])[0] + enlargeby_pt)
  89                 except ArithmeticError:
  90                     cs = 0
  91                 try:
  92                     ce = 1/(normsubpathitem.curveradius_pt([1])[0] + enlargeby_pt)
  93                 except ArithmeticError:
  94                     ce = 0
  95
  96                 # this should be a function (in path?), not a method in a class
  97                 # the parameter convention is quite different from other places ...
  98                 bezierparams = deformer.normcurve_from_endgeometry_pt((nxs, nys), (nxe, nye),
  99                                                                       (ts.matrix[0][0], ts.matrix[1][0]),
 100                                                                       (te.matrix[0][0], te.matrix[1][0]),
 101                                                                       cs, ce)
 102                 c.fill(path.circle_pt(bezierparams.x0_pt, bezierparams.y0_pt, 1), [color.rgb.blue])
 103                 c.fill(path.circle_pt(bezierparams.x3_pt, bezierparams.y3_pt, 1), [color.rgb.blue])
 104                 newnormpathitems = normpath.normcurve_pt(bezierparams.x0_pt, bezierparams.y0_pt,
 105                                                bezierparams.x1_pt, bezierparams.y1_pt,
 106                                                bezierparams.x2_pt, bezierparams.y2_pt,
 107                                                bezierparams.x3_pt, bezierparams.y3_pt)
 108                 showtangent(newnormpathitems) # line alignment of bezier curves
 109                 showcircle(newnormpathitems)  # circle alignment of bezier curves
 110             else:
 111                 # line
 112                 newnormpathitems = path.normline_pt(nxs, nys, nxe, nye)
 113
 114             if len(newnormsubpathitems):
 115                 newnormsubpathitems.extend(connect(newnormsubpathitems[-1], newnormpathitems, round=round))
 116             newnormsubpathitems.append(newnormpathitems)
 117         if normsubpath.closed:
 118             newnormsubpathitems.extend(connect(newnormsubpathitems[-1], newnormsubpathitems[0], round=round))
 119         newnormsubpaths.append(path.normsubpath(newnormsubpathitems, normsubpath.closed))
 120     return normpath.normpath(newnormsubpaths)
 121
 122
 123 def showtangent(normcurve_pt):
 124     dx = 2
 125     dy = 1
 126
 127     cx =                                       3*normcurve_pt.x1_pt - 3*normcurve_pt.x0_pt
 128     bx =                   3*normcurve_pt.x2_pt - 6*normcurve_pt.x1_pt + 3*normcurve_pt.x0_pt
 129     ax = normcurve_pt.x3_pt - 3*normcurve_pt.x2_pt + 3*normcurve_pt.x1_pt -   normcurve_pt.x0_pt
 130
 131     cy =                                       3*normcurve_pt.y1_pt - 3*normcurve_pt.y0_pt
 132     by =                   3*normcurve_pt.y2_pt - 6*normcurve_pt.y1_pt + 3*normcurve_pt.y0_pt
 133     ay = normcurve_pt.y3_pt - 3*normcurve_pt.y2_pt + 3*normcurve_pt.y1_pt -   normcurve_pt.y0_pt
 134
 135     try:
 136         x = - (dy*bx-dx*by)/(3.0*dy*ax-3.0*dx*ay)
 137         s = math.sqrt(x*x-(dy*cx-dx*cy)/(3.0*dy*ax-3.0*dx*ay))
 138     except (ArithmeticError, ValueError):
 139         return
 140     t1 = x-s
 141     t2 = x+s
 142     ts = []
 143     if 0 < t1 < 1:
 144         ts.append(t1)
 145     if 0 < t2 < 1:
 146         ts.append(t2)
 147     for t in ts:
 148         trafo = normcurve_pt.trafo([t])[0]
 149         x, y = trafo.apply_pt(0, 0)
 150         c.fill(path.circle_pt(x, y, 1), [color.rgb.red])
 151         c.stroke(path.line_pt(*list(trafo.apply_pt(-100, 0))+list(trafo.apply_pt(100, 0))), [color.rgb.red])
 152
 153 def showcircle(normcurve_pt):
 154     gx = 350
 155     gy = 200
 156     hx = -1
 157     hy = 2
 158     cos = hx/math.sqrt(hx*hx+hy*hy)
 159     sin = hy/math.sqrt(hx*hx+hy*hy)
 160
 161     r = 16
 162
 163     dx =                                                             normcurve_pt.x0_pt
 164     cx =                                       3*normcurve_pt.x1_pt - 3*normcurve_pt.x0_pt
 165     bx =                   3*normcurve_pt.x2_pt - 6*normcurve_pt.x1_pt + 3*normcurve_pt.x0_pt
 166     ax = normcurve_pt.x3_pt - 3*normcurve_pt.x2_pt + 3*normcurve_pt.x1_pt -   normcurve_pt.x0_pt
 167
 168     dy =                                                             normcurve_pt.y0_pt
 169     cy =                                       3*normcurve_pt.y1_pt - 3*normcurve_pt.y0_pt
 170     by =                   3*normcurve_pt.y2_pt - 6*normcurve_pt.y1_pt + 3*normcurve_pt.y0_pt
 171     ay = normcurve_pt.y3_pt - 3*normcurve_pt.y2_pt + 3*normcurve_pt.y1_pt -   normcurve_pt.y0_pt
 172
 173     def x(t, gx=gx):
 174         return ax*t*t*t+bx*t*t+cx*t+normcurve_pt.x0_pt
 175     def y(t, gy=gy):
 176         return ay*t*t*t+by*t*t+cy*t+normcurve_pt.y0_pt
 177     def xdot(t):
 178         return 3*ax*t*t+2*bx*t+cx
 179     def ydot(t):
 180         return 3*ay*t*t+2*by*t+cy
 181
 182     # we need to find roots of the following polynom:
 183     # (xdot(t)**2+ydot(t)**2)*((x(t)-gx)*hy-(y(t)-gy)*hx)**2-r*r*(xdot(t)*hx+ydot(t)*hy)**2 == 0
 184
 185     # the polynom is of order 10. its coefficients are:
 186
 187     coeffs = [9*(ax**2+ay**2)*(hy*ax-ay*hx)**2, # * t**10,
 188
 189               -6*(hy*ax-ay*hx)*(-5*ax**2*bx*hy+3*ax**2*hx*by-2*ax*by*ay*hy+
 190               2*ax*hx*bx*ay-3*ay**2*bx*hy+5*hx*ay**2*by), # * t**9, etc.
 191
 192               -30*hy*hx*ay*ax**2*cx-32*hy*hx*ay*bx**2*ax-42*hy*hx*bx*ax**2*by-18*hy*hx*ax**3*cy+
 193               24*hy**2*ax**3*cx+37*hy**2*bx**2*ax**2+24*hx**2*ay**3*cy+37*hx**2*by**2*ay**2-
 194               30*hy*hx*ay**2*ax*cy-32*hy*hx*ay*ax*by**2-42*hy*hx*by*ay**2*bx-18*hy*hx*ay**3*cx+
 195               6*hy**2*ax**2*cy*ay+4*hy**2*ax**2*by**2+24*hy**2*bx*ax*by*ay+18*hy**2*ay**2*cx*ax+
 196               9*hy**2*ay**2*bx**2+6*hx**2*ay**2*cx*ax+4*hx**2*ay**2*bx**2+24*hx**2*bx*ax*by*ay+
 197               18*hx**2*ax**2*cy*ay+9*hx**2*ax**2*by**2,
 198
 199               20*hx**2*by**3*ay+8*hx**2*by*ay*bx**2+4*hx**2*cx*bx*ay**2-18*hx**2*ax**2*ay*gy+
 200               18*hx**2*ax**2*ay*dy+20*hy**2*bx**3*ax+18*hy**2*ax**3*dx+12*hx**2*by*ay*cx*ax+18*
 201               hx**2*cy*by*ax**2+12*hx**2*bx*ax*by**2+8*hy**2*bx*ax*by**2+4*hy**2*cy*by*ax**2-
 202               18*hy**2*ay**2*ax*gx+18*hy**2*ay**2*ax*dx+18*hy**2*cx*bx*ay**2+12*hy**2*by*ay*bx**2
 203               +58*hy**2*cx*bx*ax**2-8*hy*hx*ay*bx**3+18*hy*hx*ax**3*gy-18*hy*hx*ax**3*dy-
 204               18*hy*hx*ay**3*dx-8*hy*hx*ax*by**3+58*hx**2*cy*by*ay**2+24*hx**2*bx*ax*cy*ay+
 205               12*hy**2*bx*ax*cy*ay+24*hy**2*by*ay*cx*ax-18*hy**2*ax**3*gx-44*hy*hx*ay*cx*bx*ax+
 206               18*hy*hx*ay*ax**2*gx-18*hy*hx*ay*ax**2*dx-30*hy*hx*by*ax**2*cx-32*hy*hx*by*bx**2*ax-
 207               42*hy*hx*bx*ax**2*cy+18*hy*hx*ay**3*gx-44*hy*hx*ay*cy*by*ax-30*hy*hx*bx*cy*ay**2-
 208               32*hy*hx*ay*bx*by**2-42*hy*hx*by*ay**2*cx+18*hy*hx*ay**2*ax*gy-
 209               18*hy*hx*ay**2*ax*dy-18*hx**2*ay**3*gy+18*hx**2*ay**3*dy,
 210
 211               hx**2*cx**2*ay**2+46*hy**2*cx*bx**2*ax-42*hy**2*bx*ax**2*gx+42*hy**2*bx*ax**2*dx+
 212               46*hx**2*cy*by**2*ay-42*hx**2*by*ay**2*gy+42*hx**2*by*ay**2*dy-8*hy*hx*bx*by**3+
 213               6*hy**2*cy*ay*bx**2+8*hy**2*by**2*cx*ax-18*hy**2*ay**2*bx*gx+18*hy**2*ay**2*bx*dx+
 214               6*hx**2*by**2*cx*ax+8*hx**2*cy*ay*bx**2-18*hx**2*ax**2*by*gy+18*hx**2*ax**2*by*dy-
 215               8*hy*hx*by*bx**3+4*hy**2*bx**4+22*hy**2*cx**2*ax**2+22*hx**2*cy**2*ay**2+hy**2*cy**2*ax**2+
 216               4*hy**2*by**2*bx**2+9*hy**2*cx**2*ay**2+4*hx**2*by**2*bx**2+9*hx**2*cy**2*ax**2+
 217               4*hx**2*by**4-14*hy*hx*ay*cy**2*ax-44*hy*hx*ay*cy*by*bx-30*hy*hx*cx*cy*ay**2-
 218               32*hy*hx*ay*cx*by**2+42*hy*hx*by*ay**2*gx-42*hy*hx*by*ay**2*dx-16*hy*hx*cy*by**2*ax+
 219               24*hy*hx*by*ay*ax*gy-24*hy*hx*by*ay*ax*dy+18*hy*hx*ay**2*bx*gy-18*hy*hx*ay**2*bx*dy+
 220               8*hy**2*cy*by*bx*ax+12*hy**2*cy*ay*cx*ax-24*hy**2*by*ay*ax*gx+24*hy**2*by*ay*ax*dx+
 221               24*hy**2*cx*bx*by*ay+8*hx**2*cx*bx*by*ay+12*hx**2*cy*ay*cx*ax-24*hx**2*bx*ax*ay*gy+
 222               24*hx**2*bx*ax*ay*dy+24*hx**2*cy*by*bx*ax-14*hy*hx*ay*cx**2*ax-16*hy*hx*ay*cx*bx**2+
 223               24*hy*hx*ay*bx*ax*gx-24*hy*hx*ay*bx*ax*dx-44*hy*hx*by*cx*bx*ax+18*hy*hx*by*ax**2*gx-
 224               18*hy*hx*by*ax**2*dx-30*hy*hx*cy*ax**2*cx-32*hy*hx*cy*bx**2*ax+42*hy*hx*bx*ax**2*gy-
 225               42*hy*hx*bx*ax**2*dy,
 226
 227               12*hx**2*cy*by**3+12*hy**2*cx*bx**3+34*hx**2*cy**2*by*ay-30*hx**2*cy*ay**2*gy+
 228               30*hx**2*cy*ay**2*dy-32*hx**2*by**2*ay*gy+32*hx**2*by**2*ay*dy+2*hx**2*cx**2*by*ay+
 229               4*hx**2*cx*bx*by**2-8*hx**2*bx**2*ay*gy+8*hx**2*bx**2*ay*dy+8*hx**2*bx**2*cy*by+
 230               12*hx**2*cy**2*bx*ax-18*hx**2*ax**2*cy*gy+8*hx**2*cx*bx*cy*ay-12*hx**2*cx*ax*ay*gy+
 231               12*hx**2*cx*ax*ay*dy+12*hx**2*cx*ax*cy*by-24*hx**2*bx*ax*by*gy+24*hx**2*bx*ax*by*dy+
 232               18*hx**2*ax**2*cy*dy-8*hy*hx*cy*bx**3-8*hy*hx*cx*by**3+34*hy**2*cx**2*bx*ax-
 233               30*hy**2*cx*ax**2*gx+30*hy**2*cx*ax**2*dx-32*hy**2*bx**2*ax*gx+32*hy**2*bx**2*ax*dx+
 234               2*hy**2*cy**2*bx*ax+4*hy**2*bx**2*cy*by-8*hy**2*by**2*ax*gx+8*hy**2*by**2*ax*dx+
 235               8*hy**2*cx*bx*by**2+12*hy**2*cx**2*by*ay-18*hy**2*ay**2*cx*gx+18*hy**2*ay**2*cx*dx-
 236               10*hy*hx*ay*cx**2*bx+12*hy*hx*ay*cx*ax*gx-12*hy*hx*ay*cx*ax*dx+
 237               8*hy*hx*ay*bx**2*gx-8*hy*hx*ay*bx**2*dx-14*hy*hx*by*cx**2*ax-16*hy*hx*by*cx*bx**2+
 238               24*hy*hx*by*bx*ax*gx-24*hy*hx*by*bx*ax*dx-44*hy*hx*cy*cx*bx*ax+18*hy*hx*cy*ax**2*gx-
 239               18*hy*hx*cy*ax**2*dx+30*hy*hx*ax**2*cx*gy-30*hy*hx*ax**2*cx*dy+32*hy*hx*bx**2*ax*gy-
 240               32*hy*hx*bx**2*ax*dy-32*hy*hx*ay*by**2*dx-16*hy*hx*cy*by**2*bx+8*hy*hx*ax*by**2*gy+
 241               32*hy*hx*ay*by**2*gx-44*hy*hx*ay*cy*by*cx+12*hy*hx*ax*cy*ay*gy-12*hy*hx*ax*cy*ay*dy+
 242               24*hy*hx*by*ay*bx*gy-24*hy*hx*by*ay*bx*dy+30*hy*hx*cy*ay**2*gx-30*hy*hx*cy*ay**2*dx-
 243               14*hy*hx*ay*cy**2*bx-10*hy*hx*by*cy**2*ax-8*hy*hx*ax*by**2*dy+18*hy*hx*ay**2*cx*gy-
 244               18*hy*hx*ay**2*cx*dy+8*hy**2*cx*ax*cy*by-12*hy**2*cy*ay*ax*gx+12*hy**2*cy*ay*ax*dx+
 245               12*hy**2*cx*bx*cy*ay-24*hy**2*by*ay*bx*gx+24*hy**2*by*ay*bx*dx,
 246
 247               hy**2*bx**2*cy**2+13*hy**2*cx**2*bx**2+9*hx**2*ax**2*gy**2-9*r**2*ax**2*hx**2-
 248               9*r**2*ay**2*hy**2+9*hx**2*ax**2*dy**2+hx**2*cx**2*by**2+4*hx**2*bx**2*cy**2+8*hy**2*cx**3*ax+
 249               8*hx**2*cy**3*ay+13*hx**2*cy**2*by**2-8*hx**2*by**3*gy+8*hx**2*by**3*dy+9*hx**2*ay**2*gy**2+
 250               9*hx**2*ay**2*dy**2+4*hy**2*cx**2*by**2+9*hy**2*ay**2*gx**2+9*hy**2*ay**2*dx**2-8*hy**2*bx**3*gx+
 251               8*hy**2*bx**3*dx+9*hy**2*ax**2*gx**2+9*hy**2*ax**2*dx**2-18*hx**2*ay**2*gy*dy+2*hx**2*cx**2*cy*ay+
 252               6*hx**2*cx*ax*cy**2-8*hx**2*bx**2*by*gy+8*hx**2*bx**2*by*dy-18*hx**2*ax**2*gy*dy+8*hy*hx*by**3*gx-
 253               8*hy*hx*by**3*dx-2*hy*hx*cy**3*ax+2*hy**2*cx*ax*cy**2+6*hy**2*cx**2*cy*ay-8*hy**2*by**2*bx*gx+
 254               8*hy**2*by**2*bx*dx-18*hy**2*ay**2*gx*dx-18*hy**2*ax**2*gx*dx-18*r**2*ax*ay*hx*hy-
 255               44*hx**2*cy*by*ay*gy+44*hx**2*cy*by*ay*dy-8*hx**2*cx*bx*ay*gy+8*hx**2*cx*bx*ay*dy+
 256               8*hx**2*cx*bx*cy*by-12*hx**2*cx*ax*by*gy+12*hx**2*cx*ax*by*dy-24*hx**2*bx*ax*cy*gy+
 257               24*hx**2*bx*ax*cy*dy+8*hy*hx*bx*by**2*gy-8*hy*hx*bx*by**2*dy+18*hy*hx*ay**2*dx*gy-
 258               14*hy*hx*ay*cy**2*cx-18*hy*hx*ay**2*gx*gy-18*hy*hx*ay**2*dx*dy-16*hy*hx*cy*by**2*cx+
 259               18*hy*hx*ay**2*gx*dy-10*hy*hx*by*cy**2*bx+44*hy*hx*ay*cy*by*gx-44*hy*hx*ay*cy*by*dx+
 260               8*hy*hx*cy*by*ax*gy-8*hy*hx*cy*by*ax*dy+12*hy*hx*bx*cy*ay*gy-12*hy*hx*bx*cy*ay*dy+
 261               24*hy*hx*by*ay*cx*gy-24*hy*hx*by*ay*cx*dy-8*hy**2*cy*by*ax*gx+8*hy**2*cy*by*ax*dx+
 262               8*hy**2*cx*bx*cy*by-12*hy**2*cy*ay*bx*gx+12*hy**2*cy*ay*bx*dx-24*hy**2*by*ay*cx*gx+
 263               24*hy**2*by*ay*cx*dx-44*hy**2*cx*bx*ax*gx+44*hy**2*cx*bx*ax*dx-2*hy*hx*ay*cx**3+
 264               8*hy*hx*bx**3*gy-8*hy*hx*bx**3*dy+8*hy*hx*ay*cx*bx*gx-8*hy*hx*ay*cx*bx*dx-
 265               10*hy*hx*by*cx**2*bx+12*hy*hx*by*cx*ax*gx-12*hy*hx*by*cx*ax*dx+8*hy*hx*by*bx**2*gx-
 266               8*hy*hx*by*bx**2*dx-14*hy*hx*cy*cx**2*ax-16*hy*hx*cy*cx*bx**2+24*hy*hx*cy*bx*ax*gx-
 267               24*hy*hx*cy*bx*ax*dx+44*hy*hx*cx*bx*ax*gy-44*hy*hx*cx*bx*ax*dy-18*hy*hx*ax**2*gx*gy+
 268               18*hy*hx*ax**2*gx*dy+18*hy*hx*ax**2*dx*gy-18*hy*hx*ax**2*dx*dy,
 269
 270               6*hx**2*cy**3*by-24*hx**2*by*ay*gy*dy-8*hy**2*cy*by*bx*gx-14*hy**2*cx**2*ax*gx-
 271               12*hy*hx*r**2*bx*ay+2*hy*hx*cx**2*ay*gx-12*r**2*bx*ax*hx**2-12*r**2*by*ay*hy**2+
 272               16*hy**2*cx*bx**2*dx+6*hy**2*cx**3*bx-16*hy**2*cx*bx**2*gx+14*hy**2*cx**2*ax*dx+
 273               12*hy**2*bx*ax*gx**2+12*hy**2*bx*ax*dx**2-2*hy*hx*by*cx**3+4*hx**2*cx*bx*cy**2-
 274               2*hx**2*cx**2*ay*gy+2*hx**2*cx**2*ay*dy+2*hx**2*cx**2*cy*by+12*hx**2*bx*ax*gy**2+
 275               12*hx**2*bx*ax*dy**2-8*hx**2*cy*bx**2*gy+8*hx**2*cy*bx**2*dy+16*hx**2*cy*by**2*dy-
 276               14*hx**2*cy**2*ay*gy+14*hx**2*cy**2*ay*dy+12*hx**2*by*ay*gy**2+12*hx**2*by*ay*dy**2+
 277               4*hy**2*cx**2*cy*by-2*hy**2*cy**2*ax*gx+2*hy**2*cy**2*ax*dx+2*hy**2*cx*bx*cy**2+
 278               12*hy**2*by*ay*gx**2+12*hy**2*by*ay*dx**2-8*hy**2*cx*by**2*gx+8*hy**2*cx*by**2*dx-
 279               24*hy**2*bx*ax*gx*dx-2*hy*hx*cx**2*ay*dx+8*hy*hx*by*cx*bx*gx-8*hy*hx*by*cx*bx*dx-
 280               10*hy*hx*cy*cx**2*bx+12*hy*hx*cy*cx*ax*gx-12*hy*hx*cy*cx*ax*dx+8*hy*hx*cy*bx**2*gx-
 281               8*hy*hx*cy*bx**2*dx+14*hy*hx*cx**2*ax*gy-14*hy*hx*cx**2*ax*dy+16*hy*hx*cx*bx**2*gy-
 282               16*hy*hx*cx*bx**2*dy-24*hy*hx*bx*ax*gx*gy+24*hy*hx*bx*ax*gx*dy+24*hy*hx*bx*ax*dx*gy-
 283               24*hy*hx*bx*ax*dx*dy-2*hy*hx*cy**3*bx+14*hy*hx*cy**2*ay*gx-14*hy*hx*cy**2*ay*dx-
 284               10*hy*hx*by*cy**2*cx+16*hy*hx*cy*by**2*gx-16*hy*hx*cy*by**2*dx+2*hy*hx*cy**2*ax*gy-
 285               2*hy*hx*cy**2*ax*dy+8*hy*hx*cy*by*bx*gy-8*hy*hx*cy*by*bx*dy+12*hy*hx*cx*cy*ay*gy-
 286               12*hy*hx*cx*cy*ay*dy+8*hy*hx*cx*by**2*gy-8*hy*hx*cx*by**2*dy-24*hy*hx*by*ay*gx*gy+
 287               24*hy*hx*by*ay*gx*dy+24*hy*hx*by*ay*dx*gy-24*hy*hx*by*ay*dx*dy-8*hx**2*cx*bx*by*gy+
 288               8*hx**2*cx*bx*by*dy-24*hx**2*bx*ax*gy*dy-12*hx**2*cy*cx*ax*gy+12*hx**2*cy*cx*ax*dy-
 289               16*hx**2*cy*by**2*gy+8*hy**2*cy*by*bx*dx-24*hy**2*by*ay*gx*dx-12*hy**2*cx*cy*ay*gx+
 290               12*hy**2*cx*cy*ay*dx-12*hy*hx*r**2*ax*by,
 291
 292               hy**2*cy**2*cx**2-4*r**2*hy**2*by**2+4*hx**2*bx**2*dy**2+4*hx**2*by**2*gy**2-
 293               4*r**2*hx**2*bx**2+4*hy**2*by**2*dx**2+4*hy**2*by**2*gx**2+4*hx**2*by**2*dy**2+
 294               hx**2*cy**2*cx**2+4*hy**2*bx**2*dx**2+8*hy**2*cy*by*cx*dx+4*hx**2*bx**2*gy**2+
 295               2*hy*hx*cy**2*bx*gy+4*hy**2*bx**2*gx**2-2*hy*hx*cy**2*bx*dy-10*hy*hx*cy**2*by*dx-
 296               8*hy**2*cy*by*cx*gx+10*hy*hx*cy**2*by*gx-6*hy*hx*r**2*cx*ay-10*hx**2*cy**2*by*gy+
 297               10*hx**2*cy**2*by*dy+6*hx**2*cy*ay*gy**2+6*hx**2*cy*ay*dy**2-8*hx**2*by**2*gy*dy-
 298               10*hy**2*cx**2*bx*gx+10*hy**2*cx**2*bx*dx+6*hy**2*cx*ax*gx**2+6*hy**2*cx*ax*dx**2-
 299               8*hy**2*bx**2*gx*dx-6*r**2*hx**2*cx*ax-2*hy**2*cy**2*bx*gx+2*hy**2*cy**2*bx*dx+
 300               6*hy**2*cy*ay*gx**2+6*hy**2*cy*ay*dx**2-8*hy**2*by**2*gx*dx-2*hy*hx*cy**3*cx-
 301               6*r**2*hy**2*cy*ay-2*hy*hx*cy*cx**3-2*hx**2*cx**2*by*gy+2*hx**2*cx**2*by*dy+
 302               6*hx**2*cx*ax*gy**2+6*hx**2*cx*ax*dy**2-8*hx**2*bx**2*gy*dy+hx**2*cy**4+
 303               hy**2*cx**4-12*hx**2*cy*ay*gy*dy-12*hy**2*cx*ax*gx*dx-8*hy*hx*r**2*bx*by-
 304               6*hy*hx*r**2*ax*cy-12*hy**2*cy*ay*gx*dx+8*hy*hx*cy*by*cx*gy-8*hy*hx*cy*by*cx*dy-
 305               12*hy*hx*cy*ay*gx*gy+12*hy*hx*cy*ay*gx*dy+12*hy*hx*cy*ay*dx*gy-
 306               12*hy*hx*cy*ay*dx*dy-8*hy*hx*by**2*gx*gy+8*hy*hx*by**2*gx*dy+8*hy*hx*by**2*dx*gy-
 307               8*hy*hx*by**2*dx*dy+2*hy*hx*cx**2*by*gx-2*hy*hx*cx**2*by*dx+8*hy*hx*cy*cx*bx*gx-
 308               8*hy*hx*cy*cx*bx*dx+10*hy*hx*cx**2*bx*gy-10*hy*hx*cx**2*bx*dy-12*hy*hx*cx*ax*gx*gy+
 309               12*hy*hx*cx*ax*gx*dy+12*hy*hx*cx*ax*dx*gy-12*hy*hx*cx*ax*dx*dy-8*hy*hx*bx**2*gx*gy+
 310               8*hy*hx*bx**2*gx*dy+8*hy*hx*bx**2*dx*gy-8*hy*hx*bx**2*dx*dy-8*hx**2*cx*bx*cy*gy+
 311               8*hx**2*cx*bx*cy*dy-12*hx**2*cx*ax*gy*dy,
 312
 313               -2*hx**2*cy**3*gy-2*hy**2*cx**3*gx+2*hy**2*cx**3*dx-4*hy*hx*r**2*cx*by+
 314               2*hy*hx*cy*cx**2*gx+2*hx**2*cy**3*dy-4*r**2*cx*bx*hx**2-4*r**2*cy*by*hy**2-
 315               2*hy**2*cy**2*cx*gx+2*hy**2*cy**2*cx*dx+4*hy**2*cy*by*gx**2+4*hy**2*cy*by*dx**2+
 316               4*hx**2*cy*by*gy**2+4*hx**2*cy*by*dy**2+2*hy*hx*cx**3*gy-2*hy*hx*cx**3*dy+
 317               4*hy**2*cx*bx*gx**2+4*hy**2*cx*bx*dx**2-2*hx**2*cy*cx**2*gy+2*hx**2*cy*cx**2*dy+
 318               4*hx**2*cx*bx*gy**2+4*hx**2*cx*bx*dy**2+2*hy*hx*cy**3*gx-2*hy*hx*cy**3*dx-
 319               8*hy**2*cy*by*gx*dx-4*hy*hx*r**2*bx*cy-8*hx**2*cy*by*gy*dy-2*hy*hx*cy*cx**2*dx-
 320               8*hy*hx*cx*bx*gx*gy+8*hy*hx*cx*bx*gx*dy+8*hy*hx*cx*bx*dx*gy-
 321               8*hy*hx*cx*bx*dx*dy-8*hy**2*cx*bx*gx*dx-8*hx**2*cx*bx*gy*dy+2*hy*hx*cy**2*cx*gy-
 322               2*hy*hx*cy**2*cx*dy-8*hy*hx*cy*by*gx*gy+8*hy*hx*cy*by*gx*dy+
 323               8*hy*hx*cy*by*dx*gy-8*hy*hx*cy*by*dx*dy,
 324
 325               cx**2*hy**2*gx**2-2*cx**2*hy**2*gx*dx+cx**2*hy**2*dx**2-2*r**2*hx*hy*cx*cy+
 326               cy**2*hx**2*gy**2-2*cy**2*hx**2*gy*dy+cy**2*hx**2*dy**2+cy**2*hy**2*gx**2-
 327               2*cy**2*hy**2*gx*dx+cy**2*hy**2*dx**2-2*cy**2*hx*hy*gy*gx+2*cy**2*hx*hy*dy*gx+
 328               2*cy**2*hx*hy*gy*dx-2*cy**2*hx*hy*dy*dx+cx**2*hx**2*gy**2-2*cx**2*hx**2*gy*dy+
 329               cx**2*hx**2*dy**2-r**2*hy**2*cy**2-2*cx**2*hx*hy*gy*gx+2*cx**2*hx*hy*dy*gx+
 330               2*cx**2*hx*hy*gy*dx-2*cx**2*hx*hy*dy*dx-r**2*hx**2*cx**2]
 331
 332     # def z(t):
 333     #     res = 0
 334     #     for x in coeffs:
 335     #         res = res*t + x
 336     #     return res
 337
 338     # def zdot(t):
 339     #     res = 0
 340     #     i = len(coeffs)
 341     #     for x in coeffs[:-1]:
 342     #         res = res*t + (i-1)*x
 343     #         i -= 1
 344     #     return res
 345
 346     # def newton(t):
 347     #     # newton iteration to find the zero of z with t being the starting value for t
 348     #     if t < 0: t = 0
 349     #     if t > 1: t = 1
 350     #     slow = 1
 351     #     try:
 352     #         isz=z(t)
 353     #     except (ValueError, ArithmeticError):
 354     #         return
 355     #     loop = 0
 356     #     while abs(isz) > 1e-6:
 357     #         try:
 358     #             nt = t - slow * isz/zdot(t)
 359     #             newisz = z(nt)
 360     #         except (ValueError, ArithmeticError):
 361     #             slow *= 0.5 # we might have slow down the newton iteration
 362     #         else:
 363     #             t = nt
 364     #             isz = newisz
 365     #             if slow <= 0.5:
 366     #                 slow *= 2
 367     #         if loop > 100:
 368     #             break
 369     #         loop += 1
 370     #     else:
 371     #         return t
 372
 373     # def ts():
 374     #     epsilon = 1e-5
 375     #     res = []
 376     #     count = 100
 377     #     for x in range(count+1):
 378     #        t = newton(x/float(count))
 379     #        if t is not None and 0-0.5*epsilon < t < 1+0.5*epsilon:
 380     #            res.append(t)
 381     #     res.sort()
 382     #     i = 1
 383     #     while len(res) > i:
 384     #         if res[i] > res[i-1] + epsilon:
 385     #             i += 1
 386     #         else:
 387     #             del res[i]
 388     #     return res
 389
 390     # use Numeric to find the roots (via an equivalent eigenvalue problem)
 391     import Numeric, LinearAlgebra
 392     mat = Numeric.zeros((10, 10), Numeric.Float)
 393     for i in range(9):
 394         mat[i+1][i] = 1
 395     for i in range(10):
 396         mat[0][i] = -coeffs[i+1]/coeffs[0]
 397     ists = [zero.real for zero in LinearAlgebra.eigenvalues(mat) if -1e-10 < zero.imag < 1e-10 and 0 <= zero.real <= 1]
 398
 399     for t in ists:
 400         isl = (xdot(t)*(x(t)-gx)+ydot(t)*(y(t)-gy))/(xdot(t)*cos+ydot(t)*sin)
 401         c.fill(path.circle_pt(x(t), y(t), 1), [color.rgb.green])
 402         c.stroke(path.circle_pt(gx+isl*cos, gy+isl*sin, r), [color.rgb.green])
 403
 404
 405 # parameters for the enlargement:
 406 enlargeby_pt = 10
 407 round = 0 # a boolean to turn round corners on and off
 408
 409 # some examples:
 410 for anormpath in [path.rect(0, 0, 5, 5).normpath(),
 411                  path.circle(10, 3, 3).normpath(),
 412                  path.path(path.moveto(0, 8),
 413                            path.curveto(1, 8, 1, 9, 2, 11),
 414                            path.curveto(1, 12, 1, 12, 0, 11),
 415                            path.lineto(0, 8),
 416                            path.closepath()).normpath(),
 417                  path.path(path.moveto(5, 8),
 418                            path.curveto(20, 8, 6, 9, 7, 11),
 419                            path.curveto(6, 12, 6, 12, 5, 11),
 420                            path.lineto(5, 8),
 421                            path.closepath()).normpath(),
 422                  path.path(path.moveto(16, 0),
 423                            path.curveto(18, 0, 18, 4, 20, 4),
 424                            path.lineto(20, 12),
 425                            path.curveto(12, 12, 20, 8, 16, 8),
 426                            path.closepath()).normpath()]:
 427     c.stroke(anormpath)
 428     c.stroke(enlarged(anormpath, enlargeby_pt, round), [color.rgb.blue])
 429
 430 c.writeEPSfile("newbox")