Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- '''
- vpython_spring_centrifugal_coriolis_force_2.py
- 30_07_21
- tested on: Linux Mint 20.2, Geany 1.36, python 3.8.10,
- Web Epiphany, i5-7260U CPU @ 2.20GHz × 2, Iris Plus Graphics 640
- references:
- An Introduction to the Coriolis Force (1989)
- by Henry M. Stommel & Dennis W. Moore, CHAPTER III
- p.60: 3.4 EXPERTS, NOVICES AND HOOKE SPRINGS
- "Now we attach a particle of unit mass to the origin by
- means of a massless spring whose tension is equal to
- -k²rⁿ where k² is the so-called spring constant..."
- Classical Mechanics by John R. Taylor (2005), pp. 342-
- An Introduction to Dynamic Meteorology, Holton & Hakim, Fifth Ed.
- p.17: "The Coriolis force can only change the
- direction of motion, not the speed."
- #todo how to test p.17? graph v.mag?
- vpython vector v = v.hat * v.mag
- copying:
- # v must be a vector to start with
- v = vector(...)
- then:
- v.value = vector(...) "ensures a copy"
- DON'T FALL INTO THIS TRAP (LIKE I DID):
- v = u
- '''
- from math import radians, sin
- import sys
- from vpython import *
- I = vec(1, 0, 0)
- J = vec(0, 1, 0)
- K = vec(0, 0, 1)
- def main():
- impulse_shape = [sin(radians(x))*0.07 for x in range(0, 180, 10)]
- init()
- g.spring.constant = 0.03
- # mass of particle increased to bring out Coriolis effect
- g.ball.m = 5.0
- g.ball.pos = vec(1, 0, 0)
- last_pos = vec(g.ball.pos)
- spring_null = 1.0
- dt = 0.01
- idt2 = 1.0 / dt**2
- r = g.ball.pos - g.center
- v = vec(0, 0, 0)
- g.impulse = False
- impulse_index = 0
- while 1:
- rate(30)
- if not g.run: continue
- ω = vec(0.0, float(g.slider_text.text), 0.0)
- force = centrifugal_force(g.ball.m, ω, r)
- if g.coriolis:
- force += coriolis_force(g.ball.m, ω, v)
- extension = (r - (spring_null * r.hat))
- force += (-g.spring.constant * extension)
- if g.impulse:
- force += g.sign * impulse_shape[impulse_index] * r.hat
- impulse_index += 1
- if impulse_index == len(impulse_shape):
- impulse_index = 0
- g.impulse = 0
- if 1:
- a = force * idt2
- v += a * dt
- g.ball.pos += v * dt
- else:
- g.ball.pos += force
- g.spring.axis = g.ball.pos + g.spring.pos
- v = g.ball.pos - last_pos
- g.general_text.text = 'velocity ' + str(v.mag)
- r = g.ball.pos - g.center
- last_pos.value = g.ball.pos # "ensures a copy"
- if g.rotate:
- angle = g.slider.value / 2
- g.disk.rotate(angle=angle, axis=J, origin=g.center)
- g.ball.rotate(angle=angle, axis=J, origin=g.center)
- g.spring.rotate(angle=angle, axis=J, origin=g.center)
- #scene.waitfor('redraw')
- def init():
- scene.visible = False
- g.center = vec(0, 0, 0)
- scene.userspin = False
- scene.userzoom = False
- autoscale = False
- autozoom = False
- # these vary with browser for some reason
- scene.width, scene.height = 1900, 910
- # THE "PLATFORM"
- g.disk = cylinder(axis=J,
- size=vec(0.01, 6, 6),
- pos=vec(0, -0.2, 0),
- shininess=0.2,
- texture=textures.wood)
- # THE "HOOK SPRING"
- g.spring = helix(pos=g.center,
- axis=I,
- radius=0.07,
- #size=vec(1, 0.2, 0.2),
- coils=10,
- thickness=0.02,
- color=color.blue,
- shininess=0.4)
- # THE "PARTICLE"
- g.ball = sphere(size=vec(0.2, 0.2, 0.2), color=vec(0,0,0.4))
- # SPINDLE
- cylinder(axis=J, size=vec(0.2, 0.2, 0.2), shininess=0.4,
- pos=vec(0, -0.1, 0), color=vec(0.33, 0.35, 0.37))
- # DIRECTION
- arrow(pos=vec(0.5,0,-3.05), axis=vec(-1,0,0), shaftwidth=0.04,
- color=color.red, length=1)
- # ANGULAR VELOCITY SLIDER
- def setomega(s):
- g.slider_text.text = '{:1.3f}'.format(s.value)
- scene.append_to_caption('\n', ' '*110)
- wtext(text='ANGULAR VELOCITY ')
- g.slider = slider(min=0.0, max=.08, value=0.04, length=600, bind=setomega)
- # SLIDER TEXT
- scene.append_to_caption(' '*4)
- g.slider_text = wtext(text='{:1.2f}'.format(g.slider.value))
- scene.append_to_caption(' radians/s\n')
- text = 'You are rotating with the disk.'
- r_text = label(text=text, pos=vec(-2.2, 3, 0),
- visible=True, height=25, box=False, color=vec(0,0,0.5))
- def set_view_text():
- r_text.visible = False if g.rotate else True
- # ROTATE
- g.rotate = False
- def toggle_rotate():
- g.rotate = not g.rotate
- set_view_text()
- scene.append_to_caption('\n'+' '*112)
- button(text="TOGGLE ROTATE", bind=toggle_rotate)
- # VIEW
- g.view = 0
- def toggle_view():
- g.view = not g.view
- set_view_text()
- scene.range = 5.0
- if g.view: #elevation
- scene.camera.pos = vec(0.0, 5.5, 0.0)
- scene.camera.axis = vec(0.0, -1.7, -1e-9)
- else: # plan
- scene.camera.pos = vec(0.0, 1.05, 5.5)
- scene.camera.axis = vec(0.0, 0.0, -1.7)
- toggle_view()
- scene.append_to_caption(' '*5)
- button(text="TOGGLE VIEW", bind=toggle_view)
- # TOGGLE CORIOLIS FORCE
- g.coriolis = True
- def togglecoriolis():
- g.coriolis = not g.coriolis
- toggle_coriolis.text = 'CORIOLIS OFF' if g.coriolis else ' CORIOLIS ON '
- scene.append_to_caption(' '*5)
- toggle_coriolis = button(bind=togglecoriolis)
- togglecoriolis()
- # PAUSE, RUN
- g.run = False
- def runpause():
- g.run = not g.run
- run_pause.text = 'PAUSE' if g.run else ' RUN '
- scene.append_to_caption(' '*5)
- run_pause = button(bind=runpause)
- runpause()
- #APPLY FORCE OUT
- def applyforceout():
- g.impulse = True
- g.sign = 1
- scene.append_to_caption(' '*5)
- button(text='APPLY FORCE OUT', bind=applyforceout)
- #APPLY FORCE IN
- def applyforcein():
- g.impulse = True
- g.sign = -1
- scene.append_to_caption(' '*5)
- button(text='APPLY FORCE IN', bind=applyforcein)
- scene.append_to_caption(' '*5)
- g.general_text = wtext(text='####')
- scene.waitfor('textures')
- scene.visible = True
- def centrifugal_force(m, 𝛀, r):
- """Find the vector of the centrifugal force given
- mass (m), angular velocity vector (𝛀), and position vector (r).
- Fcf = m(𝛀 x r) x 𝛀 # Taylor p.344 """
- return m * cross(cross(𝛀, r), 𝛀)
- def coriolis_force(m, 𝛀, v):
- """Find the vector of the coriolis force given
- mass (m), angular velocity vector (𝛀), and velocity vector (v).
- Fcor = 2mṙ x 𝛀 = 2mv x 𝛀 # Taylor p.348 """
- return cross(2.0*m*v, 𝛀)
- class g: ...
- if __name__ == '__main__':
- sys.exit(main())
Add Comment
Please, Sign In to add comment