%Deflection of FOSF

%Feasible Optimum Shape Finding

%Seeks the shape optimised beam design with the minimum possible embodied
%carbon when depth profile is restricted to be parabolic

%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 'FinalMin' in the memory for design details of each beam
%Check 'Deflections' in the memory for deflections of each beam

global dflange bf fc Ec fs Es epcp;
dflange=160; %mm
%flange width=web +1680
%fc=30; %N/mm2
%Ec=33000; %N/mm2
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

F=[12 20 30 40 50];
E=[27000 30000 33000 35000 37000];

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


%Analyse for this beam
Hmid=600;
L=0:1:10;
LM=[0 5 10];

beamno=1;

%this is used in shear design
%sin2(teta)
sin2t=sind(2*acotd(2.5));

for Hmid=300:10:900

    row=1;

    for Hend=dflange:2:Hmid
        DM=[Hend Hmid Hend];
        p=polyfit(LM,DM,2);
        D=polyval(p,L);
        %plot(L,D);
        %hold on;

        %Calculate flexure



        for s=0:5

            bw=200;
            bf=bw+2240;
            Wu=67.67; % kN/m for ULS
            
            %Flexural performance

            d=Hmid-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

            

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

            %Change the H now
            h=polyval(p,L(s+1));
            d=h-40;
            H(row,s+1)=h;

            %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)>round(d)
                Asf=NaN;
            end

            Ans=['For h= ', num2str(h), 'at l= ', num2str(s+1), ' Asflex= ', num2str(Asf)];
            disp (Ans);


            
            %ASM(N,i+1)=As;


            %Shear design

            Ved=Wu*L0/2-Wu*L0*s/10;
            %Desin shear stress
            Ved=Ved*1000;

            %bshear=Ved/(2*d);
            bshear=(Ved)/(0.18*d*fc*(1-fc/250)*sin2t);
            bmin=200;
            basf=Asf*100/(4*h);
            bw=max([bshear bmin basf]);

            %B(row,s+1)=bw;
            B(s+1)=bw;
            
            Asmax=0.04*h*bw;
            if Asf>Asmax
                %Asf=NaN;
            end
            ASM{row,3}=Asf;

            %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
            %Ans=['for h= ', num2str(h),' Aw/s = ',num2str(AwS)];
            %disp (Ans);

            %Save l v AwS in a matrix
            %AwSM(x+1,1)=x*L0/10;
            %AwSM(x+1,2+(h-200)/10)=AwS;

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


            % Calculate EC

            if isnan(ASM{row,3,s+1}) | isnan(ASM{row,7,s+1})
                EC=NaN;
                StlV=NaN;
                ConV=NaN;
            else
                StlV=L0/10*((Asf/10^6)+ASM{row,7,s+1}*(h+bw-20*4)/10^6);
                %StlV=L0/10*(ASM{row,7,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,9,s+1}=EC;
            ASM{row,10,s+1}=StlV;
            ASM{row,11,s+1}=ConV;

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




            %{

            D=cell2mat(ASM(:,1,s+1));
            B=cell2mat(ASM(:,2,s+1));
            A=cell2mat(ASM(:,9,s+1));

            [EM,ED]=min(A);
            DProf(s+1)=D(ED);
            BProf(s+1)=B(ED);
            EEE(s+1)=EM;

            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 Hmid]);
            view(2);
            caxis([50 170]);

            %}


        end



        EEE=cell2mat(squeeze(ASM(row,9,:)));
        EEC=sum(EEE)*2-EEE(1)-EEE(6);
        CCC=cell2mat(squeeze(ASM(row,11,:)));
        CCV=sum(CCC)*2-CCC(1)-CCC(6);
        Ans=['for Hmax= ', num2str(Hmid),' EC = ', num2str(EEC)];
        disp (Ans);

        ECM(row)=EEC;
        END(row)=Hend;
        CVM(row)=CCV;
        AllBM(row,:)=B;
        AllDM(row,:)=D;

        row=row+1;

    end
    
    [MinEC, MinBeam]=min(ECM);
    MinEnd=END(MinBeam);
    FinalMin(beamno,1)=Hmid;
    FinalMin(beamno,2)=MinEnd;
    FinalMin(beamno,3)=MinEC;
    FinalMin(beamno,4)=CVM(MinBeam);
    
    %Now Save the shapes in the matrix
    
    %Depth Matrix
    BestDepths(beamno,:)=AllDM(MinBeam,:);
    %Width Matrix
    BMatrix=AllBM(MinBeam,:);
    RBMatrix=flip(BMatrix);
    FBMatrix=[BMatrix RBMatrix(2:end)];
    BestWidths(beamno,:)=FBMatrix;

    %Caclulate Deflections
    
    %Charactersitic Deflection

    W=48.22; %kN/m
    As=Asf;
    
    if isnan(Asf)
        vact=NaN;
    else
        for x=0:10
            lx=L0/10;
            d=BestDepths(beamno,x+1)-40;
            bw=BestWidths(beamno,x+1);
            m(x+1)=W*(L0*lx*x-lx*lx*x^2)/2; %kNm
            l(x+1)=x*lx; %m
            M=m(x+1);
            k(x+1)=Section2;
        end
        k(1)=0;
        k(11)=0;

        %interation of curvature

        k1(1)=0;
        for N=1:10
            li=l(1:N+1);
            ki=k(1:N+1);
            k1(N+1)=trapz(li,ki);
        end

        k2(1)=0;
        for N=1:10
            li=l(1:N+1);
            k1i=k1(1:N+1);
            k2(N+1)=trapz(li,k1i);
        end
        c=-k2(11)/L0;
        for N=0:10
            k2c(N+1)=k2(N+1)+c*L0*N/10;
        end

        vact=k2c(6)*-1;

    end
    
    %Frequent Deflection

    W=36.16; %kN/m for quasi-permanent loads
    As=Asf;
    
    if isnan(Asf)
        vfrq=NaN;
    else
        for x=0:10
            lx=L0/10;
            d=BestDepths(beamno,x+1)-40;
            bw=BestWidths(beamno,x+1);
            m(x+1)=W*(L0*lx*x-lx*lx*x^2)/2; %kNm
            l(x+1)=x*lx; %m
            M=m(x+1);
            k(x+1)=Section2;
        end
        k(1)=0;
        k(11)=0;

        %interation of curvature

        k1(1)=0;
        for N=1:10
            li=l(1:N+1);
            ki=k(1:N+1);
            k1(N+1)=trapz(li,ki);
        end

        k2(1)=0;
        for N=1:10
            li=l(1:N+1);
            k1i=k1(1:N+1);
            k2(N+1)=trapz(li,k1i);
        end
        c=-k2(11)/L0;
        for N=0:10
            k2c(N+1)=k2(N+1)+c*L0*N/10;
        end

        vfrq=k2c(6)*-1;

    end
    
    Deflections(beamno,1)=Hmid;
    Deflections(beamno,2)=vact*1000;
    Deflections(beamno,3)=vfrq*1000;
    
    %Plot shapes
    
    if isnan(MinEC)
        
    else
        %Length=[0 1 2 3 4 5];
        figure(1);
        HMatrix=AllDM(MinBeam,:);
        %RHMatrix=flip(HMatrix);
        %FHMatrix=[HMatrix RHMatrix(2:end)];
        %plot(L, FHMatrix,'-r', 'LineWidth', 2);
        plot(L, HMatrix,'-r', 'LineWidth', 2);
        hold on;
        ylim([0 Hmid+50]);
        figure(2);
        plot(L, FBMatrix/2,'-g', 'LineWidth', 2);
        hold on;
        plot(L, -FBMatrix/2,'-g', 'LineWidth', 2);
        hold on;
        ylim([-500 500]);
    end
    
    beamno=beamno+1;

end

HHH=FinalMin(:,1);
CCC=FinalMin(:,3);

figure(3);
plot(HHH,CCC,'-b', 'LineWidth', 2);
xlabel('Depth (mm)');
ylabel('EC (kgCO2)');