%% read geometry information
geo_data=xlsread('geo_info');
bc = geo_data(1); % width of the section
hc = geo_data(2); % height of the section
cover = geo_data(3:6); % thickness of the covering layer
A_margin = geo_data(7:10); % area of rebars located at the margin of the section
n_margin = geo_data(11:14); % number of rebars located at the margin of the section
A_corner = geo_data(15:18); % area of rebars located at the corner of the section
global fc;
global fy;
fc = geo_data(19); % peak stress of concrete
fy = geo_data(20); % yielding stress of rebar
global yps_cu;
global Es;
yps_cu = 0.0033;
Es = 200000;
%% initiate the coordinates and areas of each rebar
% the order is [corner(lt,rt,lb,rb), margin_t(l2r),margin_b(l2r),margin_l(t2b),margin_r(t2b)]
coor_rebar = ones(2,4+sum(n_margin));
Area_rebar = ones(1,4+sum(n_margin));
Area_rebar(1:4)=A_corner;
Area_rebar(5:4+n_margin(1))=A_margin(1)*ones(1,n_margin(1));
Area_rebar(5+n_margin(1):4+n_margin(1)+n_margin(2)) = A_margin(2)*ones(1,n_margin(2));
Area_rebar(5+n_margin(1)+n_margin(2):4+sum(n_margin)-n_margin(4))=A_margin(3)*ones(1,n_margin(3));
Area_rebar(5+sum(n_margin)-n_margin(4):4+sum(n_margin))=A_margin(4)*ones(1,n_margin(4));
coor_rebar(:,1:4) = [-bc/2+cover(3),bc/2-cover(4),-bc/2+cover(3),bc/2-cover(4);
hc/2-cover(1),hc/2-cover(1),-hc/2+cover(2),-hc/2+cover(2)];
assistant = linspace(coor_rebar(1,1),coor_rebar(1,2),n_margin(1)+2);
coor_rebar(:,5:4+n_margin(1)) = [assistant(2:n_margin(1)+1);coor_rebar(2,1)*ones(1,n_margin(1))];
assistant = linspace(coor_rebar(1,3),coor_rebar(1,4),n_margin(2)+2);
coor_rebar(:,5+n_margin(1):4+n_margin(1)+n_margin(2)) = [assistant(2:n_margin(2)+1);coor_rebar(2,3)*ones(1,n_margin(2))];
assistant = linspace(coor_rebar(2,1),coor_rebar(2,3),n_margin(3)+2);
coor_rebar(:,5+n_margin(1)+n_margin(2):4+sum(n_margin)-n_margin(4))=[coor_rebar(1,1)*ones(1,n_margin(3));
assistant(2:n_margin(3)+1)];
assistant = linspace(coor_rebar(2,2),coor_rebar(2,4),n_margin(4)+2);
coor_rebar(:,5+sum(n_margin)-n_margin(4):4+sum(n_margin))=[coor_rebar(1,2)*ones(1,n_margin(4));
assistant(2:n_margin(4)+1)];
clear assistant;
% the sequence of [lb, rb, lt, rt]
%% cut the concrete section into small parts and initate the coordinates
n_width = 100; % number of parts in the width direction
n_height = 100; % number of parts in the height direction
top = [-bc/2;hc/2]; % coordinate of the top point of the section
coor_con = zeros(2,n_width*n_height);
coor_con(1,:) = kron(ones(1,n_height),linspace(bc/2/n_width-bc/2,bc/2-bc/2/n_width,n_width));
coor_con(2,:) = kron(linspace(hc/2/n_height-hc/2,hc/2-hc/2/n_height,n_height),ones(1,n_width));
A_part = bc*hc/(n_width*n_height);
%% canculate the ultimate strength of the section, assuming that compression point is always on the positive y-axis
n_theta = 100; % number of slices of theta
n_yps = 100; % number of slices of ypsilon
Capicity = zeros(2,n_theta*n_yps);
Theta = linspace(1e-8,pi/2,n_theta);
Yps_rebar = linspace(-8*fy/Es,yps_cu,n_yps);
i=1;
global yps_rebar;
while i <= n_theta
disp(i);
theta = Theta(i);
T = Tran(theta);
rebar_cur = T * coor_rebar;
con_cur = T * coor_con;
top_cur = T * top;
j=1;
while j <= n_yps
disp(j);
yps_rebar = Yps_rebar(j);
flaga = sum(ssr_con(yps(con_cur(2,:)/top_cur(2))));
Nu = sum(ssr_con(yps(con_cur(2,:)/top_cur(2))))*A_part ...
+ ssr_rebar(yps(rebar_cur(2,:)/top_cur(2)))*Area_rebar';
Mu = ssr_con(yps(con_cur(2,:)/top_cur(2)))*con_cur(2,:)'*A_part ...
+(ssr_rebar(yps(rebar_cur(2,:)/top_cur(2))).*rebar_cur(2,:))*Area_rebar';
% if y_neu>max(con_cur(2,1),con_cur(2,4)) % if the neutral axis lies above point 1&4
% Nu = integral(@(y)ssr_con(yps(y/con_cur(2,3))).*(y-con_cur(2,3))*((con_cur(1,3)-con_cur(1,4))/(con_cur(2,3)-con_cur(2,4))-(con_cur(1,3)-con_cur(1,1))/(con_cur(2,3)-con_cur(2,1))),y_neu,con_cur(2,3)) ...
% + ssr_rebar(yps(rebar_cur(2,:)/con_cur(2,3)))*Area_rebar';
% Mu = integral(@(y)y.*ssr_con(yps(y/con_cur(2,3))).*(y-con_cur(2,3))*((con_cur(1,3)-con_cur(1,4))/(con_cur(2,3)-con_cur(2,4))-(con_cur(1,3)-con_cur(1,1))/(con_cur(2,3)-con_cur(2,1))),y_neu,con_cur(2,3)) ...
% + (ssr_rebar(yps(rebar_cur(2,:)/con_cur(2,3))).*rebar_cur(2,:))*Area_rebar';
%
% elseif min(con_cur(2,1),con_cur(2,4))<=y_neu && y_neu<=max(con_cur(2,1),con_cur(2,4)) % if the neutral axis lies between point 1&4
% Nu = integral(@(y)ssr_con(yps(y/con_cur(2,3))).*(y-con_cur(2,3))*((con_cur(1,3)-con_cur(1,4))/(con_cur(2,3)-con_cur(2,4))-(con_cur(1,3)-con_cur(1,1))/(con_cur(2,3)-con_cur(2,1))),max(con_cur(2,1),con_cur(2,4)),con_cur(2,3)) ...
% + integral(@(y)ssr_con(yps(y/con_cur(2,3))).*min(bc/cos(theta),hc/sin(theta)),y_neu,max(con_cur(2,1),con_cur(2,4))) ...
% + ssr_rebar(yps(rebar_cur(2,:)/con_cur(2,3)))*Area_rebar';
% Mu = integral(@(y)y.*ssr_con(yps(y/con_cur(2,3))).*(y-con_cur(2,3))*((con_cur(1,3)-con_cur(1,4))/(con_cur(2,3)-con_cur(2,4))-(con_cur(1,3)-con_cur(1,1))/(con_cur(2,3)-con_cur(2,1))),max(con_cur(2,1),con_cur(2,4)),con_cur(2,3)) ...
% + integral(@(y)y.*ssr_con(yps(y/con_cur(2,3))).*min(bc/cos(theta),hc/sin(theta)),y_neu,max(con_cur(2,1),con_cur(2,4))) ...
% + (ssr_rebar(yps(rebar_cur(2,:)/con_cur(2,3))).*rebar_cur(2,:))*Area_rebar';
%
% elseif y_neu<min(con_cur(2,1),con_cur(2,4)) && y_neu>= con_cur(2,2)
% Nu = integral(@(y)ssr_con(yps(y/con_cur(2,3))).*(y-con_cur(2,3))*((con_cur(1,3)-con_cur(1,4))/(con_cur(2,3)-con_cur(2,4))-(con_cur(1,3)-con_cur(1,1))/(con_cur(2,3)-con_cur(2,1))),max(con_cur(2,1),con_cur(2,4)),con_cur(2,3)) ...
% + integral(@(y)ssr_con(yps(y/con_cur(2,3))).*min(bc/cos(theta),hc/sin(theta)),min(con_cur(2,1),con_cur(2,4)),max(con_cur(2,1),con_cur(2,4))) ...
% + integral(@(y)ssr_con(yps(y/con_cur(2,3))).*(y-con_cur(2,2))*((con_cur(1,2)-con_cur(1,4))/(con_cur(2,2)-con_cur(2,4))-(con_cur(1,2)-con_cur(1,1))/(con_cur(2,2)-con_cur(2,1))),y_neu,min(con_cur(2,1),con_cur(2,4))) ...
% + ssr_rebar(yps(rebar_cur(2,:)/con_cur(2,3)))*Area_rebar';
% Mu = integral(@(y)y.*ssr_con(yps(y/con_cur(2,3))).*(y-con_cur(2,3))*((con_cur(1,3)-con_cur(1,4))/(con_cur(2,3)-con_cur(2,4))-(con_cur(1,3)-con_cur(1,1))/(con_cur(2,3)-con_cur(2,1))),max(con_cur(2,1),con_cur(2,4)),con_cur(2,3)) ...
% + integral(@(y)y.*ssr_con(yps(y/con_cur(2,3))).*min(bc/cos(theta),hc/sin(theta)),min(con_cur(2,1),con_cur(2,4)),max(con_cur(2,1),con_cur(2,4))) ...
% + integral(@(y)y.*ssr_con(yps(y/con_cur(2,3))).*(y-con_cur(2,2))*((con_cur(1,2)-con_cur(1,4))/(con_cur(2,2)-con_cur(2,4))-(con_cur(1,2)-con_cur(1,1))/(con_cur(2,2)-con_cur(2,1))),y_neu,min(con_cur(2,1),con_cur(2,4))) ...
% + (ssr_rebar(yps(rebar_cur(2,:)/con_cur(2,3))).*rebar_cur(2,:))*Area_rebar';
%
% else
% Nu = integral(@(y)ssr_con(yps(y/con_cur(2,3))).*(y-con_cur(2,3))*((con_cur(1,3)-con_cur(1,4))/(con_cur(2,3)-con_cur(2,4))-(con_cur(1,3)-con_cur(1,1))/(con_cur(2,3)-con_cur(2,1))),max(con_cur(2,1),con_cur(2,4)),con_cur(2,3)) ...
% + integral(@(y)ssr_con(yps(y/con_cur(2,3)))*min(bc/cos(theta),hc/sin(theta)),min(con_cur(2,1),con_cur(2,4)),max(con_cur(2,1),con_cur(2,4)),'ArrayValued',true) ...
% + integral(@(y)ssr_con(yps(y/con_cur(2,3))).*(y-con_cur(2,2))*((con_cur(1,2)-con_cur(1,4))/(con_cur(2,2)-con_cur(2,4))-(con_cur(1,2)-con_cur(1,1))/(con_cur(2,2)-con_cur(2,1))),con_cur(2,2),min(con_cur(2,1),con_cur(2,4))) ...
% + ssr_rebar(yps(rebar_cur(2,:)/con_cur(2,3)))*Area_rebar';
% Mu = integral(@(y)y.*ssr_con(yps(y/con_cur(2,3))).*(y-con_cur(2,3))*((con_cur(1,3)-con_cur(1,4))/(con_cur(2,3)-con_cur(2,4))-(c