In our last lesson we saw how we could take pretty much any complex number, plot its powers, and we would get a spiral. But these spirals all kind of looked the same and they flew out into space really fast. Let's try to make something that doesn't fly away so fast.
We're going to do something that might seem a little random at first, but be patient. Here is our plan. We will plot points so that:
The n-th point is distance n away from zero.
We rotate by the "Golden Ratio" each time - 2πφ radians.
Before we do that, we need to get better at plotting.
We'll be plotting a lot, so let's write a function to plot complex numbers and test it out on something simple:
importmathimportnumpyasnpimportmatplotlib.pyplotaspltdefplotComplex(cplxNums):pointArray=forcurrentPointincplxNums:pointArray.append([currentPoint.real,currentPoint.imag])points=np.array(pointArray)# the next lines look a little magic. LEarn more here: https://docs.scipy.org/doc/numpy/reference/arrays.indexing.htmlxPoints=points[:,0]yPoints=points[:,1]plt.figure(figsize=(8,8))plt.plot(xPoints,yPoints,'bo')plt.show()# now let's make some points and plot somethinggamma=-0.5-math.sqrt(3)*1j/2points=foriinrange(-10,10):forjinrange(-10,10):points.append(i+j*gamma)plotComplex(points)
The numbers we plotted here are Eisentstein Integers. Look them up! Can you plot Gaussian Integers? Any other quadratic integers?
OK, back to work. Let's rotate by Golden Ratio and move out by n and see what happens.
Does this look familiar? It's not a spiral, but do you see spirals?
Do the following exercises:
Plot 50 points and mark the clockwise spirals you see. Count them. Mark the counterclockwise spirals. Count them. (You can save the picture and open it in an image editor to do this). How many spirals were there in each direction?
Do the same for 100 points. Is it different? If so, are the old spirals still there?
Now 1000 points. Is it different?
It makes it a little easier to see these if we draw the Voronoi diagram for the points. The Voronoi diagram for a set of points comes from asking a simple question: for any point in the plane, what is the closest point from the set? We can think that a point in our set "owns" all of the points that are closer to it that to the other points. Now we can break the plane into regions owned by points in our list. Understanding these constructions is part of the field of Computational Geometry.
Let's compute the same point set, but now use the scipy Voronoi function to compute a Voronoi diagram for us.
fromscipy.spatialimportVoronoi,voronoi_plot_2ddefcomplexToArray(z):return[z.real,z.imag]# map our array of complex number into an array of pairs of pointspointArray=list(map(complexToArray,points))vor=Voronoi(pointArray)fig=voronoi_plot_2d(vor,show_vertices=False)plt.show()
To understand what the Voronoi function does for us a little better, let's take a simpler example with just a few points ant print the whole data structure
smallPoints=currentRotation=1foriinrange(1,9):nextPoint=currentRotation*math.sqrt(i)smallPoints.append(nextPoint)currentRotation=currentRotation*rotation# map our array of complex number into an array of pairs of pointssmallPointArray=list(map(complexToArray,smallPoints))vor=Voronoi(smallPointArray)print(vor.points)print(vor.regions)print(vor.vertices)print(vor.point_region)
The points are just the points we put in. What are the vertices and regions? We'll talk about this in class.
Let's use these structures to color our diagram.
plt.figure(figsize=(8,8))# let's make a color map with standard HTML color namescolors=['aqua','aquamarine','azure','bisque','blue','burlywood','cadetblue','coral','cornflowerblue','darkorange','darksalmon','darkseagreen','darkslategray','deepskyblue','dimgray','floralwhite','gainsboro','goldenrod','indianred','khaki','lavender','lightseagreen','lightskyblue','lightslategray','lightsteelblue','mediumseagreen','mediumslateblue','midnightblue','Moccasin','navajowhite','navy','oldlace','orange','peru','saddlebrown']vor=Voronoi(pointArray)colorModulus=21foriinrange(1,1490):region=vor.regions[vor.point_region[i]]polygon=[vor.vertices[i]foriinregion]ifnot-1inregion:plt.fill(*zip(*polygon),colors[i%colorModulus])plt.show()