//Copyright (C) 2006  Mihai Pruna

//This program is free software; you can redistribute it and/or
//modify it under the terms of the GNU General Public License
//as published by the Free Software Foundation; either version 2
//of the License, or (at your option) any later version.

//This program is distributed in the hope that it will be useful,
//but WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//GNU General Public License for more details.
//http://www.gnu.org/licenses/gpl.html

//You should have received a copy of the GNU General Public License
//along with this program; if not, write to the Free Software
//Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

//Calculation of moments of inertia, including around the principal axes, of a cross section that can be constructed from rectangular components.


rectangles=evstr(input("Number of rectangular components?"));
x=zeros(10);
y=zeros(10);
b=zeros(10);
h=zeros(10);
ycg=zeros(10);
xcg=zeros(10);
area=zeros(10);
Ax=zeros(10);
Ay=zeros(10);
Ax2=zeros(10);
Ay2=zeros(10);
I0x=zeros(10);
I0y=zeros(10);
sectionxcg=0;
sectionycg=0;
sumarea=0;
sumay=0;
sumay2=0;
sumax=0;
sumax2=0;
sumI0x=0;
sumI0y=0;
sectionIx=0;
sectionIy=0;
sumaxy=0;
hole=zeros(10);
for i=1:rectangles
printf("for rectangle %d \n",i);
x(i)=evstr(input("Input x base coordinate of rectangle"));
y(i)=evstr(input("Input y base coordinate of rectangle"));
b(i)=evstr(input("Input length base of rectangle"));
h(i)=evstr(input("Input height of rectangle"));
hole(i)=evstr(input("Solid or hole (1 or 0)"));
if hole(i)==0
hole(i)=-1;
end

ycg(i)=y(i)+h(i)/2;
xcg(i)=x(i)+b(i)/2;
area(i)=hole(i)*b(i)*h(i);
Ax(i)=area(i)*xcg(i);
Ay(i)=area(i)*ycg(i);
Ax2(i)=Ax(i)*xcg(i);
Ay2(i)=Ay(i)*ycg(i);
I0x(i)=hole(i)*b(i)*(h(i)^3)/12;
I0y(i)=hole(i)*h(i)*(b(i)^3)/12;
sumarea=sumarea+area(i);
sumax=sumax+Ax(i);
sumay=sumay+Ay(i);
sumax2=sumax2+Ax2(i);
sumay2=sumay2+Ay2(i);
sumI0x=sumI0x+I0x(i);
sumI0y=sumI0y+I0y(i);
sumaxy=sumaxy+area(i)*xcg(i)*ycg(i);
end

for i=1:rectangles
printf("Element %d ",i);
printf(" base=%f",b(i));
printf(" height=%f",h(i));
printf(" x centroid=%f",xcg(i));
printf(" y centroid=%f",ycg(i));
printf(" area=%f",area(i));
printf(" Ax=%f",Ax(i));
printf(" Ay=%f",Ay(i));
printf(" Ax2=%f",Ax2(i));
printf(" Ay2=%f",Ay2(i));
printf(" I0x=%f",I0x(i));
printf(" I0y=%f",I0y(i));
printf("\n");
end

sectionxcg=sumax/sumarea;
sectionycg=sumay/sumarea;
sectionIx=sumI0x+sumay2-sectionycg*sumay;
sectionIy=sumI0y+sumax2-sectionxcg*sumax;
Ixy=sumaxy-sumarea*sectionxcg*sectionycg;
printf("Section x centroid=%f\n",sectionxcg);
printf("Section y centroid=%f\n",sectionycg);
printf("Section x moment of inertia around cg(Ix)=%f\n",sectionIx);
printf("Section y moment of inertia around cg(Iy)=%f\n",sectionIy);
printf("Section Ixy=%f\n",Ixy);
alpha=atan(2*Ixy/(sectionIy-sectionIx))/2;
Ixp=sectionIx*(cos(alpha))^2+sectionIy*(sin(alpha))^2-Ixy*sin(2*alpha);
Iyp=sectionIx*(sin(alpha))^2+sectionIy*(cos(alpha))^2+Ixy*sin(2*alpha);
alphadeg=alpha*180/3.1415926;
printf("Alpha(angle of principal x axis,degrees)=%f\n",alphadeg);
printf("Ixp=%f\n",Ixp);
printf("Iyp=%f\n",Iyp);