//mass moment of inertia of 2D slice rotating about axis with tilt
ll=5;
q=1.732051;
x=ll/q;
py=linspace(-ll,ll,300);
px=x*ones(1,length(py));
pm=ones(1,length(px))/length(px);
points=[//X[mm];Y[mm];mass[g] of points to evaluate
    px;
    py;
    pm];
rps=5; //rotational speed [rev/s] about axis
w2=(rps*2*%pi)^2; //squared angular velocity
ra_theta=linspace(0,%pi); //rotation axes to evaluate [rad], 0 = X axis to the right
lr=length(ra_theta);

if(%f)
//method 1: add forces
ms=zeros(1,lr);
for i=1:lr
    t=ra_theta(i);
    rmat=[cos(t),sin(t);-sin(t),cos(t)];
    pr=rmat*points(1:2,:); //rotated points in frame where rotated axis is X+
    //centrifugal force F=m*w^2*r
    ctf=points(3,:).*w2.*pr(2,:); // [uN]
    mom=ctf.*pr(1,:); // [uN*mm] moment of inertia, +=CCW torque
    ms(i)=sum(mom); // [uN*mm] summed moment of inertia
end
plot(ra_theta,ms);
end

//method 2: add phasors
R=sqrt(points(1,:).^2+points(2,:).^2);
phi=atan(points(2,:),points(1,:));
phs=points(3,:).*(R.^2).*exp(%i*2*phi);
phss=sum(phs);
phs2=2*(x^3*q-(x^3*q^3)/3)/(2*ll);
pht=w2.*(0.5*exp(%i*-2*ra_theta));
md=imag(phss.*pht);
//plot(ra_theta,md,'r');
disp([phss;phs2]);
