%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%               File Name: Strain_CompatibilityAnalysis26mAASHTOTypeIVMidspanCSAS6_14rev17               Written by: Koral Eren               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;
clear all;
clc;

%Imput Parameters

Eps = 200000;
Ec = (3000*sqrt(40)+6900)*(2500/2300)^1.5;

Aps = 98.7*32;
Ag = 509031.24;
Ag_T = 528487.6385;

y_total = 1371.6;
y_total_composite = 1371.6 + 200;

d = y_total - 93.75;
d_composite = y_total_composite - 93.75;

yt = 743.3605;
yb = y_total - yt;

ec = 534.4895;

Ig = 1.0853*10^11;
Ig_T = 1.1391*10^11;

Ppe = 3464.485816857033 * 10^3;
Epsilon_ctop = 0.0035;
Epsilon_0 = 0.002;
fc_prime = 40;

%Initial Stage "Epsilon_si"

Epsilon_si = Ppe / (Aps * Eps);

%Initial Stage "Epsilon_ci"

ft = -(Ppe / Ag_T) + ( Ppe * ec * yt / Ig_T);
fb = -(Ppe / Ag_T) - ( Ppe * ec * yb / Ig_T);
Epsilon_t = ft / Eps;
Epsilon_b = fb / Eps;
Epsilon_ci_p = Epsilon_t + (Epsilon_b - Epsilon_t) * d / y_total;

Mg = 1053.82192461 * 10^6;
Ms = 1014 * 10^6;

Epsilon_t_DL = -(Mg + Ms) * yt / (Ig_T * Eps);
Epsilon_b_DL = (Mg + Ms) * yb / (Ig_T * Eps);
Epsilon_ci_DL = (ec / yb) * Epsilon_b_DL;

Epsilon_ci = abs(Epsilon_ci_p) + Epsilon_ci_DL;

%Calculation of area segments: Rectangle width 0.01 m or 1 cm%

Area = zeros(y_total_composite/0.01,1);

i=1;
for y=0.01:0.01:y_total_composite
    
    	if (y <= 200)
        	Area(i)=2500*0.01;
    	elseif (y <= 403.2)
        	Area(i)=508*0.01;
    	elseif (y <= 555.6)
        	Area(i)=(-2*y+1314.4508)*0.01;
    	elseif (y <= 1139.8)
        	Area(i)=203.2*0.01;
    	elseif (y <= 1368.4)
        	Area(i)=(2*y-2076.4)*0.01;
    	else
        	Area(i)=660.4*0.01;
    	end
    
    	i=i+1;
    
end

i=1;
j=1;

%Calculation part

for c = 0.01:0.01:y_total_composite
    
    	% Epsilon_s = Strain at prestressing steel at the level of the CGS.
    	Epsilon_s = Epsilon_si + Epsilon_ci + Epsilon_ctop * (d_composite-c)/c;

    	if (Epsilon_s <= 0.008)
        	fps = Eps * Epsilon_s;
    	elseif (Epsilon_s > 0.008)
        	fps = 1848 - 0.517 / (Epsilon_s - 0.005915);
    	end

    	if (fps > 1843.38)
        	fps = -10^15;
    	end

    	integration_limit = -(0.002-Epsilon_ctop)*c/Epsilon_ctop;

    	C = 0;
    	i = 1;

    	T = fps * Aps;

    	for y=0.01:0.01:c

        	Epsilon_c = Epsilon_ctop - Epsilon_ctop*y/c;

        	if (y <= integration_limit)
            		C = C + (-3333.3333333333*Epsilon_c+46.6666666667)*Area(i);
        	else
            		C = C + fc_prime*(2*(Epsilon_c/Epsilon_0)-(Epsilon_c/Epsilon_0)^2)*Area(i);
        	end

        	i = i+1;
        
    	end
    
    	if (abs(C-T) < 500)
        	break;
    	end
    
    	j = j+1;
    
end

Sum = 0;
i=1;
for y=0.01:0.01:c
    
    	Epsilon_c = Epsilon_ctop - Epsilon_ctop*y/c;

    	if (y <= integration_limit)
        	Cdy = (-3333.3333333333*Epsilon_c+46.6666666667)*Area(i);
    	else
        	Cdy = fc_prime*(2*(Epsilon_c/Epsilon_0)-(Epsilon_c/Epsilon_0)^2)*Area(i);
    	end

    	Sum = Sum + Cdy * y;
    	i = i+1;
    
end

y_bar = Sum / C;

z = d_composite - y_bar;

M = z * T;

M_kNm = M * 10^-6;

fprintf('Moment Resistance of the given section: %4.0f +- 0.5 kNm\n\n', M_kNm);