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