mechanizmos58

vpython_animate_omega_Eulers_equations_1.py

Aug 23rd, 2021 (edited)
1,465
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # 23_08_21
  2.  
  3. '''
  4. tested on: Linux Mint 20.2, Geany 1.36, python 3.8.10,
  5. Web Epiphany, i5-7260U CPU @ 2.20GHz × 2, Iris Plus Graphics 640
  6. vpython version: ['7.6.1', 'jupyter']
  7.  
  8. '''
  9.  
  10. from vpython import *
  11. import sys
  12.  
  13. def main():
  14.  
  15.     scene.width, scene.height = 1000, 800
  16.     scene.align = 'left'
  17.     scene.resizable = False
  18.     scene.autoscale = False
  19.     scene.autozoom = False
  20.     scene.range = 1.5
  21.  
  22.     length = 1.0
  23.  
  24.     A = arrow(pos=vec(0,0,0), axis=vec(length,0,0), color=vec(0,0,1),
  25.                 shaftwidth=0.02, shininess=0.0)
  26.  
  27.  
  28.     def arrow_tip():
  29.         return (A.axis)
  30.     g.T = attach_trail(arrow_tip, color=color.red, trail_radius=0.002, retain=170)
  31.  
  32.  
  33.     def pause_run():
  34.         g.paused = not g.paused
  35.         pause.text='  RUN   ' if g.paused else 'PAUSE'
  36.         #g.T.clear()
  37.     scene.append_to_caption(' '*10)
  38.     pause = button(bind=pause_run)
  39.     g.paused = True
  40.     pause_run()
  41.  
  42.  
  43.     Omega = vec(1.8, 0.0, 0.06)
  44.     Omega_dot = vec(0, 0, 0)
  45.  
  46.     Ixx, Iyy, Izz = get_inertia()
  47.  
  48.     e1 = vec(length, 0, 0)
  49.     e2 = vec(0, length, 0)
  50.     e3 = vec(0, 0, length)
  51.  
  52.     dt = 0.02
  53.  
  54.     while 1:
  55.  
  56.         if g.paused: continue
  57.  
  58.         rate(100)
  59.  
  60.         Omega_dot.x = (Iyy-Izz) * Omega.y * Omega.z / Ixx
  61.         Omega_dot.y = (Izz-Ixx) * Omega.z * Omega.x / Iyy
  62.         Omega_dot.z = (Ixx-Iyy) * Omega.x * Omega.y / Izz
  63.  
  64.         Omega += Omega_dot * dt
  65.         T = Omega / 50.0
  66.         # ~ print(Omega)
  67.  
  68.         e1.value = rotate(e1, angle=T.z, axis=e3)
  69.         e1.value = rotate(e1, angle=T.y, axis=e2)
  70.         A.axis = e1
  71.  
  72.         e2.value = rotate(e2, angle=T.z, axis=e3)
  73.         e2.value = rotate(e2, angle=T.x, axis=e1)
  74.  
  75.         e3.value = rotate(e3, angle=T.x, axis=e1)
  76.         e3.value = rotate(e3, angle=T.y, axis=e2)
  77.  
  78.         e1 /= mag(e1)
  79.         e2 /= mag(e2)
  80.         e3 /= mag(e3)
  81.  
  82.         e1.value = cross(e2, e3)
  83.         e2.value = cross(e3, e1)
  84.         e3.value = cross(e1, e2)
  85.  
  86.         # ~ scene.waitfor("redraw")
  87.         # ~ scene.waitfor("draw_complete")
  88.  
  89.  
  90. def get_inertia():
  91.  
  92.     # Measured on Taylor's Classical Mechanics
  93.     h, w, d = 0.26, 0.18, 0.042
  94.     m = 1.58
  95.  
  96.     I1 = m/12.0 * (h*h + d*d)
  97.     I2 = m/12.0 * (w*w + d*d)
  98.     I3 = m/12.0 * (w*w + h*h)
  99.  
  100.     return I1, I2, I3
  101.  
  102.  
  103. class g: ...
  104.  
  105.  
  106. if __name__ == '__main__':
  107.     sys.exit(main())
  108.  
  109.  
RAW Paste Data