1 import sys
; sys
.path
.insert(0, "../..")
5 from pyx
import normpath
7 # make math.sqrt raising an exception for negativ values
13 math
._sqrt
= math
.sqrt
16 raise ValueError("sqrt of negativ number")
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:
33 # calculate crossing point of normal
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
)]
42 # since t1 is a translation + rotation, param is the radius
43 # hence, param should be equal to enlargeby_pt
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
:
51 return [path
._arctobcurve
(x
, y
, param
, angle1
, angle2
)]
53 # calculate crossing point of tangents
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
)]
62 x
, y
= t1
.apply_pt(param
, 0)
63 return [path
.normline_pt(xs
, ys
, x
, y
), path
.normline_pt(x
, y
, xe
, ye
)]
66 def enlarged(anormpath
, enlargeby_pt
, round):
68 for normsubpath
in anormpath
.normsubpaths
:
69 splitnormsubpathitems
= normsubpath
.normsubpathitems
[:] # do splitting on a copy
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
)
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
)
85 if isinstance(normsubpathitem
, normpath
.normcurve_pt
):
86 # We should do not alter the sign. Could we do any better here?
88 cs
= 1/(normsubpathitem
.curveradius_pt([0])[0] + enlargeby_pt
)
89 except ArithmeticError:
92 ce
= 1/(normsubpathitem
.curveradius_pt([1])[0] + enlargeby_pt
)
93 except ArithmeticError:
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]),
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
112 newnormpathitems
= path
.normline_pt(nxs
, nys
, nxe
, nye
)
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
)
123 def showtangent(normcurve_pt
):
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
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
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):
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
])
153 def showcircle(normcurve_pt
):
158 cos
= hx
/math
.sqrt(hx
*hx
+hy
*hy
)
159 sin
= hy
/math
.sqrt(hx
*hx
+hy
*hy
)
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
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
174 return ax
*t
*t
*t
+bx
*t
*t
+cx
*t
+normcurve_pt
.x0_pt
176 return ay
*t
*t
*t
+by
*t
*t
+cy
*t
+normcurve_pt
.y0_pt
178 return 3*ax
*t
*t
+2*bx
*t
+cx
180 return 3*ay
*t
*t
+2*by
*t
+cy
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
185 # the polynom is of order 10. its coefficients are:
187 coeffs
= [9*(ax
**2+ay
**2)*(hy
*ax
-ay
*hx
)**2, # * t**10,
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.
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,
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
,
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
,
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
,
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
,
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
,
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
,
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
,
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]
341 # for x in coeffs[:-1]:
342 # res = res*t + (i-1)*x
347 # # newton iteration to find the zero of z with t being the starting value for t
353 # except (ValueError, ArithmeticError):
356 # while abs(isz) > 1e-6:
358 # nt = t - slow * isz/zdot(t)
360 # except (ValueError, ArithmeticError):
361 # slow *= 0.5 # we might have slow down the newton iteration
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:
383 # while len(res) > i:
384 # if res[i] > res[i-1] + epsilon:
390 # use Numeric to find the roots (via an equivalent eigenvalue problem)
391 import Numeric
, LinearAlgebra
392 mat
= Numeric
.zeros((10, 10), Numeric
.Float
)
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]
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
])
405 # parameters for the enlargement:
407 round = 0 # a boolean to turn round corners on and off
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),
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),
421 path
.closepath()).normpath(),
422 path
.path(path
.moveto(16, 0),
423 path
.curveto(18, 0, 18, 4, 20, 4),
425 path
.curveto(12, 12, 20, 8, 16, 8),
426 path
.closepath()).normpath()]:
428 c
.stroke(enlarged(anormpath
, enlargeby_pt
, round), [color
.rgb
.blue
])
430 c
.writeEPSfile("newbox")