%Theoretical Optimum Shape Finding

%Seeks the shape optimised beam design with the minimum possible embodied
%carbon when shape is not restricted in any form

%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 EEEE 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

%finding the moments of the section and calling for curvature
global L0 d h bw 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;


%Big loop
beamno=1;

for Hmax=300:10:900

    %Analyse for this beam
    

    %Calculate flexure

    for s=0:5

        row=1;

        for h=dflange:5:Hmax

            %N1=2+(h-250)/10;

            for bw=200:5:700
                bf=bw+2240;
                Wu=67.67; %kN/m for ULS


                %Flexural performance

                d=Hmax-40;
                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

                Asmax=0.04*h*bw;
                if Asf>Asmax
                    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,s+1}=h;
                ASM{row,2,s+1}=bw;

                %Change the d now
                d=h-40;

                %Check whether it satisfies apploed moment
                Lx=L0*s/10;
                Mx=Wu/2*Lx*(L0-Lx);

                zx=Mx*1000000/(0.87*fs*Asf);

                xx=1.908*fs*Asf/(fc*bf);

                dx=zx+xx*0.4;

                if dx<zx/0.95
                    dx=zx/0.95;
                end

                if round(dx)>d
                    Asf=NaN;
                end

                ASM{row,3}=Asf;
                %ASM(N,i+1)=As;


                %Shear design

                Ved=Wu*L0/2-Wu*L0*s/10;
                %Desin shear stress
                Ved=Ved*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
                
                stress=Ved/(d*bw);

                ASM{row,4,s+1}=AwS;
                ASM{row,5,s+1}=shearcon;
                Ans=['for h= ', num2str(h),'and b= ', num2str(bw),' Aw/s avg = ', num2str(ASM{row,4,s+1}), ' with ', ASM{row,5,s+1}];
                disp (Ans);


                % Calculate EC

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

                Ans=['for h= ', num2str(h),'and b= ', num2str(bw),' EC = ', num2str(ASM{row,6,s+1})];
                disp (Ans);

                row=row+1;

            end

        end


        D=cell2mat(ASM(:,1,s+1));
        B=cell2mat(ASM(:,2,s+1));
        A=cell2mat(ASM(:,6,s+1));
        CV=cell2mat(ASM(:,8,s+1));

        [EM,ED]=min(A);
        DProf(s+1)=D(ED);
        BProf(s+1)=B(ED);
        EProf(s+1)=EM;
        CVProf(s+1)=CV(ED);
        
        %{
        dt=delaunayTriangulation(D,B);
        tri=dt.ConnectivityList;
        Di=dt.Points(:,1);
        Bi=dt.Points(:,2);
        F=scatteredInterpolant(D,B,A);
        Ai=F(Di,Bi);
        figure(s+1);
        trisurf(tri,Di,Bi,Ai);
        shading interp
        colormap('jet');
        xlabel('Depth(mm)');
        ylabel('Width(mm)');
        zlabel('EC (kgCO2)');
        colorbar;
        %caxis([30 40]);
        xlim([250 Hmax]);
        view(2);
        %}
        

    end
    
    FinalEC=sum(EProf)*2-EProf(1)-EProf(6);
    FinalCV=sum(CVProf)*2-CVProf(1)-CVProf(6);
    FDEH(beamno)=Hmax;
    FDEC(beamno)=FinalEC;
    EEEE(beamno,1)=Hmax;
    EEEE(beamno,2)=FinalEC;
    EEEE(beamno,3)=FinalCV;

    
    Length=[0 1 2 3 4 5];
    figure(7);
    plot(Length, DProf);
    hold on;
    ylim([0 Hmax+100]);
    figure(8);
    plot(Length, BProf);
    hold on;
    ylim([0 500]);
    
    
    beamno=beamno+1;
end

figure(9);
plot(FDEH,FDEC,'-r','LineWidth',2);
xlabel('Depth(mm)');
ylabel('EC(kgCO2)');