mechanizmos58

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:
  32. # v must be a vector to start with
  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,
  143.                         radius=0.07,
  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))
  171.     scene.append_to_caption(' radians/s\n')
  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 = 2mṙ 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