Advertisement
Guest User

Untitled

a guest
Sep 13th, 2013
88
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 1.54 KB | None | 0 0
  1. from math import *
  2. from numpy import *
  3. from numpy.linalg import *
  4. from pylab import *
  5.  
  6. R = 20.0 # Radius of center of gravity
  7. w = 5.5  # Weight at center of gravity
  8. k = 1.0  # Spring constant
  9. D = 16.0 # Distance at which spring is connected
  10. x = 1.0  # Initial spring displacement
  11.  
  12. N = 200  # Number of incremental steps in curve simulation
  13.  
  14. # Vary any of the parameters by placing here
  15. for w in [4,5,6,7]:
  16.  
  17.     Q = array([0.0, 0.0])
  18.     D0 = norm(Q - array([D,0]))
  19.     x0 = x
  20.    
  21.     cam_pts = []
  22.     cir_pts = []
  23.  
  24.     for i in range(N):
  25.         ang0 = i * pi/2/N
  26.         ang1 = (i+1) * pi/2/N
  27.         P0 = D * array([cos(ang0), sin(ang0)])
  28.         P1 = D * array([cos(ang1), sin(ang1)])
  29.    
  30.         # Vector magic... http://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Vector_formulation
  31.         r = dot(P1-Q,P0-Q) / norm(P0-Q)**2
  32.         h = norm(r*(P0-Q) - (P1-Q))
  33.         l = r*norm(P0-Q)
  34.  
  35.         # Potential energy transferred from gravity into spring
  36.         K = sqrt( 2 * w*R/k * (cos(ang0)-cos(ang1)) + (norm(P0 - Q) + (x0-D0) )**2 ) - (x0-D0)
  37.  
  38.         if not l < K < l+h: break
  39.         #assert l<K<l+h
  40.  
  41.         # Move point the distance required to extend the spring by the proper amount
  42.         d = (l + K + h*h/(l-K))/2
  43.         Q = Q + d*(P0-Q)/norm(P0-Q)
  44.         x0 += d
  45.  
  46.         #print '(%f,%f),'%(Q[0],Q[1]),
  47.         cam_pts.append(Q)
  48.         cir_pts.append(P1)
  49.  
  50.         if not i % (N/20):
  51.             plot(*zip(Q,P1),color=(0,0,0,.1))
  52.    
  53.     plot(*zip(*cam_pts))
  54.     plot(*zip(*cir_pts))
  55.  
  56. print
  57. print 'Gravitational energy:', w * R
  58. print '       Spring energy:', .5 * k * (x0 + norm(P1-Q) - D0)**2
  59.  
  60. axes().set_aspect('equal','datalim')
  61. grid(1)
  62. show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement