/
7 code plus.py
183 lines (161 loc) · 7.81 KB
/
7 code plus.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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
# -*- coding: utf-8 -*-
"""
Created on Wed Oct 26 20:58:53 2016
@author: Administrator
"""
#Chapter 3 Oscillatory Motions and Chaos
#Section 3.1 & 3.2
#in SI
import pylab as pl
import math
import numpy as np
#import fractions as fr
class simple_harmonic_motion(object):
# def __init__(self, time_step = float(input('time step: ')), time_duration = float(input('time duration: ')), initial_theta = float(input('initial theta: ')), length = float(input('length: ')), strength_of_damping = float(input('stength of damping: ')), amplitude = float(input('amplitude of driving force: ')), anguluar_frequency = float(input('angular frequency of driving force: '))):
def __init__(self, time_step = 0.04, time_duration = float(input('time duration: ')), initial_theta = float(input('initial theta: ')), length = 9.8, strength_of_damping = float(input('damping factor: ')), amplitude = float(input('amplitude of driving force: ')), anguluar_frequency = 0.66666666):
self.l = length
self.dt = time_step
self.T = time_duration
self.n_steps = int(self.T / self.dt + 1)
self.theta = [initial_theta]
self.omega = [0]
self.t = [0]
self.g = 9.8
self.q = strength_of_damping
self.F_D = amplitude #The unit of F_D is /s²
self.Omega_D = anguluar_frequency
def calculate(self):
for i in range(self.n_steps):
self.omega.append(self.omega[i] - self.g / self.l * math.sin(self.theta[i]) * self.dt - self.q * self.omega[i] * self.dt + self.F_D * math.sin(self.Omega_D * self.t[i]) * self.dt)
self.theta.append(self.theta[i] + self.omega[i + 1] * self.dt)
# print("12")90
while(self.theta[i + 1] > math.pi):
self.theta[i + 1] -= 2 * math.pi
# print("123")
# if self.theta[i + 1] <= math.pi:
# break
# break
while self.theta[i + 1] < -math.pi:
self.theta[i + 1] += 2 * math.pi
# print("124")
# if self.theta[i + 1] >= -math.pi:
# break
# break
self.t.append(self.t[i] + self.dt)
global omega
omega = self.omega
global time_array
time_array = np.array(self.t)
global a
a = np.array(self.theta)
# print(a)
def calculate_delta(self):
# b= simple_harmonic_motion(time_duration = float(input('time duration: ')), initial_theta = float(input('initial theta: ')), amplitude = float(input('amplitude of driving force: ')))
b= simple_harmonic_motion()
b.calculate()
# self.theta_1 = b.calculate(self).self.theta
self.theta_1 = a
# print(self.theta_1)
b= simple_harmonic_motion(time_duration = float(input('time duration: ')), initial_theta = float(input('initial theta: ')), strength_of_damping = float(input('damping factor: ')), amplitude = float(input('amplitude of driving force: ')))
# b= simple_harmonic_motion()
b.calculate()
# self.theta_2 = b.calculate().self.theta
self.theta_2 = a
# print(self.theta_2)
self.delta = [abs(self.theta_1[0] - self.theta_2[0])]
# self.log_delta = [math.log10(abs(self.theta_1[0] - self.theta_2[0]))]
# self.delta = self.theta_1 - self.theta_2
self.time_array = time_array
for i in range(self.n_steps):
self.delta.append(abs(self.theta_1[i + 1] - self.theta_2[i + 1]))
# self.log_delta.append(math.log10(abs(self.theta_1[i + 1] - self.theta_2[i + 1])))
# print(self.delta)
# print(self.time_array)
def phase(self):
self.n = int(self.T / (3 * math.pi) + math.pi / 4)
self.time_phase = [0]
self.theta_phase = [0.2]
self.omega_phase = [0]
# print(self.t)
# print(self.omega)
# print(self.theta)
for i in range(self.n):
index = int((3 * (i + 1) * math.pi - math.pi / 4) / self.dt)
if abs(self.t[index] - 3 * (i + 1)) < self.dt / 2:
self.time_phase.append(self.t[index])
self.omega_phase.append(self.omega[index])
self.theta_phase.append(self.theta[index])
else:
self.time_phase.append(self.t[index + 1])
self.omega_phase.append(self.omega[index + 1])
self.theta_phase.append(self.theta[index + 1])
def show(self):
# pl.semilogy(self.theta, self.omega)
# , label = '$L =%.1f m, $'%self.l + '$dt = %.2f s, $'%self.dt + '$\\theta_0 = %.2f radians, $'%self.theta[0] + '$q = %i, $'%self.q + '$F_D = %.2f, $'%self.F_D + '$\\Omega_D = %.1f$'%self.Omega_D)
pl.plot(self.theta_phase ,self.omega_phase, '.', label = '$t \\approx 2\\pi n / \\Omega_D$')
pl.xlabel('$\\theta$ (radians)')
pl.ylabel('$\\omega$ (radians/s)')
pl.legend()
# pl.text(-1.4, 0.3, '$\\omega$ versus $\\theta$ $F_D = 1.2$', fontsize = 'x-large')
pl.title('Chaotic Regime')
# pl.show()
# pl.semilogy(self.time_array, self.delta)
# pl.legend(loc = 'upper center', fontsize = 'small')
# pl.xlabel('$time (s)$')
# pl.ylabel('$\\Delta\\theta (radians)$')
# pl.xlim(0, self.T)
# pl.ylim(float(input('ylim-: ')),float(input('ylim+: ')))
# pl.ylim(1E-11, 0.01)
# pl.text(4, -0.15, 'nonlinear pendulum - Euler-Cromer method')
# pl.text(10, 1E-3, '$\\Delta\\theta versus time F_D = 0.5$')
# pl.title('Simple Harmonic Motion')
pl.title('Chaotic Regime')
def show_log(self):
# pl.subplot(121)
pl.semilogy(self.time_array, self.delta, 'c')
pl.xlabel('$time (s)$')
pl.ylabel('$\\Delta\\theta$ (radians)')
pl.xlim(0, self.T)
# pl.ylim(1E-11, 0.01)
pl.text(42, 1E-7, '$\\Delta\\theta$ versus time $F_D = 1.2$', fontsize = 'x-large')
pl.title('Chaotic Regime')
pl.show()
# def show_log_sub122(self):
# pl.subplot(122)
# pl.semilogy(self.time_array, self.delta, 'g')
# pl.xlabel('$time (s)$')
# pl.ylabel('$\\Delta\\theta$ (radians)')
# pl.xlim(0, self.T)
# pl.ylim(1E-6, 100)
# pl.text(20, 1E-5, '$\\Delta\\theta$ versus time $F_D = 1.2$', fontsize = 'x-large')
# pl.title('Chaotic Regime')
# pl.show()
def multi_show(self):
for i in range(2):
a = simple_harmonic_motion(time_step = float(input('time step: ')), time_duration = float(input('time duration: ')), initial_theta = float(input('initial theta: ')), length = float(input('length: ')), strength_of_damping = float(input('stength of damping: ')), amplitude = float(input('amplitude of driving force: ')), anguluar_frequency = float(input('angular frequency of driving force: ')))
a.calculate()
a.show()
pl.show()
#class please_input():
# string_input = input('xlocation ,ylocation: ')
# numbers = [float(n) for n in string_input.split()]
# x = numbers[0]
# y = numbers[1]
#b = simple_harmonic_motion()
#b.calculate_delta()
#b.show()
s = simple_harmonic_motion()
s.calculate()
#s.show()
s.phase()
s.show()
#s = simple_harmonic_motion()
#s.calculate_delta()
#s.show_log()
#c = simple_harmonic_motion()
#c.calculate_delta()
#c.show_log()
#s = simple_harmonic_motion(time_duration = float(input('time duration: ')), initial_theta = float(input('initial theta: ')), strength_of_damping = float(input('damping factor: ')), amplitude = float(input('amplitude of driving force: ')))
##s = simple_harmonic_motion()
#s.calculate_delta()
#s.show_log_sub122()