-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcomp.py
116 lines (98 loc) · 5.44 KB
/
comp.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
import scipy.integrate as integrate
from pylab import *
class Pendulum(object):
def __init__(self):
self.step = 0.01
self.Dr = 0 # sinusoidal driving coefficient
self.Da = 0 # damping coefficient
self.mr = 0.5 # mass of pendulum arm
self.md = 0.05 # mass of disk
self.l = 1 # length from disk to pivot
self.R = 0.025 # radius of the disk
self.B = 0.005 # magnetic field
self.Res = 1 #resistivity of disk metal
self.g = 9.8 #gravity
self.duration = 100
self.C = 2 * self.R * self.R * math.pi
self.time_vector = arange(0, self.duration, self.step) # time step
self.omegaD = ((1/2*pi)*((1/2)*self.mr + self.md)*self.g*self.l/(((1/3)*self.mr + self.md)*self.l**2)) #driving frequency reminder: two different frequencies
self.A = math.pi * self.R * self.R # area of disc
self.sigma = 0.5 #gaussian coefficient = gaussian variance
self.oms = (self.Res * 10 ** -8) * (self.C / self.A) #resistance of disc metal (copper)
self.mu = 2.6 * self.sigma # gaussian shift = expected value
########################################################################################################################
def update_time_vector(self):
self.time_vector = arange(0, self.duration, self.step)
########################################################################################################################
#System function
def system_function(self,r, t):
theta = r[0]
self.omega = r[1]
ftheta = self.omega
fomega = self.AngularAcceleration(theta) - self.Damped() +self.Driving(t) - self.Balpha(theta)
return [ftheta, fomega]
#######################################################################################################################
#Moment of inertia
def Moment(self):
return ((1/3)*self.mr*self.l**2) + (self.md*self.l**2)+((1/2)*self.md*self.R**2)
#######################################################################################################################
#gaussian
def GaussianPositive(self,theta):
return (1 / self.sigma * (2 * math.pi) ** 0.5) * math.e ** (-0.5 * ((theta - self.mu) / self.sigma) ** 2)
def GaussianNegative(self, theta):
return -(1 / self.sigma * (2 * math.pi) ** 0.5) * math.e ** (-0.5 * ((self.mu + theta) / self.sigma) ** 2)
def gaussian(self,theta):
if theta < 0:
return self.GaussianPositive(theta)
else:
return self.GaussianNegative(theta)
#######################################################################################################################
#induced current
def induced_current(self,theta, omega):
return -(self.B * 0.5 * self.R ** 2 * self.omega) / self.oms
########################################################################################################################
# magnetic dampening/driving: angular acceleration due to mag field
def Balpha(self,theta):
return (self.induced_current(theta, self.omega) * self.B * self.A * math.sin(theta)) / (.5 *self.md * self.R ** 2) * self.gaussian(theta)
########################################################################################################################
#damping term
def Damped(self):
return self.Da * (self.omega)
########################################################################################################################
#driving term
def Driving(self, t):
return self.Dr * math.cos(self.omegaD * t)
########################################################################################################################
#compound pendulum intrinsic angular acceleration
def AngularAcceleration(self,theta):
return (-.5 * (self.mr + self.md)) * self.g * self.l * math.sin(theta) / (((1 / 3) * self.mr + self.md) * self.l ** 2 + 0.5 * self.md * self.R * self.R)
########################################################################################################################
#momentum
def Momentum(self,theta, omega):
return self.Moment()*sin(theta)*omega
########################################################################################################################
def calculate(self):
f_int = [3, 0.0]
kinematics = integrate.odeint(self.system_function, f_int, self.time_vector)
position = kinematics[:, 0]
AngularVelocity = kinematics[:, 1]
figure(1)
suptitle('Time vs Angular Velocity')
plot(self.time_vector,AngularVelocity)
figure(2)
suptitle('Position vs Momentum')
plot(position, self.Momentum(position ,AngularVelocity))
#######################################################################################################################
#Create graphs in polar
figure(3)
plt.subplot(211, projection='polar')
plt.title("A line plot of the pendulum moving on the polar axis,\n theta=Angular Velocity versus r=Time")
plot(AngularVelocity, self.time_vector)
plt.subplot(212, projection='polar',label="Angular Velocity vs Position")
plt.title("A line plot of the pendulum moving on the polar axis,\n Angular Velocity versus Position\n\n\n")
plot(AngularVelocity, position)
plt.tight_layout()
show()
if __name__ == '__main__':
s = Pendulum()
s.calculate()