vpython_random_points_in_cylinder_moment_of_inertia

#31_08_21

'''

 CALCULATE THE MOMENT(S) OF INERTIA OF A CYLINDER ROTATING ABOUT CENTER

                        |
               +--------|--------+
               |        |        |
               +--------|--------+
                        |

 Ideas from:
 Rhett Allain's Physics Explained Channel
 "Finding the moment of inertia for a sphere by making 1000 random points"
 "Building a circle with random points - cartesian vs. polar coordinates"
 "Two Ways to Find the Moment of Inertia Tensor for a Solid Cube"

 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']

 In [1]: %timeit 2.0/3.0
 7.39 ns ± 0.0186 ns per loop

'''

from vpython import *
import sys

def main():
    #print(version, '\n')

    scene.width, scene.height = 1900, 980
    scene.background = color.black
    scene.align = 'left'
    scene.autoscale = False
    scene.autozoom = False
    scene.userspin = True
    scene.userzoom = True
    scene.userpan = False
    scene.resizable = False
    scene.range = 8.0

    '''
    For R=10:

    (x,0,0) -> 100/12 = 8.3333 # VERY THIN ROD ABOUT ANY AXIS
    (x,0,z) -> 103/12 = 8.5833 # ROD ABOUT Y-AXIS
    (x,y,z) -> 106/12 = 8.8333 # ROD ABOUT ANY AXIS

    '''

    SUM = 'xz' # 'x', 'xy', 'xz', 'xyz'

    print('\n', 'SUMMATION STRING: ', "'", SUM, "'\n", sep='')

    DRAW_CYLINDER = 0

    if DRAW_CYLINDER:
        sphere(pos=vec(0,0,0), radius=0.1, color=color.red)

    N = 5000 if DRAW_CYLINDER else 10000000

    LENGTH = 10.0
    LENGTH_2 = LENGTH / 2.0
    R = 1.0  # RADIUS
    M = 1.0  # TOTAL MASS
    m = M/N  # PART MASS

    # MOMENT OF INERTIA
    I = 0.0

    # CENTER OF MASS
    CM = vec(0.,0.,0.)

    for _ in range(N):

        while 1:
            r = vec(0.0, 2.0*random()-1.0, 2.0*random()-1.0)
            if r.mag < R: break
        r.x = random() * LENGTH - LENGTH_2

        if DRAW_CYLINDER:
            sphere(pos=r, radius=0.02, opacity=0.5)

        if SUM == 'x':              # R²/12
            # ∑ mx²
            r = vec(r.x, 0., 0.)

        elif SUM == 'xz':           # (R²+3)/12
            # ∑ m(x² + z²)
            r = vec(r.x, 0., r.z)

        elif SUM == 'xy':           # (R²+3)/12
            # ∑ m(x² + y²)
            r = vec(r.x, r.y, 0.)

        elif SUM == 'xyz':          # (R²+6)/12
            # ∑ m(x² + y² + z²)
            r = vec(r.x, r.y, r.z)

        else:
            print('INVALID SUMMATION STRING')
            print("valid strings are: 'x', 'xy', 'xz', 'xyz'")
            sys.exit()

        I += m * r.mag2

        CM += m * r

    CM /= M

    print('I = ', I, '\n', sep='')
    print('ML^2/12 = ', M * LENGTH**2 / 12.0, '\n')

    print('CM = ', CM)


class g: ...


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