disp(' ');
disp(' mass_properties_GRB.m   ver 1.2  March 23, 2010 ');
disp(' ');
disp(' by Tom Irvine   Email: tomirvine@aol.com ');
disp(' ');
%
clear m;
clear GRB;
%
disp(' ');
disp(' Enter the units system ');
disp(' 1=English  2=metric ');
iu=input(' ');
%
mass_scale=1;
%
if(iu==1)
     disp(' Select input mass unit ');
     disp('  1=lbm  2=lbf sec^2/in  ');
     imu=input(' ');
     if(imu==1)
         mass_scale=386;
         mass_unit='lbm';       
     else
         mass_unit='lbf sec^2/in'; 
     end
     loc_unit='inch';
     inertia_unit='lbf sec^2 in'
else
     disp(' mass unit = kg ');
     mass_unit='kg';
     loc_unit='meters';
     inertia_unit='kg m^2';
end
%
disp(' ');
disp(' Select file input method ');
disp('   1=file preloaded into Matlab ');
disp('   2=Excel file ');
file_choice = input('');
%
disp(' ');
disp(' Mass Matrix ');
%
if(file_choice==1)
        m = input(' Enter the matrix name:  ');
end
if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        m = xlsread(xfile);
%         
end
%
m=m/mass_scale;
%
mass=m;
%
size(m);
%
num=max(size(m));
%
disp(' ');
disp(' The mass matrix is');
m
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
out1=sprintf(' The nodal coordinates must have three columns: X, Y & Z  (%s)',loc_unit);
disp(out1);
disp(' Select the nodal coordinate input method ');
%
disp(' ')
disp(' Select data input method ');
disp('   1=external ASCII file ');
disp('   2=file preloaded into Matlab ');
disp('   3=Excel file');
file_choice = input('');
%
if(file_choice==1)
        [filename, pathname] = uigetfile('*.*');
        filename = fullfile(pathname, filename);
        fid = fopen(filename,'r');
        THM = fscanf(fid,'%g %g %g',[3 inf]);
        fclose(fid);
        THM=THM';
end
if(file_choice==2)
        THM = input(' Enter the matrix name:  ');
end
if(file_choice==3)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        THM = xlsread(xfile);
%         
end
%
nodex=THM(:,1);
nodey=THM(:,2);
nodez=THM(:,3);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
nn=round(max(size(m))/6);
%
GRB=zeros(nn,6);
%
refx=0;
refy=0;
refz=0;
%
for(i=1:nn)    
%
    dx=-refx+nodex(i);
    dy=-refy+nodey(i);
    dz=-refz+nodez(i);
%
    GRB((i-1)*6+2,4)=-dz;
    GRB((i-1)*6+3,4)= dy;  
    GRB((i-1)*6+1,5)= dz;
    GRB((i-1)*6+3,5)=-dx;  
    GRB((i-1)*6+1,6)=-dy;
    GRB((i-1)*6+2,6)= dx; 
%
    GRB((i-1)*6+1,1)=1;
    GRB((i-1)*6+2,2)=1;  
    GRB((i-1)*6+3,3)=1;
    GRB((i-1)*6+4,4)=1;  
    GRB((i-1)*6+5,5)=1;
    GRB((i-1)*6+6,6)=1;  
end
GRB
%
disp(' ');
disp(' Rigid-body Mass Matrix ');
mq=GRB'*m*GRB
%
%
disp(' ');
out1=sprintf(' The reference point (%s): ',loc_unit);
disp(out1);
out1=sprintf('\n refx=%8.4g \n refy=%8.4g \n refz=%8.4g \n',refx,refy,refz);
disp(out1);
%
disp(' ');
out1=sprintf(' CG coordinates relative to reference point (%s) \n',loc_unit);
disp(out1);
Ixx=mq(4,4);
Ixy=mq(4,5);
Ixz=mq(4,6);
Iyy=mq(5,5);
Iyz=mq(5,6);
Izz=mq(6,6);
%
mx=mq(1,1);
my=mq(2,2);
mz=mq(3,3);
%
mmass=(mx+my+mz)/3;
%
xg= mq(2,6)/mmass;
yg=-mq(1,6)/mmass;
zg= mq(1,5)/mmass;
if(abs(xg)<1.0e-08)
    xg=0;
end
if(abs(yg)<1.0e-08)
    yg=0;
end
if(abs(zg)<1.0e-08)
    zg=0;
end
out1=sprintf(' xg = %8.4g  ',xg);
out2=sprintf(' yg = %8.4g  ',yg);
out3=sprintf(' zg = %8.4g  ',zg);
disp(out1);
disp(out2);
disp(out3);
%
disp(' ');
out1=sprintf(' Mass (%s)',mass_unit);
disp(out1);
disp(' ');
out1=sprintf(' mx = %8.4g  ',mx*mass_scale);
out2=sprintf(' my = %8.4g  ',my*mass_scale);
out3=sprintf(' mz = %8.4g  ',mz*mass_scale);
disp(out1);
disp(out2);
disp(out3);
%
Ixx_cg=Ixx-mmass*(yg^2+zg^2);
Iyy_cg=Iyy-mmass*(xg^2+zg^2);
Izz_cg=Izz-mmass*(xg^2+yg^2);
%
Ixy_cg=abs(Ixy)-abs(mmass*(xg*yg));
Ixz_cg=abs(Ixz)-abs(mmass*(xg*zg));
Iyz_cg=abs(Iyz)-abs(mmass*(yg*zg));
%
Icg=[Ixx_cg Ixy_cg Ixz_cg; Ixy_cg Iyy_cg Iyz_cg; Ixz_cg Iyz_cg Izz_cg];
%
III=[Ixx_cg,Iyy_cg,Izz_cg];
minI=min(abs(III));
minI=minI/10000;
for(i=1:3)
    for(j=1:3)
        if(i~=j)
            if(abs(Icg(i,j))<minI)
                Icg(i,j)=0;
            end
        end
    end
end    
%
disp(' ');
out5=sprintf(' Mass Moments of Inertia about CG (%s) ',inertia_unit);
disp(out5);
Icg