vpython_animate_omega_Eulers_equations_1.py

# 23_08_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
 vpython version: ['7.6.1', 'jupyter']

'''

from vpython import *
import sys

def main():

    scene.width, scene.height = 1000, 800
    scene.align = 'left'
    scene.resizable = False
    scene.autoscale = False
    scene.autozoom = False
    scene.range = 1.5

    length = 1.0

    A = arrow(pos=vec(0,0,0), axis=vec(length,0,0), color=vec(0,0,1),
                shaftwidth=0.02, shininess=0.0)


    def arrow_tip():
        return (A.axis)
    g.T = attach_trail(arrow_tip, color=color.red, trail_radius=0.002, retain=170)


    def pause_run():
        g.paused = not g.paused
        pause.text='  RUN   ' if g.paused else 'PAUSE'
        #g.T.clear()
    scene.append_to_caption(' '*10)
    pause = button(bind=pause_run)
    g.paused = True
    pause_run()


    Omega = vec(1.8, 0.0, 0.06)
    Omega_dot = vec(0, 0, 0)

    Ixx, Iyy, Izz = get_inertia()

    e1 = vec(length, 0, 0)
    e2 = vec(0, length, 0)
    e3 = vec(0, 0, length)

    dt = 0.02

    while 1:

        if g.paused: continue

        rate(100)

        Omega_dot.x = (Iyy-Izz) * Omega.y * Omega.z / Ixx
        Omega_dot.y = (Izz-Ixx) * Omega.z * Omega.x / Iyy
        Omega_dot.z = (Ixx-Iyy) * Omega.x * Omega.y / Izz

        Omega += Omega_dot * dt
        T = Omega / 50.0
        # ~ print(Omega)

        e1.value = rotate(e1, angle=T.z, axis=e3)
        e1.value = rotate(e1, angle=T.y, axis=e2)
        A.axis = e1

        e2.value = rotate(e2, angle=T.z, axis=e3)
        e2.value = rotate(e2, angle=T.x, axis=e1)

        e3.value = rotate(e3, angle=T.x, axis=e1)
        e3.value = rotate(e3, angle=T.y, axis=e2)

        e1 /= mag(e1)
        e2 /= mag(e2)
        e3 /= mag(e3)

        e1.value = cross(e2, e3)
        e2.value = cross(e3, e1)
        e3.value = cross(e1, e2)

        # ~ scene.waitfor("redraw")
        # ~ scene.waitfor("draw_complete")


def get_inertia():

    # Measured on Taylor's Classical Mechanics
    h, w, d = 0.26, 0.18, 0.042
    m = 1.58

    I1 = m/12.0 * (h*h + d*d)
    I2 = m/12.0 * (w*w + d*d)
    I3 = m/12.0 * (w*w + h*h)

    return I1, I2, I3


class g: ...


if __name__ == '__main__':
    sys.exit(main())