 # vpython_spinning_and_flipping_T-handle_5.py

Aug 13th, 2021 (edited)
258
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. # 23_08_21
2.
3. '''
4. Originally typed in from from the video:
5.  "Modeling a Spinning and Flipping T-Handle With Spring Forces and Python"
6.   in Rhett Allain's Physics Explained Channel (c. May 2021)
7.
8. tested on: Linux Mint 20.2, Geany 1.36, python 3.8.10,
9. Web Epiphany, i5-7260U CPU @ 2.20GHz × 2, Iris Plus Graphics 640
10. vpython version: ['7.6.1', 'jupyter']
11.
12.
13. references to Classical Mechanics by John R. Taylor (2005)
14.
15. p. 83
16.
17.   Principle of Conservation of Momentum
18.
19.   Ṗ = Fext
20.
21.     𝐏 = ∑ mα vα   α = 1,..., N
22.
23.
24. p. 344 The Centrifugal Force
25.
26.   Fcf = m(𝛀 × r) × 𝛀
27.
28. p. 348 The Coriolis Force
29.
30.   Fcor = 2mV × 𝛀
31.
32.
33. p. 394-400 Euler's Equations
34.
35. p. 681 Continuum Mechanics
36. p. 715 Relation between Stress and Strain: Hooke's Law
37. p. 718 The Equation of Motion for an Elastic Solid
38.
39.
40. An Introduction to the Coriolis Force
41. by Henry M. Stommel & Dennis W. Moore (1989)
42. "Particles that are not at rest in the rotating reference frame...
43. find themselves acted upon by the Coriolis forces." (p.62)
44.
45. Particles that deviate from the equilibrium spring position find
46. themselves acted upon by the restoring spring forces.
47.
48. #todo:
49.
50. update omega!!
51.
52. move the object with the Coriolis forces and only use the spring force
53.      to hold the object together (and allow movement to measure velocity)
54.
55. plot Euler's equation output against the object's movement
56.
58.
59. '''
60.
61. from vpython import *
62. # vpython does this:
63. #   from math import *
64. #   from numpy import arange
65. #   from random import random
66. #   import colorsys
67. #     use colorsys directly, no need for list_to_vec(L), use vec(*L)
68.
69. import sys
70.
71. def main():
72.
73.     scene.width, scene.height = 1000, 800
74.     scene.align = 'left'
75.     scene.resizable = False
76.     scene.autoscale = False
77.     scene.autozoom = False
78.     scene.range = 0.08
79.
80.     PAUSED_AT_START = False
81.
82.     TBAR = False
83.
84.     TRAILS = False
85.
86.     GRAPH = True
87.     XMAX = 20
88.
89.     RATE = 500
90.
91.     # ~ SPRING_CONSTANT = 2000
92.     SPRING_CONSTANT = 166700 # MAX
93.
94.     DELTA_T = 0.0017
95.
96.     #todo
97.     OMEGA_LIST = ((3, 0, 0), (3, 0.2, 0), (0, 0.1, 3), (0, 3, 0.2))
98.     # ~ g.omega = vec(3, 0, 0) # stable (unperturbed) rotation
99.     g.omega = vec(3, 0.2, 0) # unstable (perturbed) rotation
100.     # ~ g.omega = vec(0, 0.1, 3) # stable circular rotation
101.     # ~ g.omega = vec(0, 3, 0.2) # stable with a wobble
102.
103.
104.     # THE MASSES
105.     posA = vec(0, 0.035, 0)
106.     posB = vec(0, -0.035, 0)
107.     posC = vec(0.05, 0, 0)
108.     # ~ posC = vec(0.07, 0, 0)
109.
110.     # m1 is the mass bottom of the T
111.     m1 = 0.3
112.     # m2 are the masses on the T-handle
113.     m2 = 0.4
114.
115.     # Center of Mass
116.     CoM = (m2*posA + m2*posB + m1*posC) / (m2 + m2 + m1)
117.     label(text='CoM', pos=vec(0,0,0), font='bold',
118.           height=11, color=color.red, box=False, opacity=0)
119.
120.
121.     retain = 17
123.     massA = sphere(pos=posA-CoM, radius=m2/60, color=color.green,
125.     massB = sphere(pos=posB-CoM, radius=m2/60, color=color.red,
127.     massC = sphere(pos=posC-CoM, radius=m1/60, color=color.blue,
129.
130.     massA.m = m2
131.     massB.m = m2
132.     massC.m = m1
133.
134.     # v = 𝝎  ⨯ r   # p. 338
135.     massA.v = cross(g.omega, massA.pos)
136.     massB.v = cross(g.omega, massB.pos)
137.     massC.v = cross(g.omega, massC.pos)
138.
139.     # p = mv   # p. 14
140.     massA.p = massA.m * massA.v
141.     massB.p = massB.m * massB.v
142.     massC.p = massC.m * massC.v
143.
144.     # THE VECTORS
145.     rAB = massB.pos - massA.pos
146.     rBC = massC.pos - massB.pos
147.     rCA = massA.pos - massC.pos
148.
149.
150.     if TBAR:
151.         posD = massA.pos + rAB/2.0
152.         cyl1 = cylinder(pos=massA.pos, axis=rAB, radius=m2/60)
153.         cyl2 = cylinder(pos=massC.pos, axis=posD-massC.pos, radius=m1/60)
154.     else:
155.         # MAKE THE SPRINGS
156.         coils = 25
157.         thick = 0.0005
159.         steel = vec(0.33, 0.33, 0.45)
160.         springAB = helix(pos=massA.pos, axis=rAB, coils=coils,
162.         springBC = helix(pos=massB.pos, axis=rBC, coils=coils,
164.         springCA = helix(pos=massC.pos, axis=rCA, coils=coils,
166.
167.     # "UN-STRETCHED" LENGTH OF THE SPRINGS (L₀)
168.     rAB_0 = rAB.mag
169.     rBC_0 = rBC.mag
170.     rCA_0 = rCA.mag
171.
172.     # BUILTIN COLORS:
173.     # black, white, red, green, blue, yellow, cyan, magenta, orange, purple
174.     if GRAPH:
175.         tgraph = graph(xtitle='Time (s)', width=880, height=800,
176.                        ytitle='',
177.                        xmin=0, xmax=XMAX,
178.                        #ymin=-8e-4, ymax=8e-4,
179.                        fast=False, align='right')
180.
181.         red_pen = gcurve(color=color.red, width=2)
182.         blue_pen = gcurve(color=color.blue, width=2)
183.         black_pen = gcurve(color=color.black, width=2)
184.         sumofforce = gcurve(color=color.black, label='Conbined force', width=2)
185.         gray_pen = gcurve(color=vec(0.5,0.5,0.5), width=2)
186.         centrifugal_pen = gcurve(color=color.black, label='centrifugal force', width=2)
187.         coriolis_pen = gcurve(color=color.red, label='Coriolis force', width=2)
188.         spring_pen = gcurve(color=vec(0.6,0.6,1), label='Spring force', width=2)
189.
190.         def reset_graph():
191.             red_pen.delete()
192.             blue_pen.delete()
193.             black_pen.delete()
194.             sumofforce.delete()
195.             gray_pen.delete()
196.             centrifugal_pen.delete()
197.             coriolis_pen.delete()
198.             spring_pen.delete()
199.             g.t = 0
200.         scene.append_to_caption(' '*25)
201.         button(text="Reset Graph", bind=reset_graph)
202.
203.
204.     def pause_run():
205.         g.paused = not g.paused
206.         pause.text='  RUN   ' if g.paused else 'PAUSE'
207.     scene.append_to_caption(' '*10)
208.     pause = button(bind=pause_run)
209.     g.paused = not PAUSED_AT_START
210.     pause_run()
211.
212.     g.t = 0
213.     dt = DELTA_T
214.     k = SPRING_CONSTANT
215.
216.     while True:
217.
218.         if g.paused: continue
219.
220.         rate(RATE)
221.
222.         omega = g.omega
223.
224.         # UPDATE VECTORS
225.         rAB = massB.pos - massA.pos
226.         rBC = massC.pos - massB.pos
227.         rCA = massA.pos - massC.pos
228.
229.         # UPDATE SPRING FORCE OF ONE MASS ON ANOTHER
230.         # Fs = -k * (|L| - L₀) * L.hat
231.         FAB = -k * (rAB.mag - rAB_0) * rAB.hat
232.         FBC = -k * (rBC.mag - rBC_0) * rBC.hat
233.         FCA = -k * (rCA.mag - rCA_0) * rCA.hat
234.
235.         # UPDATE MOMENTUM OF THE MASSES
236.         # p2 = p1 + Fnet * 𝛥t
237.         massA.p += (FCA - FAB) * dt
238.         massB.p += (FAB - FBC) * dt
239.         massC.p += (FBC - FCA) * dt
240.
241.         # UPDATE POSITIONS OF MASSES
242.         # r2 = r1 + p2/m * 𝛥t
243.         massA.pos += massA.p * dt / massA.m
244.         massB.pos += massB.p * dt / massB.m
245.         massC.pos += massC.p * dt / massC.m
246.
247.         # UPDATE VELOCITIES
248.         # v = 𝝎  ⨯ r   # p. 338
249.         massA.v = cross(omega, massA.pos)
250.         massB.v = cross(omega, massB.pos)
251.         massC.v = cross(omega, massC.pos)
252.
253.         # UPDATE CORIOLIS FORCES
254.         corA = Coriolis_force(massA.m, omega, massA.v)
255.         corB = Coriolis_force(massB.m, omega, massB.v)
256.         corC = Coriolis_force(massC.m, omega, massC.v)
257.
258.         # UPDATE CENTRIFUGAL FORCES
259.         cenA = centrifugal_force(massA.m, omega, massA.pos)
260.         cenB = centrifugal_force(massB.m, omega, massB.pos)
261.         cenC = centrifugal_force(massC.m, omega, massC.pos)
262.
263.         sum_of_forcesA = cenA + corA
264.
265.         if TBAR:
266.             posD = massA.pos + rAB/2.0
267.             cyl1.pos = massA.pos
268.             cyl1.axis = rAB
269.             cyl2.pos = massC.pos
270.             cyl2.axis = posD-massC.pos
271.         else:
272.             # UPDATE SPRING POSITIONS AND LENGTHS
273.             springAB.pos = massA.pos; springAB.axis = rAB
274.             springBC.pos = massB.pos; springBC.axis = rBC
275.             springCA.pos = massC.pos; springCA.axis = rCA
276.
277.         # ℓ = r ⨯ p  # p.90
278.         # LA = cross(massA.pos, massA.p)
279.         # LB = cross(massB.pos, massB.p)
280.         # LC = cross(massC.pos, massC.p)
281.
282.         # black_pen, red_pen, blue_pen, gray_pen, magenta_pen,
283.         # cyan_pen, orange_pen, purple_pen
284.         # centrifugal_pen, coriolis_pen
285.         if GRAPH:
286.             # PLOT SELECTED STUFF
287.             if g.t <= 20.01:
288.                 # ~ black_pen.plot(g.t, LA.z)
289.                 # ~ red_pen.plot(g.t, corA.z)
290.                 # ~ black_pen.plot(g.t, cenA.z)
291.                 spring_pen.plot(g.t, FAB.z)
292.                 sumofforce.plot(g.t, sum_of_forcesA.z)
293.                 coriolis_pen.plot(g.t, corA.z)
294.
295.         g.t += dt
296.
297.
298. def centrifugal_force(m, 𝛀, r):
299.
300.     """Find the vector of the centrifugal force given
301.    mass (m), angular velocity vector (𝛀), and position vector (r).
302.
303.    Fcf = m(𝛀 ⨯ r) ⨯ 𝛀    # Taylor p.344  """
304.
305.     return m * cross(cross(𝛀, r), 𝛀)
306.
307.
308. def Coriolis_force(m, 𝛀, v):
309.
310.     """Find the vector of the Coriolis force given
311.    mass (m), angular velocity vector (𝛀), and velocity vector (v).
312.
313.    Fcor = 2mṙ ⨯ 𝛀  = 2mv ⨯ 𝛀    # Taylor p.348  """
314.
315.     return cross(2.0*m*v, 𝛀)
316.
317.
318. class g: ...
319.
320.
321. if __name__ == '__main__':
322.     sys.exit(main())
323.
324. '''
325.    L = OMEGA_LIST #todo
326.
327.    def M(m):
328.        val = m.selected
329.        if val == str(L):
330.            g.omega = vec(*L)
331.        elif val == str(L):
332.            g.omega = vec(*L)
333.        elif val == str(L):
334.            g.omega = vec(*L)
335.        elif val == str(L):
336.            g.omega = vec(*L)
337.        print(g.omega)
338.        #todo RESET POSITIONS,...
339.        g.paused = False
340.    scene.append_to_caption(' '*10)