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

# Untitled

a guest Sep 13th, 2013 50 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
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()
RAW Paste Data
Top