windmill.py

import math, random
import matplotlib.pyplot as plt
import matplotlib.cm as cm
plt.style.use('dark_background')

def create_point():
    x = random.random()
    y = random.random()

    return x,y

def print_points(points):
    print("pts:")
    for point in points:
        print("    ",point)
    print("end")

def plot_points(points):
    for pt in points: plt.scatter(pt[0],pt[1])
    plt.axis([0, 1, 0, 1])
    plt.show()

pts_x = []
pts_y = []
points = 7

for i in range(points):
    pt_x,pt_y = create_point()
    pts_x.append(pt_x)
    pts_y.append(pt_y)

#print and show points
#plot_points(points)


# import packages
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

class windmill:
    """Gaussian Initial Condition, Fixed Boundary Condition"""
    def __init__(self,x0,y0,x_last,y_last,pts_x,pts_y,dt=0.01,length=1.5):
        self.x0=x0
        self.y0=y0
        self.x_last = x_last
        self.y_last = y_last
        self.theta0 = math.atan2(y_last-y0,x_last-x0)
        self.theta = self.theta0
        self.dt=dt
        self.length = length
        self.i=10000

        self.pts_x=pts_x
        self.pts_y=pts_y

        #find next point
        max_theta = 999
        for i in range(len(pts_x)):
            _x = pts_x[i]
            _y = pts_y[i]
            if ((_x is not self.x0) and (_y is not self.y0)) and ((_y is not self.y_last) and (_x is not self.x_last)):
                v1 = np.array([x_last-self.x0,y_last-self.y0])
                v2 = np.array([_x - self.x0, _y - self.y0])
                adotb = v1[0]*v2[0] + v1[1]*v2[1]
                _theta = math.acos(adotb / np.sqrt(v1.dot(v1)) / np.sqrt(v2.dot(v2)))

                if v1[1]*v2[0] > v1[0]*v2[1]: dir = "cw"
                else: dir = "ccw"
                if dir == "ccw": _theta = math.pi - _theta

                if _theta < max_theta:
                    max_theta = _theta
                    self.x1 = _x
                    self.y1 = _y
                    self.theta1 = self.theta0-_theta

        self.X= [np.linspace(x0-(math.cos(self.theta))*self.length,
                           x0+(math.cos(self.theta))*self.length,self.i,endpoint=True)]
        self.Y= [np.linspace(y0 - (math.sin(self.theta)) * self.length,
                           y0 + (math.sin(self.theta)) * self.length, self.i, endpoint=True)]
        return None

    def calculate(self):

        def find_end(x0,y0,x_last,y_last,theta0,pts_x,pts_y):
            max_theta = 999
            for i in range(len(pts_x)):
                _x = pts_x[i]
                _y = pts_y[i]
                if ((_x is not x0) and (_y is not y0)) and (
                        (_y is not y_last) and (_x is not x_last)):
                    v1 = np.array([x_last - x0, y_last - self.y0])
                    v2 = np.array([_x - x0, _y - self.y0])
                    adotb = v1[0] * v2[0] + v1[1] * v2[1]
                    _theta = math.acos(adotb / np.sqrt(v1.dot(v1)) / np.sqrt(v2.dot(v2)))

                    if v1[1] * v2[0] > v1[0] * v2[1]:
                        dir = "cw"
                    else:
                        dir = "ccw"
                    if dir == "ccw": _theta = math.pi - _theta

                    if _theta < max_theta:
                        max_theta = _theta
                        x1 = _x
                        y1 = _y
                        theta1 = theta0 - _theta
            return x1,y1,theta1

        for i in range(20):
            while self.theta > self.theta1:
                self.theta -= self.dt
                self.X.append(np.linspace(self.x0-(math.cos(self.theta))*self.length,
                               self.x0+(math.cos(self.theta))*self.length,self.i,endpoint=True))
                self.Y.append(np.linspace(self.y0 - (math.sin(self.theta)) * self.length,
                               self.y0 + (math.sin(self.theta)) * self.length, self.i, endpoint=True))
            self.x_last = self.x0
            self.y_last = self.y0
            self.x0 = self.x1
            self.y0 = self.y1

            #Find next point
            self.x1,self.y1,self.theta1 = find_end(self.x0,self.y0,self.x_last,self.y_last,self.theta,pts_x,pts_y)
        return None
    def plot(self,i):
        plt.plot(self.X[i],self.Y[i])
        return 0
    def movie(self):
        # New figure with white background
        fig = plt.subplots(figsize=(6,6), facecolor='white')
        # New axis over the whole figure, no frame and a 1:1 aspect ratio
        ax = fig.add_axes([0,0,1,1], frameon=False, aspect=1)
        line, = ax.plot([], [], lw=2)
        def init():
            line.set_data([], [])
            return line,
        def animate(i):
            x = self.X
            y = self.Y
            line.set_data(list(x), list(y))
            return line,

        #anim1=animation.FuncAnimation(fig, animate, init_func=init, frames=2500, interval=1,blit=True)

        #plt.show()
        return 0

A=windmill(pts_x[0],pts_y[0],pts_x[1],pts_y[1],pts_x,pts_y)
A.calculate()


# New figure with white background
fig = plt.figure(figsize=(6,6)) #, facecolor='#ffa600'
# New axis over the whole figure, no frame and a 1:1 aspect ratio
ax = plt.axes(xlim=(0, 1), ylim=(0, 1))
line, = ax.plot([], [], lw=2, zorder = 1, c='#003f5c')

def init():
    line.set_data([], [])
    return line,
def animate(i):
    x = A.X[i]
    y = A.Y[i]
    line.set_data(list(x), list(y))
    return line,

colors = ['#00876c',
'#459e70',
'#9fc878',
'#cedc80',
'#fdcd70',
'#f9ab5c',
'#e5644e',
'#d43d51']

scat = plt.scatter(pts_x,pts_y,c=colors[:len(pts_x)], zorder=2)
#scat.set_alpha(0.5)

anim1=animation.FuncAnimation(fig, animate, init_func=init, frames=len(A.X), interval=25, repeat_delay  = 100)
plt.xticks([], [])
plt.yticks([], [])
#plt.show()
anim1.save('windmill.gif', writer='pillow', fps=45)