Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement