Pastebin launched a little side project called VERYVIRAL.com, check it out ;-) Want more features on Pastebin? Sign Up, it's FREE!
Guest

Simulated pendulum

By: a guest on May 27th, 2013  |  syntax: Python  |  size: 3.01 KB  |  views: 38  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  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)