// Submitted by Dr. Triranjita Srivastava. Assistant Professor, Physics Dept., Kalindi College, University of Delhi

// Aim: Determination of the principal axes of moment of inertia through diagonalization 
// Example is a Dumbell with masses 'm1' and 'm2' situated at points, say coordinates are (1,1,0) and (-1,-1,0)

clear;
clc;
//*************************************************************
//Function for Kronecker Delta
//*************************************************************
function d=delta(i,j)
    if i==j then  
              d=1;
        else d=0;
    end
endfunction

//*************************************************************
// Input of number of particles at discrete points
//************************************************************* 
n=input('enter no. of particles  ')

r=zeros(3,3)
for i=1:n
    mprintf("Enter the mass (in kg) at point (%d):  ",i)
    M(i)=input("")
    mprintf("Enter the position (x,y,z) coordinate at point (%d):  ",i)
    for j =1:3
        r(i,j)=input("")
end
end

I=zeros(3,3)
for i=1:1:3
    for j=1:1:3
        for k=1:1:n
            I(i,j)=I(i,j)+(M(k)*(sum(r(k,:).^2)*delta(i,j)-(r(k,i).*r(k,j))))
        end
    end
end
disp("Moment of Inertia Tensor for given problem is:   ")
disp(I)     
[ab,x,bs]=bdiag(I);
disp("Moment of Inertia Tensor after diagonalization is:   ")
disp(ab)

