# vpython_spring_centrifugal_coriolis_force_2.py

Jul 27th, 2021 (edited)
1,890
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1.
2.
3. '''
4. vpython_spring_centrifugal_coriolis_force_2.py
5. 30_07_21
6.
7. tested on: Linux Mint 20.2, Geany 1.36, python 3.8.10,
8. Web Epiphany, i5-7260U CPU @ 2.20GHz × 2, Iris Plus Graphics 640
9.
10. references:
11.
12. An Introduction to the Coriolis Force (1989)
13.    by Henry M. Stommel & Dennis W. Moore, CHAPTER III
14.
15. p.60:  3.4 EXPERTS, NOVICES AND HOOKE SPRINGS
16. "Now we attach a particle of unit mass to the origin by
17.    means of a massless spring whose tension is equal to
18.    -k²rⁿ where k² is the so-called spring constant..."
19.
20. Classical Mechanics by John R. Taylor (2005), pp. 342-
21.
22. An Introduction to Dynamic Meteorology, Holton & Hakim, Fifth Ed.
23. p.17: "The Coriolis force can only change the
24.        direction of motion, not the speed."
25.
26. #todo how to test p.17? graph v.mag?
27.
28.
29. vpython vector v = v.hat * v.mag
30.
31. copying:
33. v = vector(...)
34. then:
35. v.value = vector(...) "ensures a copy"
36. DON'T FALL INTO THIS TRAP (LIKE I DID):
37. v = u
38.
39. '''
40.
41. from math import radians, sin
42. import sys
43. from vpython import *
44.
45. I = vec(1, 0, 0)
46. J = vec(0, 1, 0)
47. K = vec(0, 0, 1)
48.
49.
50. def main():
51.
52.     impulse_shape = [sin(radians(x))*0.07 for x in range(0, 180, 10)]
53.
54.     init()
55.
56.     g.spring.constant = 0.03
57.     # mass of particle increased to bring out Coriolis effect
58.     g.ball.m = 5.0
59.
60.     g.ball.pos = vec(1, 0, 0)
61.     last_pos = vec(g.ball.pos)
62.
63.     spring_null = 1.0
64.     dt = 0.01
65.     idt2 = 1.0 / dt**2
66.
67.     r = g.ball.pos - g.center
68.     v = vec(0, 0, 0)
69.
70.     g.impulse = False
71.     impulse_index = 0
72.
73.     while 1:
74.
75.         rate(30)
76.
77.         if not g.run: continue
78.
79.         ω = vec(0.0, float(g.slider_text.text), 0.0)
80.         force = centrifugal_force(g.ball.m, ω, r)
81.
82.         if g.coriolis:
83.             force += coriolis_force(g.ball.m, ω, v)
84.
85.         extension = (r - (spring_null * r.hat))
86.         force += (-g.spring.constant * extension)
87.
88.         if g.impulse:
89.             force += g.sign * impulse_shape[impulse_index] * r.hat
90.             impulse_index += 1
91.             if impulse_index == len(impulse_shape):
92.                 impulse_index = 0
93.                 g.impulse = 0
94.
95.         if 1:
96.             a = force * idt2
97.             v += a * dt
98.             g.ball.pos += v * dt
99.         else:
100.             g.ball.pos += force
101.
102.         g.spring.axis = g.ball.pos + g.spring.pos
103.
104.         v = g.ball.pos - last_pos
105.         g.general_text.text = 'velocity ' + str(v.mag)
106.
107.         r = g.ball.pos - g.center
108.
109.         last_pos.value = g.ball.pos # "ensures a copy"
110.
111.         if g.rotate:
112.             angle = g.slider.value / 2
113.             g.disk.rotate(angle=angle, axis=J, origin=g.center)
114.             g.ball.rotate(angle=angle, axis=J, origin=g.center)
115.             g.spring.rotate(angle=angle, axis=J, origin=g.center)
116.
117.         #scene.waitfor('redraw')
118.
119.
120. def init():
121.
122.     scene.visible = False
123.
124.     g.center = vec(0, 0, 0)
125.
126.     scene.userspin = False
127.     scene.userzoom = False
128.     autoscale = False
129.     autozoom = False
130.     # these vary with browser for some reason
131.     scene.width, scene.height = 1900, 910
132.
133.     # THE "PLATFORM"
134.     g.disk = cylinder(axis=J,
135.                          size=vec(0.01, 6, 6),
136.                          pos=vec(0, -0.2, 0),
137.                          shininess=0.2,
138.                          texture=textures.wood)
139.
140.     # THE "HOOK SPRING"
141.     g.spring = helix(pos=g.center,
142.                         axis=I,
144.                         #size=vec(1, 0.2, 0.2),
145.                         coils=10,
146.                         thickness=0.02,
147.                         color=color.blue,
148.                         shininess=0.4)
149.
150.     # THE "PARTICLE"
151.     g.ball = sphere(size=vec(0.2, 0.2, 0.2), color=vec(0,0,0.4))
152.
153.     # SPINDLE
154.     cylinder(axis=J, size=vec(0.2, 0.2, 0.2), shininess=0.4,
155.             pos=vec(0, -0.1, 0), color=vec(0.33, 0.35, 0.37))
156.
157.     # DIRECTION
158.     arrow(pos=vec(0.5,0,-3.05), axis=vec(-1,0,0), shaftwidth=0.04,
159.             color=color.red, length=1)
160.
161.     # ANGULAR VELOCITY SLIDER
162.     def setomega(s):
163.         g.slider_text.text = '{:1.3f}'.format(s.value)
164.     scene.append_to_caption('\n', ' '*110)
165.     wtext(text='ANGULAR VELOCITY ')
166.     g.slider = slider(min=0.0, max=.08, value=0.04, length=600, bind=setomega)
167.
168.     # SLIDER TEXT
169.     scene.append_to_caption(' '*4)
170.     g.slider_text = wtext(text='{:1.2f}'.format(g.slider.value))
172.
173.     text = 'You are rotating with the disk.'
174.     r_text = label(text=text, pos=vec(-2.2, 3, 0),
175.         visible=True, height=25, box=False, color=vec(0,0,0.5))
176.
177.     def set_view_text():
178.         r_text.visible = False if g.rotate else True
179.
180.     # ROTATE
181.     g.rotate = False
182.     def toggle_rotate():
183.         g.rotate = not g.rotate
184.         set_view_text()
185.     scene.append_to_caption('\n'+' '*112)
186.     button(text="TOGGLE ROTATE", bind=toggle_rotate)
187.
188.     # VIEW
189.     g.view = 0
190.     def toggle_view():
191.         g.view = not g.view
192.         set_view_text()
193.         scene.range = 5.0
194.         if g.view: #elevation
195.             scene.camera.pos = vec(0.0, 5.5, 0.0)
196.             scene.camera.axis = vec(0.0, -1.7, -1e-9)
197.         else: # plan
198.             scene.camera.pos = vec(0.0, 1.05, 5.5)
199.             scene.camera.axis = vec(0.0, 0.0, -1.7)
200.     toggle_view()
201.     scene.append_to_caption(' '*5)
202.     button(text="TOGGLE VIEW", bind=toggle_view)
203.
204.     # TOGGLE CORIOLIS FORCE
205.     g.coriolis = True
206.     def togglecoriolis():
207.         g.coriolis = not g.coriolis
208.         toggle_coriolis.text = 'CORIOLIS OFF' if g.coriolis else ' CORIOLIS ON '
209.     scene.append_to_caption(' '*5)
210.     toggle_coriolis = button(bind=togglecoriolis)
211.     togglecoriolis()
212.
213.     # PAUSE, RUN
214.     g.run = False
215.     def runpause():
216.         g.run = not g.run
217.         run_pause.text = 'PAUSE' if g.run else '  RUN   '
218.     scene.append_to_caption(' '*5)
219.     run_pause = button(bind=runpause)
220.     runpause()
221.
222.     #APPLY FORCE OUT
223.     def applyforceout():
224.         g.impulse = True
225.         g.sign = 1
226.     scene.append_to_caption(' '*5)
227.     button(text='APPLY FORCE OUT', bind=applyforceout)
228.
229.     #APPLY FORCE IN
230.     def applyforcein():
231.         g.impulse = True
232.         g.sign = -1
233.     scene.append_to_caption(' '*5)
234.     button(text='APPLY FORCE IN', bind=applyforcein)
235.
236.     scene.append_to_caption(' '*5)
237.     g.general_text = wtext(text='####')
238.
239.     scene.waitfor('textures')
240.     scene.visible = True
241.
242.
243. def centrifugal_force(m, 𝛀, r):
244.
245.     """Find the vector of the centrifugal force given
246.    mass (m), angular velocity vector (𝛀), and position vector (r).
247.
248.    Fcf = m(𝛀 x r) x 𝛀    # Taylor p.344  """
249.
250.     return m * cross(cross(𝛀, r), 𝛀)
251.
252.
253. def coriolis_force(m, 𝛀, v):
254.
255.     """Find the vector of the coriolis force given
256.    mass (m), angular velocity vector (𝛀), and velocity vector (v).
257.
258.    Fcor = 2mṙ x 𝛀  = 2mv x 𝛀    # Taylor p.348  """
259.
260.     return cross(2.0*m*v, 𝛀)
261.
262.
263. class g: ...
264.
265.
266. if __name__ == '__main__':
267.     sys.exit(main())
268.
RAW Paste Data