import sympy as sm
import sympy.physics.mechanics as me
me.init_vprinting()

m1, m2, m3, m4 = sm.symbols('m1, m2, m3, m4')

N = me.ReferenceFrame('N')

r1 = 3*N.x + 4*N.y
r2 = 3*N.y + 4*N.z
r3 = 4*N.x + 3*N.z
r4 = 4*N.y + 3*N.z

first_moment = m1*r1 + m2*r2 + m3*r3 + m4*r4

zeroth_moment = m1 + m2 + m3 + m4

mass_center = first_moment/zeroth_moment
mass_center

Ix = (
m1*me.cross(r1, me.cross(N.x, r1)) +
m2*me.cross(r2, me.cross(N.x, r2)) +
m3*me.cross(r3, me.cross(N.x, r3)) +
m4*me.cross(r4, me.cross(N.x, r4))
)
Ix

na = (N.x + N.y + N.z).normalize()

Ia = (
m1*me.cross(r1, me.cross(na, r1)) +
m2*me.cross(r2, me.cross(na, r2)) +
m3*me.cross(r3, me.cross(na, r3)) +
m4*me.cross(r4, me.cross(na, r4))
)
Ia

sm.simplify(Ia.dot(na))

Ia.dot(N.x)

r_r = -3*N.y + 4*N.x
r_l =  3*N.y + 4*N.x

m = sm.symbols('m')

Ix = (
    m*me.cross(r_r, me.cross(N.x, r_r)) +
    m*me.cross(r_l, me.cross(N.x, r_l))
)
Ix

m*me.cross(r_r, me.cross(N.x, r_r))

m*me.cross(r_l, me.cross(N.x, r_l))

Ix.dot(N.y)