%Optimum Prismatic Shape Finding

%Seeks the prismatic beam design with the minimum possible embodied carbon

%Output figure shows the variation of embodied carbon with design depth

%Minimum point of the curve shows the optimum depth and minimum possibel
%embodied carbon

%Check ASM in the memory for design details of each beam

global dflange bf fc Ec fs Es;
dflange=160; %mm

fs=500; %N/mm2
Es=200000; %N/mm2
%epcp=(0.57*fc)/Ec;

%finding the moments of the section and calling for curvature
global L0 W d h bw As M;
L0=8; %m

%Concrete compressive strength and elastic modulus
F=[12 20 30 40 50];
E=[27000 30000 33000 35000 37000];

G=3;
fc=F(G);
Ec=E(G)/3;

row=1;

%this is used in shear design
%sin2(teta)
sin2t=sind(2*acotd(2.5));
    
for h=300:10:900
    
    %Minimum possible width
    
    d=h-40; %mm
    bw=200;
    bf=bw+2240;
    
    Wu=67.67; %kN/m ULS load

    %Flexural performance

    M=Wu*L0*L0/8;
    k=M*1000000/(bf*d*d*fc);
    z=d*(0.5+(0.25-(k/1.14))^0.5);
    if isreal(z)
        if z>0.95*d
            z=0.95*d;
        end
        Asf=M*1000000/(0.87*fs*z);
    else
        Asf=NaN;
    end

    x=(d-z)/0.4;
    Ans=['For h= ', num2str(h), 'and b= ', num2str(bw), ' Asflex= ', num2str(Asf)];
    disp (Ans);
    %N2=2+(h-bw)/10;
    ASM{row,1}=h;
    %ASM{row,2}=bw;
    
    %ASM(N,i+1)=As;

    
    %Shear design
    
    Vt=Wu*L0*1000/2;
    bshear=(Vt)/(0.18*d*fc*(1-fc/250)*sin2t); % Minimum width required for shear performance
    bmin=200; %Minimum possible width due to fire criteron or spacing of rebar
    bw=max([bshear bmin]);
    ASM{row,2}=bw;
    
    %Maximum reinforcement ratio
    Asmax=0.04*h*bw;
    if Asf>Asmax
        Asf=NaN;
    end
    ASM{row,3}=Asf;
    
    for x=1:11
        v(x)=abs(Wu*L0/2-x*L0*Wu/10);
        %Desin shear stress
        Ved=v(x)*1000;

        %Concrete Capacity Vrdc
        k1=1+(200/d)^0.5;
        k2=2;
        k=min(k1,k2);
        ro1=Asf/(bw*d);
        ro2=0.02;
        ro=min(ro1,ro2);
        Vrdc1=0.12*k*bw*d*(100*ro*fc)^(1/3);
        Vrdc2=0.035*bw*d*(fc^0.5)*(k^(3/2));
        Vrdc=max(Vrdc1,Vrdc2);

        Vmin=0.15*d*bw*(fc)^0.5;

        %Minimum reinforcement decision
        if Ved<Vrdc
            %disp ('No need for design reinforcement');
            %Minimum RF area
            %AwS=(0.08*bw*(fc)^0.5)/fs;
        end

        if Ved<Vmin
            shearcon='Minimum RF';
            %Minimum RF area
            AwS=(0.08*bw*(fc)^0.5)/fs;
            else 
                %disp ('Design reinforcement required');
                %Caclulation of design shear RF
                teta=0.5*asin((5.56*Ved)/(bw*d*fc*(1-fc/250)));
                CT=cot(teta);
                if CT<1
                    %disp ('Resize!!');
                    shearcon='Resize';
                    AwS=NaN;
                else if CT>=2.5
                    %disp ('Use cot(teta)=2.5');
                    AwS=Ved/(0.9*2.5*d*fs*0.87);
                    shearcon='Use cot(teta)=2.5';
                    %Calculate concrete strut strength
                    Vcon=0.13*d*bw*fc*(1-fc/250);
                    if Vcon<Ved
                        %disp ('Concrete strut fails');
                        shearcon='Concrete strut fails';
                        AwS=NaN;
                    end
                    else 
                        %disp ('Check RF at zcot(teta)');
                        shearcon='Check RF at zcot(teta)';
                        AwS=Ved/(0.9*CT*d*fs*0.87);
                    end
                end

        end
        AwST(x)=AwS;
    end
    ASM{row,6}=sum(AwST(2:end))/10;
    ASM{row,7}=shearcon;
    Ans=['for h= ', num2str(h),'and b= ', num2str(bw),' Aw/s avg = ', num2str(ASM{row,6}), ' with ', ASM{row,7}];
    disp (Ans);


    % Calculate EC

    if isnan(ASM{row,3}) | isnan(ASM{row,6})
        EC=NaN;
        ConV=NaN;
    else
        StlV=L0*((Asf/10^6)+ASM{row,6}*(h+bw-20*4)/10^6);
        ConV=(L0*(h*bw)/10^6)-StlV;
        EC=StlV*7850*1.2+ConV*2400*0.132;
    end
    ASM{row,8}=EC;
    ASM{row,9}=ConV;

    row=row+1;

end
    
D1=cell2mat(ASM(:,1));
A1=cell2mat(ASM(:,8));
V1=cell2mat(ASM(:,4));
VA1=cell2mat(ASM(:,5));
q=find(~isnan(A1),1);
D=D1(q:end);
A=A1(q:end);
V=V1(q:end);
VA=VA1(q:end);

xlow=300;
xhigh=900;
figure(1);
plot(D,A,'-r','LineWidth',2);
xlabel('Depth(mm)');
ylabel('Embodied Carbon (kgCO2)');
xlim([xlow xhigh]);