1 from matplotlib.pylab import *
2
4 im1 = i-1
5 m1 = (y1[i] - y1[im1])/(x[i] - x[im1])
6 m2 = (y2[i] - y2[im1])/(x[i] - x[im1])
7 b1 = y1[im1] - m1*x[im1]
8 b2 = y2[im1] - m2*x[im1]
9 xZero = (b1 - b2)/(m2 - m1)
10 yZero = m1*xZero + b1
11 return (xZero, yZero)
12
14 diff = y2 - y1
15 pos = []
16 neg = []
17 xx1 = [x[0]]
18 xx2 = [x[0]]
19 yy1 = [y1[0]]
20 yy2 = [y2[0]]
21 oldSign = (diff[0] < 0 )
22 npts = x.shape[0]
23 for i in range(1,npts):
24 newSign = (diff[i] < 0)
25 if newSign != oldSign:
26 xz,yz = findZero(i,x,y1,y2)
27 xx1.append(xz)
28 yy1.append(yz)
29 xx2.reverse()
30 xx1.extend(xx2)
31 yy2.reverse()
32 yy1.extend(yy2)
33 if oldSign:
34 neg.append( (xx1,yy1) )
35 else:
36 pos.append( (xx1,yy1) )
37 xx1 = [xz,x[i]]
38 xx2 = [xz,x[i]]
39 yy1 = [yz,y1[i]]
40 yy2 = [yz,y2[i]]
41 oldSign = newSign
42 else:
43 xx1.append( x[i])
44 xx2.append( x[i])
45 yy1.append(y1[i])
46 yy2.append(y2[i])
47 if i == npts-1:
48 xx2.reverse()
49 xx1.extend(xx2)
50 yy2.reverse()
51 yy1.extend(yy2)
52 if oldSign :
53 neg.append( (xx1,yy1) )
54 else:
55 pos.append( (xx1,yy1) )
56 return pos,neg
57