• API
• FAQ
• Tools
• Trends
• Archive
SHARE
TWEET

# Simulated pendulum

a guest May 27th, 2013 43 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. '''
2. Problem can be adimensionalized like this
3. m d^2x/dt^2 = m g - k (|x| - lo) x/|x| - d dx/dt
4.  d^2x/dt^2 =   g - k/m (|x| - lo) x/|x| - d/m dx/dt
5.
6. x' = x/lo
7. d^2x'/dt^2 = g - k/m (|x'| - 1) x'/|x'| - dlo/m dx'/dt
8.
9. t' = t*sqrt(k/m)
10. k/m d^2x'/dt'^2 = g - k/m (|x'| - 1) x'/|x'| - d lo sqrt(k/m)/m dx'/dt
11. d^2x'/dt'^2 = mg/k - (|x'| - 1) x'/|x'| - d lo sqrt(k m) dx'/dt
12.
13. so we just get two constants (plus the initial conditions):
14.  G = m*g/k  y  D = d lo sqrt(k m)
15.  d^2x'/dt'^2 = G - (|x'| - 1) x'/|x'| - D dx'/dt
16. '''
17. from __future__ import division
18. from visual import vector, dot
19. from math import *
20.
21. L0 = 1
22. k = 25
23. g = vector(0,-9.8,0)
24. mass = 0.2
25. vec_scale = 0.1  # How long is the vector for a 1N force
26. gamma = 100
27. tmax = 2000
28. dt = 0.0001
29.
30.
31. def get_components(vec, axis):
32.     axis_90 = (-axis[1], axis[0], 0)
33.     return dot(vec, axis), dot(vec, axis_90)
34.
35. # Create objects
36. class Point_Mass(object):
37.     def __init__(self, pos, velocity):
38.         self.pos = vector(pos)
39.         self.velocity = vector(velocity)
40.
41. def simulate(G, D, initial_length, theta0, end):
42.     pivot = Point_Mass(vector(0,L0,0), (0,0,0))
43.     bob = Point_Mass(pos=(initial_length * sin(theta0),
44.                           L0 - initial_length * cos(theta0),
45.                           0),
46.                      velocity=(0,0,0))
47.
48.     # First iteration step
49.     # Get forces
50.     fgrav = G
51.     axis = bob.pos - pivot.pos
52.     axis_n = axis.norm()
53.     deltal = axis.mag - 1
54.     if deltal > 0:
55.         fstring = -deltal*axis_n
56.     else:
57.         fstring = vector(0,0,0)
58.     v_rad, v_trans = get_components(bob.velocity, axis_n)
59.     fdamp = -D * v_rad * axis_n
60.     olddeltal = deltal
61.
62.     t = 0
63.     done = False
64.     final_pos = final_vel = None
65.     while t < tmax and not done:
66.         # Update objects
67.         # m*v = m*v0 + F_net * dt
68.         oldvelocity = bob.velocity
69.         bob.velocity = bob.velocity + (fgrav+fstring+fdamp)*dt
70.         # x = x0 + (m*v - m*v0)*dt/m/2
71.         bob.pos = bob.pos + 0.5*(bob.velocity + oldvelocity)*dt
72.
73.         # Update counters and indicators
74.         axis = bob.pos - pivot.pos
75.         axis_n = axis.norm()
76.         deltal = axis.mag - 1
77.         old_vtrans = v_trans
78.         v_rad, v_trans = get_components(bob.velocity, axis_n)
79.         if v_trans < 0 and old_vtrans > 0:
80.             print deltal
81.             if abs(deltal-olddeltal) < end:
82.                 done = True
83.                 final_pos = bob.pos
84.                 final_vel = bob.velocity
85.             olddeltal = deltal
86.         # Update force(s)
87.         if deltal > 0:
88.             fstring = -deltal*axis_n
89.         else:
90.             fstring = vector(0,0,0)
91.         fdamp = -D * v_rad * axis_n
92.         t += dt
93.     return deltal, final_pos, final_vel
94.
95.
96. G = mass*g/k
97. D = gamma*L0*sqrt(k*mass)
98. initial_length = (1+0.08476)*L0
99. angle = atan(0.60461/(L0-0.099364))
100. deltal, final_pos, final_vel =  simulate(G, D, initial_length, angle, 1e-5)
101. print (deltal, final_pos, final_vel)
102. print ("initial length: ", L0+deltal, "initial angle: ", angle)
RAW Paste Data
Top