% Design of Flat Slab Reinforcement

% 3x3 bay flat slab design
% Grid line 1,2,3,4 in X & A,B,C,D in Y

% Design Spans (m), X is the longer span
%XSpan=5; 
%YSpan=5;
CL=3.5; %Column height
% Slab Depth (mm)
%h=300;

% Column size (mm)
ColX=400;
ColY=400;

%Slab Depth range
hmin=200;
hmax=400;

%Column Spacing range
SpanMin=4;
SpanMax=10;

% Material Properties
fck=30;
%fck=20;
%fck=40;
fsy=500;
fckcol=30;
%fckcol=40;

%Carbon Coefficients kgCO2/kg
ecSlabCon=0.132;
%ecSlabCon=0.121;
%ecSlabCon=0.159;
%ecColumnCon=0.159;
ecColumnCon=0.132;
ecSteel=1.20;

%Provide more RF for deflection and shear
RFincrease=1.00;

% Loads (kN/m2)

%SWC=25; % Concrete slef weight kN/m3
%Finishes=1.8; % kN/m2
%Services=0.5;
%SImposed=2.5;
%Partitions=1.0;
%Cladding=0.0; %kN/m

%For Concrete Centre Loads
SWC=25; % Concrete slef weight kN/m3
Finishes=1.5; % kN/m2
Services=0.0;
SImposed=2.5;
Partitions=0.0;
Cladding=10.0; %kN/m

%loop for a range of h

ECM{1,1}='Slab Depth (mm)';
ECM{1,2}='Concrete EC (kgCO2)';
ECM{1,3}='Flex Steel EC (kgCO2)';
ECM{1,4}='Shear Steel EC (kgCO2)';
ECM{1,5}='Column EC (kgCO2)';
ECM{1,6}='Total EC (kgCO2)';
ECM{1,7}='Under RF Condition';
ECM{1,8}='Max RF Condition';
ECM{1,9}='Shear column perimiter Condition';
ECM{1,10}='Shear control perimiter condition';
ECM{1,11}='Deflection Condition';
ECM{1,12}='Max B RF';
ECM{1,13}='Max T RF';
ECM{1,14}='Average B RF';
ECM{1,15}='Average T RF';
ECM{1,16}='Min B RF';
ECM{1,17}='Min T RF';

%Column Design Chart

%Assume d2/h=0.15 on average
Mch=0:0.05:0.40;
Nch=0:0.1:1.3;
%      0  0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40       
Ach=[0.00 0.13 0.30 0.46 0.63 0.80 0.95 1.00 1.00 ;  %0
     0.00 0.03 0.20 0.36 0.53 0.70 0.87 1.00 1.00 ;  %0.1
     0.00 0.00 0.12 0.28 0.45 0.62 0.78 0.96 1.00 ;  %0.2
     0.00 0.00 0.12 0.31 0.50 0.68 0.86 1.00 1.00 ;  %0.3
     0.00 0.00 0.19 0.39 0.59 0.77 0.96 1.00 1.00 ;  %0.4
     0.00 0.05 0.27 0.47 0.57 0.86 1.00 1.00 1.00 ;  %0.5
     0.05 0.20 0.36 0.57 0.77 0.96 1.00 1.00 1.00 ;  %0.6
     0.20 0.31 0.49 0.69 0.89 1.00 1.00 1.00 1.00 ;  %0.7
     0.35 0.45 0.61 0.80 1.00 1.00 1.00 1.00 1.00 ;  %0.8
     0.51 0.57 0.74 0.91 1.00 1.00 1.00 1.00 1.00 ;  %0.9
     0.67 0.70 0.87 1.00 1.00 1.00 1.00 1.00 1.00 ;  %1.0
     0.83 0.85 1.00 1.00 1.00 1.00 1.00 1.00 1.00 ;  %1.1
     0.98 0.98 1.00 1.00 1.00 1.00 1.00 1.00 1.00 ;  %1.2
     1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 ]; %1.3

 [MMCH,NNCH]=meshgrid(Mch,Nch);

 %Matrixing Concrete Grades
 FC=[20 30 40];
 ECCON=[0.121 0.132 0.159];
 
 AllMinECsMat=[];
 AllMinECDepthsMat=[];
 
 for grade=2:2
 
 ecSlabCon=ECCON(grade);
 fck=FC(grade);
 
 AllECMat=[];
 
 AllCompRFMat=[];
 AllColumnMat=[];
 AllDefMat=[];
 AllShearFMat=[];
 AllShearCMat=[];
 
 SpanN=1;
 
for SP=7:0.25:7
    
    XSpan=SP; 
    YSpan=SP;

    ECSlab=[];

    AllRFDetails={};

    beam=1;

    for h=hmin:5:hmax

        beam=beam+1;

        dro0=((fck)^0.5)/1000;

        RFDesign={}; ShearDesign={};

        status=1;
        condition1='Viable';
        condition2='Viable';
        condition3='Viable';
        condition4='Viable';
        
        StatCompRF=1;
        StatColumn=1;
        StatDef=1;
        StatShearF=1;
        StatShearC=0;

        SelfWeight=SWC*h/1000;

        Dead = SelfWeight + Finishes + Services;
        Imposed = SImposed + Partitions;

        ULSLoad = 1.25*Dead + 1.5*Imposed;
        SLSQPLoad =Dead + 0.3*Imposed;
        fac0=310/(500*SLSQPLoad/(1.15*ULSLoad))*1.03*RFincrease;
        fac=min(fac0,1.5);
        

        cover=30;

        % ANALYSE & DESIGN OF GRID LINE 2(=3)

            StripWidth=YSpan; % m

            EffSpan = (XSpan*1000 - ColX + h)/1000; %m

            MEndSpan = (1.25*Dead*0.090 + 1.5*Imposed*0.10)*StripWidth*EffSpan*EffSpan; %kNm

            MIntSup = ULSLoad*0.094*StripWidth*EffSpan*EffSpan; %kNm

            MIntSpan = (1.25*Dead*0.066 + 1.5*Imposed*0.086)*StripWidth*EffSpan*EffSpan; %kNm

            AspectR=XSpan/YSpan;

            k1=75-(AspectR-1)*(75-62)/(2-1);
            k2=25+(AspectR-1)*(35-25)/(2-1);

            WColStrip2=StripWidth/2;
            WMidStrip2=StripWidth/2;

            %Column Stip Moments

                MColIntSup2=(k1/100)*MIntSup/WColStrip2; % kNm/m
                MColEndSpan2=0.55*MEndSpan/WColStrip2; %kNm/m
                MColIntSpan2=0.55*MIntSpan/WColStrip2; %kNm/m

            %Middle Stip Moments

                MMidIntSup2=(k2/100)*MIntSup/WMidStrip2; % kNm/m
                MMidEndSpan2=0.45*MEndSpan/WMidStrip2; %kNm/m
                MMidIntSpan2=0.45*MIntSpan/WMidStrip2; %kNm/m

            dx=h-cover-20/2;

            % Minimum area of Reinforcement

            fctm=0.3*(fck^0.666);    
            Asmin=0.26*fctm*1000*dx/fsy;
            AsminC=Asmin;
            Asmax=0.04*1000*h;


            %Column Strip Renforcement

                %Column Strip End Span

                M=MColEndSpan2;
                K=(M*10^6)/(1000*dx*dx*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at ColEndSpan2'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dx*(1+(1-3.53*K)^0.5);
                zmax=0.95*dx;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                if As>Asmax
                    disp('Maximum Reinforcement Reached');
                    status=0;
                    condition2='Maximum Reinforcement Reached';
                end
                As=max(As,Asmin);
                AsColEndSpan2=As;

                dro=As/(h*1000);
                
                %fac=min(facM*1.03,1.5);

                if dro<=dro0
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) + 3.2*(fck^0.5)*(((dro0/dro)-1)^1.5))*fac;
                else
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) )*fac;
                end

                LDM(1)=LD;

                %Column Strip Internal Span

                M=MColIntSpan2;
                K=(M*10^6)/(1000*dx*dx*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at ColIntSpan2'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dx*(1+(1-3.53*K)^0.5);
                zmax=0.95*dx;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                if As>Asmax
                    disp('Maximum Reinforcement Reached');
                    status=0;
                    condition2='Maximum Reinforcement Reached';
                end
                As=max(As,Asmin);
                AsColIntSpan2=As;

                dro=As/(h*1000);
                
                %fac=min(facM*1.09,1.5);

                if dro<=dro0
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) + 3.2*(fck^0.5)*(((dro0/dro)-1)^1.5))*fac;
                else
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) )*fac;
                end

                LDM(2)=LD;

                %Column Strip Internal Support

                M=MColIntSup2;
                K=(M*10^6)/(1000*dx*dx*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at ColIntSup2'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dx*(1+(1-3.53*K)^0.5);
                zmax=0.95*dx;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                if As>Asmax
                    disp('Maximum Reinforcement Reached');
                    status=0;
                    condition2='Maximum Reinforcement Reached';
                end
                As=max(As,Asmin);
                AsColIntSup2=As;

            %Middle Strip Reinforcement

                %Middle Strip End Span

                M=MMidEndSpan2;
                K=(M*10^6)/(1000*dx*dx*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at MidEndSpan2'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dx*(1+(1-3.53*K)^0.5);
                zmax=0.95*dx;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                if As>Asmax
                    disp('Maximum Reinforcement Reached');
                    status=0;
                    condition2='Maximum Reinforcement Reached';
                end
                As=max(As,Asmin);
                AsMidEndSpan2=As;

                dro=As/(h*1000);
                %fac=min(facM*1.03,1.5);

                if dro<=dro0
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) + 3.2*(fck^0.5)*(((dro0/dro)-1)^1.5))*fac;
                else
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) )*fac;
                end

                LDM(3)=LD;

                %Middle Strip Internal Span

                M=MMidIntSpan2;
                K=(M*10^6)/(1000*dx*dx*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at MidIntSpan2'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dx*(1+(1-3.53*K)^0.5);
                zmax=0.95*dx;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                if As>Asmax
                    disp('Maximum Reinforcement Reached');
                    status=0;
                    condition2='Maximum Reinforcement Reached';
                end
                As=max(As,Asmin);
                AsMidIntSpan2=As;

                dro=As/(h*1000);
                %fac=min(facM*1.09,1.5);

                if dro<=dro0
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) + 3.2*(fck^0.5)*(((dro0/dro)-1)^1.5))*fac;
                else
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) )*fac;
                end

                LDM(4)=LD;

                %Middle Strip Internal Support

                M=MMidIntSup2;
                K=(M*10^6)/(1000*dx*dx*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at MidIntSup2'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dx*(1+(1-3.53*K)^0.5);
                zmax=0.95*dx;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                if As>Asmax
                    disp('Maximum Reinforcement Reached');
                    status=0;
                    condition2='Maximum Reinforcement Reached';
                end
                As=max(As,Asmin);
                AsMidIntSup2=As;

            % Column internal Strip Adjustment

                TotPanelAs=AsColIntSup2*WColStrip2 + AsMidIntSup2*WMidStrip2;
                ColMidAs=TotPanelAs/2;
                WColMid2=0.125*StripWidth*2;
                AsColMidSup2=ColMidAs/WColMid2;

                ColOutAs=WColStrip2*AsColIntSup2-ColMidAs;
                WColOut2=WColStrip2-WColMid2;
                AsColOutSup2=ColOutAs/WColOut2;

            % At 2A column edge: perpendicular to Edge

                be=ColY+ColX;
                Mt=0.17*be*dx*dx*fck/(10^6);
                Mout=(k1/100)*0.04*ULSLoad*StripWidth*EffSpan*EffSpan/WColStrip2;
                bz=1000;
                if Mout>Mt
                    disp('Moment redistribution at corners required');
                    Mout=Mt;
                    bz=be;
                end
                K=(Mout*10^6)/(bz*dx*dx*fck);
                z0=0.5*dx*(1+(1-3.53*K)^0.5);
                zmax=0.95*dx;
                z=min(z0,zmax);
                As=(Mout*10^6)/(0.87*fsy*z);
                AsColEndSup2=As;
                WColEndSup2=ColX+2*ColY;

            % At edges generally: 25% of span

                AsMidEndSup2=AsColEndSpan2*0.25;



        % ANALYSE & DESIGN OF GRID LINE B(=C)

           StripWidth=XSpan; % m

            EffSpan = (YSpan*1000 - ColY + h)/1000; %m

            MEndSpan = (1.25*Dead*0.090 + 1.5*Imposed*0.10)*StripWidth*EffSpan*EffSpan; %kNm

            MIntSup = ULSLoad*0.094*StripWidth*EffSpan*EffSpan; %kNm

            MIntSpan = (1.25*Dead*0.066 + 1.5*Imposed*0.086)*StripWidth*EffSpan*EffSpan; %kNm

            AspectR=XSpan/YSpan;

            k1=75;
            k2=25;

            WColStripB=YSpan/2;
            WMidStripB=StripWidth-WColStripB;

            %Column Stip Moments

                MColIntSupB=(k1/100)*MIntSup/WColStripB; % kNm/m
                MColEndSpanB=0.55*MEndSpan/WColStripB; %kNm/m
                MColIntSpanB=0.55*MIntSpan/WColStripB; %kNm/m

            %Middle Stip Moments

                MMidIntSupB=(k2/100)*MIntSup/WMidStripB; % kNm/m
                MMidEndSpanB=0.45*MEndSpan/WMidStripB; %kNm/m
                MMidIntSpanB=0.45*MIntSpan/WMidStripB; %kNm/m

            dy=h-cover-20-20/2;

            % Minimum area of Reinforcement

            fctm=0.3*(fck^0.666);    
            Asmin=0.26*fctm*1000*dy/fsy;


            %Column Strip Renforcement

                %Column Strip End Span

                M=MColEndSpanB;
                K=(M*10^6)/(1000*dy*dy*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at ColEndSpanB'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dy*(1+(1-3.53*K)^0.5);
                zmax=0.95*dy;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                if As>Asmax
                    disp('Maximum Reinforcement Reached');
                    status=0;
                    condition2='Maximum Reinforcement Reached';
                end
                As=max(As,Asmin);
                AsColEndSpanB=As;

                dro=As/(h*1000);
                %fac=min(facM*1.03,1.5);

                if dro<=dro0
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) + 3.2*(fck^0.5)*(((dro0/dro)-1)^1.5))*fac;
                else
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) )*fac;
                end

                LDM(5)=LD;

                %Column Strip Internal Span

                M=MColIntSpanB;
                K=(M*10^6)/(1000*dy*dy*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at ColIntSpanB'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dy*(1+(1-3.53*K)^0.5);
                zmax=0.95*dy;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                if As>Asmax
                    disp('Maximum Reinforcement Reached');
                    status=0;
                    condition2='Maximum Reinforcement Reached';
                end
                As=max(As,Asmin);
                AsColIntSpanB=As;

                dro=As/(h*1000);
                %fac=min(facM*1.09,1.5);

                if dro<=dro0
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) + 3.2*(fck^0.5)*(((dro0/dro)-1)^1.5))*fac;
                else
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) )*fac;
                end

                LDM(6)=LD;

                %Column Strip Internal Support

                M=MColIntSupB;
                K=(M*10^6)/(1000*dy*dy*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at ColIntSupB'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dy*(1+(1-3.53*K)^0.5);
                zmax=0.95*dy;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                if As>Asmax
                    disp('Maximum Reinforcement Reached');
                    status=0;
                    condition2='Maximum Reinforcement Reached';
                end
                As=max(As,Asmin);
                AsColIntSupB=As;

            %Middle Strip Reinforcement

                %Middle Strip End Span

                M=MMidEndSpanB;
                K=(M*10^6)/(1000*dy*dy*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at MidEndSpanB'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dy*(1+(1-3.53*K)^0.5);
                zmax=0.95*dy;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                if As>Asmax
                    disp('Maximum Reinforcement Reached');
                    status=0;
                    condition2='Maximum Reinforcement Reached';
                end
                As=max(As,Asmin);
                AsMidEndSpanB=As;

                dro=As/(h*1000);
                %fac=min(facM*1.03,1.5);

                if dro<=dro0
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) + 3.2*(fck^0.5)*(((dro0/dro)-1)^1.5))*fac;
                else
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) )*fac;
                end

                LDM(7)=LD;

                %Middle Strip Internal Span

                M=MMidIntSpanB;
                K=(M*10^6)/(1000*dy*dy*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at MidIntSpanB'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dy*(1+(1-3.53*K)^0.5);
                zmax=0.95*dy;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                if As>Asmax
                    disp('Maximum Reinforcement Reached');
                    status=0;
                    condition2='Maximum Reinforcement Reached';
                end
                As=max(As,Asmin);
                AsMidIntSpanB=As;

                dro=As/(h*1000);
                %fac=min(facM*1.09,1.5);

                if dro<=dro0
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) + 3.2*(fck^0.5)*(((dro0/dro)-1)^1.5))*fac;
                else
                    LD=1.2*( 11 + 1.5*(fck^0.5)*(dro0/dro) )*fac;
                end

                LDM(8)=LD;

                %Middle Strip Internal Support

                M=MMidIntSupB;
                K=(M*10^6)/(1000*dy*dy*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at MidIntSupB'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dy*(1+(1-3.53*K)^0.5);
                zmax=0.95*dy;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                if As>Asmax
                    disp('Maximum Reinforcement Reached');
                    status=0;
                    condition2='Maximum Reinforcement Reached';
                end
                As=max(As,Asmin);
                AsMidIntSupB=As;

            % Column internal Strip Adjustment

                TotPanelAs=AsColIntSupB*WColStripB + AsMidIntSupB*WMidStripB;
                ColMidAs=TotPanelAs/2;
                WColMidB=0.125*StripWidth*2;
                AsColMidSupB=ColMidAs/WColMidB;

                ColOutAs=WColStripB*AsColIntSupB-ColMidAs;
                WColOutB=WColStripB-WColMidB;
                AsColOutSupB=ColOutAs/WColOutB;

            % At 2A column edge: perpendicular to Edge

                be=ColY+ColX;
                Mt=0.17*be*dy*dy*fck/(10^6);
                Mout=(k1/100)*0.04*ULSLoad*StripWidth*EffSpan*EffSpan/WColStripB;
                bz=1000;
                if Mout>Mt
                    disp('Moment redistribution at corners required');
                    Mout=Mt;
                    bz=be;
                end
                K=(Mout*10^6)/(bz*dy*dy*fck);
                z0=0.5*dy*(1+(1-3.53*K)^0.5);
                zmax=0.95*dy;
                z=min(z0,zmax);
                As=(Mout*10^6)/(0.87*fsy*z);
                AsColEndSupB=As;
                WColEndSupB=ColY+2*ColX;

            % At edges generally: 25% of span

                AsMidEndSupB=AsColEndSpanB*0.25;


        % ANALYSE & DESIGN OF GRID LINE 1(=4)


            StripWidth = YSpan/2; % m

            EffSpan = (XSpan*1000 - ColX + h)/1000; %m

            Mareaload = ULSLoad*0.086*StripWidth*EffSpan*EffSpan; %kNm
            Mcladding = Cladding*0.086*EffSpan*EffSpan; %kNm

            AspectR=XSpan/YSpan;

            k1=75-(AspectR-1)*(75-62)/(2-1);
            k2=25+(AspectR-1)*(35-25)/(2-1);

            WColStrip1=YSpan/4;
            WMidStrip1=StripWidth - WColStrip1;

            %Column Stip Moments

                MColIntSup1=(k1*Mareaload/100+Mcladding)/WColStrip1; % kNm/m
                MColEndSpan1=(0.55*Mareaload+Mcladding)/WColStrip1; %kNm/m
                MColIntSpan1=(0.55*Mareaload+Mcladding)/WColStrip1; %kNm/m

            %Middle Stip Moments

                MMidIntSup1=(k2/100)*Mareaload/WMidStrip1; % kNm/m
                MMidEndSpan1=0.45*Mareaload/WMidStrip1; %kNm/m
                MMidIntSpan1=0.45*Mareaload/WMidStrip1; %kNm/m

            dx=h-cover-20/2;

            % Minimum area of Reinforcement

            fctm=0.3*(fck^0.666);    
            Asmin=0.26*fctm*1000*dx/fsy;


            %Column Strip Renforcement

                %Column Strip End Span

                M=MColEndSpan1;
                K=(M*10^6)/(1000*dx*dx*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at ColEndSpan1'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dx*(1+(1-3.53*K)^0.5);
                zmax=0.95*dx;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                As=max(As,Asmin);
                AsColEndSpan1=As;

                %Column Strip Internal Span

                M=MColIntSpan1;
                K=(M*10^6)/(1000*dx*dx*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at ColIntSpan1'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dx*(1+(1-3.53*K)^0.5);
                zmax=0.95*dx;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                As=max(As,Asmin);
                AsColIntSpan1=As;

                %Column Strip Internal Support

                M=MColIntSup1;
                K=(M*10^6)/(1000*dx*dx*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at ColIntSup1'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dx*(1+(1-3.53*K)^0.5);
                zmax=0.95*dx;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                As=max(As,Asmin);
                AsColIntSup1=As;

            %Middle Strip Reinforcement

                %Middle Strip End Span

                M=MMidEndSpan1;
                K=(M*10^6)/(1000*dx*dx*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at MidEndSpan1'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dx*(1+(1-3.53*K)^0.5);
                zmax=0.95*dx;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                As=max(As,Asmin);
                AsMidEndSpan1=As;

                %Middle Strip Internal Span

                M=MMidIntSpan1;
                K=(M*10^6)/(1000*dx*dx*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at MidIntSpan1'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dx*(1+(1-3.53*K)^0.5);
                zmax=0.95*dx;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                As=max(As,Asmin);
                AsMidIntSpan1=As;

                %Middle Strip Internal Support

                M=MMidIntSup1;
                K=(M*10^6)/(1000*dx*dx*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at MidIntSup1'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dx*(1+(1-3.53*K)^0.5);
                zmax=0.95*dx;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                As=max(As,Asmin);
                AsMidIntSup1=As;

            % Column internal Strip Adjustment

                TotPanelAs=AsColIntSup1*WColStrip1 + AsMidIntSup1*WMidStrip1;
                ColMidAs=TotPanelAs/2;
                WColMid1=0.125*StripWidth*2;
                AsColMidSup1=ColMidAs/WColMid1;

                ColOutAs=WColStrip1*AsColIntSup1-ColMidAs;
                WColOut1=WColStrip1-WColMid1;
                AsColOutSup1=ColOutAs/WColOut1;

            % At 1A column edge: perpendicular to Edge

                be=ColX/2+ColY;
                Mt=0.17*be*dx*dx*fck/(10^6);
                Mout=(k1/100)*0.04*ULSLoad*StripWidth*EffSpan*EffSpan/WColStrip1;
                bz=1000;
                if Mout>Mt
                    disp('Moment redistribution at corners required');
                    Mout=Mt;
                    bz=be;
                end
                K=(Mout*10^6)/(bz*dx*dx*fck);
                z0=0.5*dx*(1+(1-3.53*K)^0.5);
                zmax=0.95*dx;
                z=min(z0,zmax);
                As=(Mout*10^6)/(0.87*fsy*z);
                AsColEndSup1=As;
                WColEndSup1=ColY+ColX;


        % ANALYSE & DESIGN OF GRID LINE A(=D)

            StripWidth = XSpan/2; % m

            EffSpan = (YSpan*1000 - ColY + h)/1000; %m

            Mareaload = ULSLoad*0.086*StripWidth*EffSpan*EffSpan; %kNm
            Mcladding = Cladding*0.086*EffSpan*EffSpan; %kNm

            AspectR=XSpan/YSpan;

            k1=75;
            k2=25;

            WColStripA = YSpan/4;
            WMidStripA = StripWidth - WColStripA;

            %Column Stip Moments

                MColIntSupA=(k1*Mareaload/100+Mcladding)/WColStripA; % kNm/m
                MColEndSpanA=(0.55*Mareaload+Mcladding)/WColStripA; %kNm/m
                MColIntSpanA=(0.55*Mareaload+Mcladding)/WColStripA; %kNm/m

            %Middle Stip Moments

                MMidIntSupA=(k2/100)*Mareaload/WMidStripA; % kNm/m
                MMidEndSpanA=0.45*Mareaload/WMidStripA; %kNm/m
                MMidIntSpanA=0.45*Mareaload/WMidStripA; %kNm/m

            dy=h-cover-20-20/2;

            % Minimum area of Reinforcement

            fctm=0.3*(fck^0.666);    
            Asmin=0.26*fctm*1000*dy/fsy;


            %Column Strip Renforcement

                %Column Strip End Span

                M=MColEndSpanA;
                K=(M*10^6)/(1000*dy*dy*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at ColEndSpanA'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dy*(1+(1-3.53*K)^0.5);
                zmax=0.95*dy;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                As=max(As,Asmin);
                AsColEndSpanA=As;

                %Column Strip Internal Span

                M=MColIntSpanA;
                K=(M*10^6)/(1000*dy*dy*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at ColIntSpanA'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dy*(1+(1-3.53*K)^0.5);
                zmax=0.95*dy;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                As=max(As,Asmin);
                AsColIntSpanA=As;

                %Column Strip Internal Support

                M=MColIntSupA;
                K=(M*10^6)/(1000*dy*dy*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at ColIntSupA'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dy*(1+(1-3.53*K)^0.5);
                zmax=0.95*dy;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                As=max(As,Asmin);
                AsColIntSupA=As;

            %Middle Strip Reinforcement

                %Middle Strip End Span

                M=MMidEndSpanA;
                K=(M*10^6)/(1000*dy*dy*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at MidEndSpanA'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dy*(1+(1-3.53*K)^0.5);
                zmax=0.95*dy;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                As=max(As,Asmin);
                AsMidEndSpanA=As;

                %Middle Strip Internal Span

                M=MMidIntSpanA;
                K=(M*10^6)/(1000*dy*dy*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at MidIntSpanA'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dy*(1+(1-3.53*K)^0.5);
                zmax=0.95*dy;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                As=max(As,Asmin);
                AsMidIntSpanA=As;

                %Middle Strip Internal Support

                M=MMidIntSupA;
                K=(M*10^6)/(1000*dy*dy*fck);
                if K > 0.168
                    rfk=['Compression Reinforcement Reqired since K= ' num2str(K) ' at MidIntSupA'];
                    disp(rfk);
                    status=0;
                    StatCompRF=0;
                    condition1='Compression Reinforcement Reqired';
                end
                z0=0.5*dy*(1+(1-3.53*K)^0.5);
                zmax=0.95*dy;
                z=min(z0,zmax);
                As=(M*10^6)/(0.87*fsy*z);
                As=max(As,Asmin);
                AsMidIntSupA=As;

            % Column internal Strip Adjustment

                TotPanelAs=AsColIntSupA*WColStripA + AsMidIntSupA*WMidStripA;
                ColMidAs=TotPanelAs/2;
                WColMidA=0.125*StripWidth*2;
                AsColMidSupA=ColMidAs/WColMidA;

                ColOutAs=WColStripA*AsColIntSupA-ColMidAs/2;
                WColOutA=WColStripA-WColMidA;
                AsColOutSupA=ColOutAs/WColOutA;

            % At 1A column edge: perpendicular to Edge

                be=ColY/2+ColX;
                Mt=0.17*be*dy*dy*fck/(10^6);
                Mout=(k1/100)*0.04*ULSLoad*StripWidth*EffSpan*EffSpan/WColStripA;
                bz=1000;
                if Mout>Mt
                    disp('Moment reistribution at corners required');
                    Mout=Mt;
                    bz=be;
                end
                K=(Mout*10^6)/(bz*dy*dy*fck);
                z0=0.5*dy*(1+(1-3.53*K)^0.5);
                zmax=0.95*dy;
                z=min(z0,zmax);
                As=(Mout*10^6)/(0.87*fsy*z);
                AsColEndSupA=As;
                WColEndSupA=ColY+ColX;


        % TABULATE RF DISTRIBUTION

        RFDesign{1,1}='Layer';
        RFDesign{1,2}='Position';
        RFDesign{1,3}='As';
        RFDesign{1,4}='Width';
        RFDesign{1,5}='Length';
        RFDesign{1,6}='Count';
        RFDesign{1,7}='Slab Area';
        RFDesign{1,8}='Steel m3';
        RFDesign{1,9}='T12 spacing';
        RFDesign{1,10}='T20 spacing';


        % Bottom 1st layer (B1)

        RFDesign{2,1}='B1';
        RFDesign{2,2}='ColEndSpan1';
        RFDesign{2,3}=AsColEndSpan1*RFincrease;
        RFDesign{2,4}=WColStrip1;
        RFDesign{2,5}=XSpan;
        RFDesign{2,6}=4;
        RFDesign{2,7}=RFDesign{2,4}*RFDesign{2,5}*RFDesign{2,6};
        RFDesign{2,8}=RFDesign{2,7}*RFDesign{2,3}/10^6;

        RFDesign{3,1}='B1';
        RFDesign{3,2}='ColIntSpan1';
        RFDesign{3,3}=AsColIntSpan1*RFincrease;
        RFDesign{3,4}=WColStrip1;
        RFDesign{3,5}=XSpan;
        RFDesign{3,6}=2;
        RFDesign{3,7}=RFDesign{3,4}*RFDesign{3,5}*RFDesign{3,6};
        RFDesign{3,8}=RFDesign{3,7}*RFDesign{3,3}/10^6;

        RFDesign{4,1}='B1';
        RFDesign{4,2}='MidEndSpan1';
        RFDesign{4,3}=AsMidEndSpan1*RFincrease;
        RFDesign{4,4}=WMidStrip1;
        RFDesign{4,5}=XSpan;
        RFDesign{4,6}=4;
        RFDesign{4,7}=RFDesign{4,4}*RFDesign{4,5}*RFDesign{4,6};
        RFDesign{4,8}=RFDesign{4,7}*RFDesign{4,3}/10^6;

        RFDesign{5,1}='B1';
        RFDesign{5,2}='MidIntSpan1';
        RFDesign{5,3}=AsMidIntSpan1*RFincrease;
        RFDesign{5,4}=WMidStrip1;
        RFDesign{5,5}=XSpan;
        RFDesign{5,6}=2;
        RFDesign{5,7}=RFDesign{5,4}*RFDesign{5,5}*RFDesign{5,6};
        RFDesign{5,8}=RFDesign{5,7}*RFDesign{5,3}/10^6;

        RFDesign{6,1}='B1';
        RFDesign{6,2}='ColEndSpan2';
        RFDesign{6,3}=AsColEndSpan2*RFincrease;
        RFDesign{6,4}=WColStrip2;
        RFDesign{6,5}=XSpan;
        RFDesign{6,6}=4;
        RFDesign{6,7}=RFDesign{6,4}*RFDesign{6,5}*RFDesign{6,6};
        RFDesign{6,8}=RFDesign{6,7}*RFDesign{6,3}/10^6;

        RFDesign{7,1}='B1';
        RFDesign{7,2}='ColIntSpan2';
        RFDesign{7,3}=AsColIntSpan2*RFincrease;
        RFDesign{7,4}=WColStrip2;
        RFDesign{7,5}=XSpan;
        RFDesign{7,6}=2;
        RFDesign{7,7}=RFDesign{7,4}*RFDesign{7,5}*RFDesign{7,6};
        RFDesign{7,8}=RFDesign{7,7}*RFDesign{7,3}/10^6;

        RFDesign{8,1}='B1';
        RFDesign{8,2}='MidEndSpan2';
        RFDesign{8,3}=AsMidEndSpan2*RFincrease;
        RFDesign{8,4}=WMidStrip2;
        RFDesign{8,5}=XSpan;
        RFDesign{8,6}=4;
        RFDesign{8,7}=RFDesign{8,4}*RFDesign{8,5}*RFDesign{8,6};
        RFDesign{8,8}=RFDesign{8,7}*RFDesign{8,3}/10^6;

        RFDesign{9,1}='B1';
        RFDesign{9,2}='MidIntSpan2';
        RFDesign{9,3}=AsMidIntSpan2*RFincrease;
        RFDesign{9,4}=WMidStrip2;
        RFDesign{9,5}=XSpan;
        RFDesign{9,6}=2;
        RFDesign{9,7}=RFDesign{9,4}*RFDesign{9,5}*RFDesign{9,6};
        RFDesign{9,8}=RFDesign{9,7}*RFDesign{9,3}/10^6;

        B1FloorArea=sum(cell2mat(RFDesign(2:9,7)));
        B1SteelVolume=sum(cell2mat(RFDesign(2:9,8)));

        %RF Map
        if status==1

            xx=[];yy=[];xxx=[];yyy=[];xxxx=[];yyyy=[];

            xx=[XSpan XSpan XSpan];
            y1(1)=RFDesign{2,4};
            y1(2)=RFDesign{4,4};
            y1(3)=RFDesign{8,4}/2;
            y1(4)=RFDesign{6,4};
            y1(5)=RFDesign{8,4}/2;
            yy=[y1 y1(3:5) flip(y1(1:2))];

            limxx=length(xx);
            xxx(1)=0;
            for i=1:limxx
                xxx(i+1)=xx(i)+xxx(i);
            end
            limxxx=length(xxx);
            xxxx(1)=0;
            for i=1:limxxx-2
                xxxx(i*2)=xxx(i+1)-0.001;
                xxxx(i*2+1)=xxx(i+1)+0.001;
            end
            xxxx(length(xxxx)+1)=xxx(end);

            limyy=length(yy);
            yyy(1)=0;
            for i=1:limyy
                yyy(i+1)=yy(i)+yyy(i);
            end
            limyyy=length(yyy);
            yyyy(1)=0;
            for i=1:limyyy-2
                yyyy(i*2)=yyy(i+1)-0.001;
                yyyy(i*2+1)=yyy(i+1)+0.001;
            end
            yyyy(length(yyyy)+1)=yyy(end);

            [X,Y]=meshgrid(xxxx,yyyy);

            ZB1=zeros(length(yyyy),length(xxxx));
            ZB1(1:2,1:2)=RFDesign{2,3};
            ZB1(1:2,3:4)=RFDesign{3,3};
            ZB1(3:4,1:2)=RFDesign{4,3};
            ZB1(3:4,3:4)=RFDesign{5,3};
            ZB1(5:6,1:2)=RFDesign{8,3};
            ZB1(5:6,3:4)=RFDesign{9,3};
            ZB1(7:8,1:2)=RFDesign{6,3};
            ZB1(7:8,3:4)=RFDesign{7,3};
            ZB1(9:10,1:4)=ZB1(5:6,1:4);
            ZB1(11:16,1:4)=ZB1(5:10,1:4);
            ZB1(17:20,1:4)=flip(ZB1(1:4,1:4),1);
            ZB1(1:20,5:6)=flip(ZB1(1:20,1:2),2);


            figure(1);
            surf(X,Y,ZB1);
            colormap 'parula';
            colorbar;
            xlim([0 XSpan*3]);
            ylim([0 YSpan*3]);
            %zlim([0 1500]);
            xticks(xxx);
            yticks([0 YSpan YSpan*2 YSpan*3]);
            view(2);

        end

            % Bottom 2nd layer (B2)

            RFDesign{10,1}='B2';
            RFDesign{10,2}='ColEndSpanA';
            RFDesign{10,3}=AsColEndSpanA*RFincrease;
            RFDesign{10,4}=WColStripA;
            RFDesign{10,5}=YSpan;
            RFDesign{10,6}=4;
            RFDesign{10,7}=RFDesign{10,4}*RFDesign{10,5}*RFDesign{10,6};
            RFDesign{10,8}=RFDesign{10,7}*RFDesign{10,3}/10^6;

            RFDesign{11,1}='B2';
            RFDesign{11,2}='ColIntSpanA';
            RFDesign{11,3}=AsColIntSpanA*RFincrease;
            RFDesign{11,4}=WColStripA;
            RFDesign{11,5}=YSpan;
            RFDesign{11,6}=2;
            RFDesign{11,7}=RFDesign{11,4}*RFDesign{11,5}*RFDesign{11,6};
            RFDesign{11,8}=RFDesign{11,7}*RFDesign{11,3}/10^6;

            RFDesign{12,1}='B2';
            RFDesign{12,2}='MidEndSpanA';
            RFDesign{12,3}=AsMidEndSpanA*RFincrease;
            RFDesign{12,4}=WMidStripA;
            RFDesign{12,5}=YSpan;
            RFDesign{12,6}=4;
            RFDesign{12,7}=RFDesign{12,4}*RFDesign{12,5}*RFDesign{12,6};
            RFDesign{12,8}=RFDesign{12,7}*RFDesign{12,3}/10^6;

            RFDesign{13,1}='B2';
            RFDesign{13,2}='MidIntSpanA';
            RFDesign{13,3}=AsMidIntSpanA*RFincrease;
            RFDesign{13,4}=WMidStripA;
            RFDesign{13,5}=YSpan;
            RFDesign{13,6}=2;
            RFDesign{13,7}=RFDesign{13,4}*RFDesign{13,5}*RFDesign{13,6};
            RFDesign{13,8}=RFDesign{13,7}*RFDesign{13,3}/10^6;

            RFDesign{14,1}='B2';
            RFDesign{14,2}='ColEndSpanB';
            RFDesign{14,3}=AsColEndSpanB*RFincrease;
            RFDesign{14,4}=WColStripB;
            RFDesign{14,5}=YSpan;
            RFDesign{14,6}=4;
            RFDesign{14,7}=RFDesign{14,4}*RFDesign{14,5}*RFDesign{14,6};
            RFDesign{14,8}=RFDesign{14,7}*RFDesign{14,3}/10^6;

            RFDesign{15,1}='B2';
            RFDesign{15,2}='ColIntSpanB';
            RFDesign{15,3}=AsColIntSpanB*RFincrease;
            RFDesign{15,4}=WColStripB;
            RFDesign{15,5}=YSpan;
            RFDesign{15,6}=2;
            RFDesign{15,7}=RFDesign{15,4}*RFDesign{15,5}*RFDesign{15,6};
            RFDesign{15,8}=RFDesign{15,7}*RFDesign{15,3}/10^6;

            RFDesign{16,1}='B2';
            RFDesign{16,2}='MidEndSpanB';
            RFDesign{16,3}=AsMidEndSpanB*RFincrease;
            RFDesign{16,4}=WMidStripB;
            RFDesign{16,5}=YSpan;
            RFDesign{16,6}=4;
            RFDesign{16,7}=RFDesign{16,4}*RFDesign{16,5}*RFDesign{16,6};
            RFDesign{16,8}=RFDesign{16,7}*RFDesign{16,3}/10^6;

            RFDesign{17,1}='B2';
            RFDesign{17,2}='MidIntSpanB';
            RFDesign{17,3}=AsMidIntSpanB*RFincrease;
            RFDesign{17,4}=WMidStripB;
            RFDesign{17,5}=YSpan;
            RFDesign{17,6}=2;
            RFDesign{17,7}=RFDesign{17,4}*RFDesign{17,5}*RFDesign{17,6};
            RFDesign{17,8}=RFDesign{17,7}*RFDesign{17,3}/10^6;

            B2FloorArea=sum(cell2mat(RFDesign(10:17,7)));
            B2SteelVolume=sum(cell2mat(RFDesign(10:17,8)));

            %RF Map
            if status==1
                xx=[];yy=[];xxx=[];yyy=[];xxxx=[];yyyy=[];

                yy=[YSpan YSpan YSpan];
                x1(1)=RFDesign{10,4};
                x1(2)=RFDesign{12,4};
                x1(3)=RFDesign{16,4}/2;
                x1(4)=RFDesign{14,4};
                x1(5)=RFDesign{16,4}/2;
                xx=[x1 x1(3:5) flip(x1(1:2))];

                limyy=length(yy);
                yyy(1)=0;
                for i=1:limyy
                    yyy(i+1)=yy(i)+yyy(i);
                end
                limyyy=length(yyy);
                yyyy(1)=0;
                for i=1:limyyy-2
                    yyyy(i*2)=yyy(i+1)-0.001;
                    yyyy(i*2+1)=yyy(i+1)+0.001;
                end
                yyyy(length(yyyy)+1)=yyy(end);

                limxx=length(xx);
                xxx(1)=0;
                for i=1:limxx
                    xxx(i+1)=xx(i)+xxx(i);
                end
                limxxx=length(xxx);
                xxxx(1)=0;
                for i=1:limxxx-2
                    xxxx(i*2)=xxx(i+1)-0.001;
                    xxxx(i*2+1)=xxx(i+1)+0.001;
                end
                xxxx(length(xxxx)+1)=xxx(end);

                [X,Y]=meshgrid(xxxx,yyyy);

                ZB2=zeros(length(yyyy),length(xxxx));
                ZB2(1:2,1:2)=RFDesign{10,3};
                ZB2(3:4,1:2)=RFDesign{11,3};
                ZB2(1:2,3:4)=RFDesign{12,3};
                ZB2(3:4,3:4)=RFDesign{13,3};
                ZB2(1:2,5:6)=RFDesign{16,3};
                ZB2(3:4,5:6)=RFDesign{17,3};
                ZB2(1:2,7:8)=RFDesign{14,3};
                ZB2(3:4,7:8)=RFDesign{15,3};
                ZB2(1:4,9:10)=ZB2(1:4,5:6);
                ZB2(1:4,11:16)=ZB2(1:4,5:10);
                ZB2(1:4,17:20)=flip(ZB2(1:4,1:4),2);
                ZB2(5:6,1:20)=flip(ZB2(1:2,1:20),1);

                figure(2);
                surf(X,Y,ZB2);
                colormap 'parula';
                colorbar;
                xlim([0 XSpan*3]);
                ylim([0 YSpan*3]);
                %zlim([0 1500]);
                xticks([0 XSpan XSpan*2 XSpan*3]);
                yticks(yyy);
                view(2);
            end


        % Top 1st layer (T1)

        RFDesign{18,1}='T1';
        RFDesign{18,2}='ColEndSup1';
        RFDesign{18,3}=AsColEndSup1*RFincrease;
        RFDesign{18,4}=WColStrip1;
        RFDesign{18,5}=XSpan/2;
        RFDesign{18,6}=4;
        RFDesign{18,7}=RFDesign{18,4}*RFDesign{18,5}*RFDesign{18,6};
        RFDesign{18,8}=RFDesign{18,7}*RFDesign{18,3}/10^6;

        RFDesign{19,1}='T1';
        RFDesign{19,2}='MidEndSup2';
        RFDesign{19,3}=AsMidEndSup2*RFincrease;
        RFDesign{19,4}=WMidStrip1+WMidStrip2/2;
        RFDesign{19,5}=XSpan/2;
        RFDesign{19,6}=6;
        RFDesign{19,7}=RFDesign{19,4}*RFDesign{19,5}*RFDesign{19,6};
        RFDesign{19,8}=RFDesign{19,7}*RFDesign{19,3}/10^6;

        RFDesign{20,1}='T1';
        RFDesign{20,2}='ColEndSup2';
        RFDesign{20,3}=AsColEndSup2*RFincrease;
        RFDesign{20,4}=WColStrip2;
        RFDesign{20,5}=XSpan/2;
        RFDesign{20,6}=4;
        RFDesign{20,7}=RFDesign{20,4}*RFDesign{20,5}*RFDesign{20,6};
        RFDesign{20,8}=RFDesign{20,7}*RFDesign{20,3}/10^6;

        RFDesign{21,1}='T1';
        RFDesign{21,2}='ColMidSup1';
        RFDesign{21,3}=AsColMidSup1*RFincrease;
        RFDesign{21,4}=WColMid1;
        RFDesign{21,5}=XSpan;
        RFDesign{21,6}=4;
        RFDesign{21,7}=RFDesign{21,4}*RFDesign{21,5}*RFDesign{21,6};
        RFDesign{21,8}=RFDesign{21,7}*RFDesign{21,3}/10^6;

        RFDesign{22,1}='T1';
        RFDesign{22,2}='ColOutSup1';
        RFDesign{22,3}=AsColOutSup1*RFincrease;
        RFDesign{22,4}=WColOut1;
        RFDesign{22,5}=XSpan;
        RFDesign{22,6}=4;
        RFDesign{22,7}=RFDesign{22,4}*RFDesign{22,5}*RFDesign{22,6};
        RFDesign{22,8}=RFDesign{22,7}*RFDesign{22,3}/10^6;

        RFDesign{23,1}='T1';
        RFDesign{23,2}='MidIntSup1';
        RFDesign{23,3}=AsMidIntSup1*RFincrease;
        RFDesign{23,4}=WMidStrip1;
        RFDesign{23,5}=XSpan;
        RFDesign{23,6}=4;
        RFDesign{23,7}=RFDesign{23,4}*RFDesign{23,5}*RFDesign{23,6};
        RFDesign{23,8}=RFDesign{23,7}*RFDesign{23,3}/10^6;

        RFDesign{24,1}='T1';
        RFDesign{24,2}='MidIntSup2';
        RFDesign{24,3}=AsMidIntSup2*RFincrease;
        RFDesign{24,4}=WMidStrip2;
        RFDesign{24,5}=XSpan;
        RFDesign{24,6}=4;
        RFDesign{24,7}=RFDesign{24,4}*RFDesign{24,5}*RFDesign{24,6};
        RFDesign{24,8}=RFDesign{24,7}*RFDesign{24,3}/10^6;

        RFDesign{25,1}='T1';
        RFDesign{25,2}='ColOutSup2';
        RFDesign{25,3}=AsColOutSup2*RFincrease;
        RFDesign{25,4}=WColOut2;
        RFDesign{25,5}=XSpan;
        RFDesign{25,6}=4;
        RFDesign{25,7}=RFDesign{25,4}*RFDesign{25,5}*RFDesign{25,6};
        RFDesign{25,8}=RFDesign{25,7}*RFDesign{25,3}/10^6;

        RFDesign{26,1}='T1';
        RFDesign{26,2}='ColMidSup2';
        RFDesign{26,3}=AsColMidSup2*RFincrease;
        RFDesign{26,4}=WColMid2;
        RFDesign{26,5}=XSpan;
        RFDesign{26,6}=4;
        RFDesign{26,7}=RFDesign{26,4}*RFDesign{26,5}*RFDesign{26,6};
        RFDesign{26,8}=RFDesign{26,7}*RFDesign{26,3}/10^6;


        T1FloorArea=sum(cell2mat(RFDesign(18:26,7)));
        T1SteelVolume=sum(cell2mat(RFDesign(18:26,8)));

        %RF Map

        if status==1

            xx=[];yy=[];xxx=[];yyy=[];xxxx=[];yyyy=[];

            xx=[XSpan/2 XSpan XSpan XSpan/2];
            y1(1)=RFDesign{21,4};
            y1(2)=RFDesign{22,4};
            y1(3)=RFDesign{23,4};
            y1(4)=RFDesign{24,4}/2;
            y1(5)=RFDesign{25,4}/2;
            y1(6)=RFDesign{26,4};
            yy=[y1 flip(y1(4:5))];
            yy=[yy yy(4:8) flip(yy(1:3))];

            limxx=length(xx);
            xxx(1)=0;
            for i=1:limxx
                xxx(i+1)=xx(i)+xxx(i);
            end
            limxxx=length(xxx);
            xxxx(1)=0;
            for i=1:limxxx-2
                xxxx(i*2)=xxx(i+1)-0.001;
                xxxx(i*2+1)=xxx(i+1)+0.001;
            end
            xxxx(length(xxxx)+1)=xxx(end);

            limyy=length(yy);
            yyy(1)=0;
            for i=1:limyy
                yyy(i+1)=yy(i)+yyy(i);
            end
            limyyy=length(yyy);
            yyyy(1)=0;
            for i=1:limyyy-2
                yyyy(i*2)=yyy(i+1)-0.001;
                yyyy(i*2+1)=yyy(i+1)+0.001;
            end
            yyyy(length(yyyy)+1)=yyy(end);

            [X,Y]=meshgrid(xxxx,yyyy);

            ZT1=zeros(length(yyyy),length(xxxx));
            ZT1(1:4,1:2)=RFDesign{18,3};
            ZT1(5:8,1:2)=RFDesign{19,3};
            ZT1(9:14,1:2)=RFDesign{20,3};
            ZT1(1:2,3:4)=RFDesign{21,3};
            ZT1(3:4,3:4)=RFDesign{22,3};
            ZT1(5:6,3:4)=RFDesign{23,3};
            ZT1(7:8,3:4)=RFDesign{24,3};
            ZT1(9:10,3:4)=RFDesign{25,3};
            ZT1(10:12,3:4)=RFDesign{26,3};
            ZT1(13:14,3:4)=RFDesign{25,3};
            ZT1(15:16,1:4)=ZT1(7:8,1:4);
            ZT1(17:26,1:4)=ZT1(7:16,1:4);
            ZT1(27:32,1:4)=flip(ZT1(1:6,1:4),1);
            ZT1(1:32,5:8)=flip(ZT1(1:32,1:4),2);


            figure(3);
            surf(X,Y,ZT1);
            colormap 'parula';
            colorbar;
            xlim([0 XSpan*3]);
            ylim([0 YSpan*3]);
            %zlim([0 1500]);
            xticks([0 XSpan XSpan*2 XSpan*3]);
            yticks([0 YSpan YSpan*2 YSpan*3]);
            view(2);
        end


            % Top 2nd layer (T2)

            RFDesign{27,1}='T2';
            RFDesign{27,2}='ColEndSupA';
            RFDesign{27,3}=AsColEndSupA*RFincrease;
            RFDesign{27,4}=WColStripA;
            RFDesign{27,5}=YSpan/2;
            RFDesign{27,6}=4;
            RFDesign{27,7}=RFDesign{27,4}*RFDesign{27,5}*RFDesign{27,6};
            RFDesign{27,8}=RFDesign{27,7}*RFDesign{27,3}/10^6;

            RFDesign{28,1}='T2';
            RFDesign{28,2}='MidEndSupB';
            RFDesign{28,3}=AsMidEndSupB*RFincrease;
            RFDesign{28,4}=WMidStripA+WMidStripB/2;
            RFDesign{28,5}=YSpan/2;
            RFDesign{28,6}=6;
            RFDesign{28,7}=RFDesign{28,4}*RFDesign{28,5}*RFDesign{28,6};
            RFDesign{28,8}=RFDesign{28,7}*RFDesign{28,3}/10^6;

            RFDesign{29,1}='T2';
            RFDesign{29,2}='ColEndSupB';
            RFDesign{29,3}=AsColEndSupB*RFincrease;
            RFDesign{29,4}=WColStripB;
            RFDesign{29,5}=YSpan/2;
            RFDesign{29,6}=4;
            RFDesign{29,7}=RFDesign{29,4}*RFDesign{29,5}*RFDesign{29,6};
            RFDesign{29,8}=RFDesign{29,7}*RFDesign{29,3}/10^6;

            RFDesign{30,1}='T2';
            RFDesign{30,2}='ColMidSupA';
            RFDesign{30,3}=AsColMidSupA*RFincrease;
            RFDesign{30,4}=WColMidA;
            RFDesign{30,5}=YSpan;
            RFDesign{30,6}=4;
            RFDesign{30,7}=RFDesign{30,4}*RFDesign{30,5}*RFDesign{30,6};
            RFDesign{30,8}=RFDesign{30,7}*RFDesign{30,3}/10^6;

            RFDesign{31,1}='T2';
            RFDesign{31,2}='ColOutSupA';
            RFDesign{31,3}=AsColOutSupA*RFincrease;
            RFDesign{31,4}=WColOutA;
            RFDesign{31,5}=YSpan;
            RFDesign{31,6}=4;
            RFDesign{31,7}=RFDesign{31,4}*RFDesign{31,5}*RFDesign{31,6};
            RFDesign{31,8}=RFDesign{31,7}*RFDesign{31,3}/10^6;

            RFDesign{32,1}='T2';
            RFDesign{32,2}='MidIntSupA';
            RFDesign{32,3}=AsMidIntSupA*RFincrease;
            RFDesign{32,4}=WMidStripA;
            RFDesign{32,5}=YSpan;
            RFDesign{32,6}=4;
            RFDesign{32,7}=RFDesign{32,4}*RFDesign{32,5}*RFDesign{32,6};
            RFDesign{32,8}=RFDesign{32,7}*RFDesign{32,3}/10^6;

            RFDesign{33,1}='T2';
            RFDesign{33,2}='MidIntSupB';
            RFDesign{33,3}=AsMidIntSupB*RFincrease;
            RFDesign{33,4}=WMidStripB;
            RFDesign{33,5}=YSpan;
            RFDesign{33,6}=4;
            RFDesign{33,7}=RFDesign{33,4}*RFDesign{33,5}*RFDesign{33,6};
            RFDesign{33,8}=RFDesign{33,7}*RFDesign{33,3}/10^6;

            RFDesign{34,1}='T2';
            RFDesign{34,2}='ColOutSupB';
            RFDesign{34,3}=AsColOutSupB*RFincrease;
            RFDesign{34,4}=WColOutB;
            RFDesign{34,5}=YSpan;
            RFDesign{34,6}=4;
            RFDesign{34,7}=RFDesign{34,4}*RFDesign{34,5}*RFDesign{34,6};
            RFDesign{34,8}=RFDesign{34,7}*RFDesign{34,3}/10^6;

            RFDesign{35,1}='T2';
            RFDesign{35,2}='ColMidSupB';
            RFDesign{35,3}=AsColMidSupB*RFincrease;
            RFDesign{35,4}=WColMidB;
            RFDesign{35,5}=YSpan;
            RFDesign{35,6}=4;
            RFDesign{35,7}=RFDesign{35,4}*RFDesign{35,5}*RFDesign{35,6};
            RFDesign{35,8}=RFDesign{35,7}*RFDesign{35,3}/10^6;


            T2FloorArea=sum(cell2mat(RFDesign(27:35,7)));
            T2SteelVolume=sum(cell2mat(RFDesign(27:35,8)));

            %RF Map

            if status==1

                xx=[];yy=[];xxx=[];yyy=[];xxxx=[];yyyy=[];

                yy=[YSpan/2 YSpan YSpan YSpan/2];
                x1(1)=RFDesign{30,4};
                x1(2)=RFDesign{31,4};
                x1(3)=RFDesign{32,4};
                x1(4)=RFDesign{33,4}/2;
                x1(5)=RFDesign{34,4}/2;
                x1(6)=RFDesign{35,4};
                xx=[x1 flip(x1(4:5))];
                xx=[xx xx(4:8) flip(xx(1:3))];

                limyy=length(yy);
                yyy(1)=0;
                for i=1:limyy
                    yyy(i+1)=yy(i)+yyy(i);
                end
                limyyy=length(yyy);
                yyyy(1)=0;
                for i=1:limyyy-2
                    yyyy(i*2)=yyy(i+1)-0.001;
                    yyyy(i*2+1)=yyy(i+1)+0.001;
                end
                yyyy(length(yyyy)+1)=yyy(end);

                limxx=length(xx);
                xxx(1)=0;
                for i=1:limxx
                    xxx(i+1)=xx(i)+xxx(i);
                end
                limxxx=length(xxx);
                xxxx(1)=0;
                for i=1:limxxx-2
                    xxxx(i*2)=xxx(i+1)-0.001;
                    xxxx(i*2+1)=xxx(i+1)+0.001;
                end
                xxxx(length(xxxx)+1)=xxx(end);

                [X,Y]=meshgrid(xxxx,yyyy);

                ZT2=zeros(length(yyyy),length(xxxx));
                ZT2(1:2,1:4)=RFDesign{27,3};
                ZT2(1:2,5:8)=RFDesign{28,3};
                ZT2(1:2,9:14)=RFDesign{29,3};
                ZT2(3:4,1:2)=RFDesign{30,3};
                ZT2(3:4,3:4)=RFDesign{31,3};
                ZT2(3:4,5:6)=RFDesign{32,3};
                ZT2(3:4,7:8)=RFDesign{33,3};
                ZT2(3:4,9:10)=RFDesign{34,3};
                ZT2(3:4,10:12)=RFDesign{35,3};
                ZT2(3:4,13:14)=RFDesign{34,3};
                ZT2(1:4,15:16)=ZT2(1:4,7:8);
                ZT2(1:4,17:26)=ZT2(1:4,7:16);
                ZT2(1:4,27:32)=flip(ZT2(1:4,1:6),2);
                ZT2(5:8,1:32)=flip(ZT2(1:4,1:32),1);

                figure(4);
                surf(X,Y,ZT2);
                colormap 'parula';
                colorbar;
                xlim([0 XSpan*3]);
                ylim([0 YSpan*3]);
                %zlim([0 1500]);
                xticks([0 XSpan XSpan*2 XSpan*3]);
                yticks([0 YSpan YSpan*2 YSpan*3]);
                view(2);
            end

            for i=2:35
                RFDesign{i,9}=1000*(3.142*6*6)/RFDesign{i,3};
                RFDesign{i,10}=1000*(3.142*10*10)/RFDesign{i,3};

            end


        % SHEAR REINFORCEMENT DESIGN

        %Three column types, internal, corner and x and y sides
        %Int=B2,B3,C2,C3
        %Corn=A1,A4,D1,D4
        %Side=B1,B4,C1,C4,A2,A3,D2,D3

        ShearDesign{1,1}='Column Type';
        ShearDesign{1,2}='Design Load';
        ShearDesign{1,3}='Column Perimiter';
        ShearDesign{1,4}='Control Perimiter';
        ShearDesign{1,5}='ved';
        ShearDesign{1,6}='uout';
        ShearDesign{1,7}='Ruout';
        ShearDesign{1,8}='Number of Perimiters';
        ShearDesign{1,9}='Number of Links';
        ShearDesign{1,10}='Steel volume';

        VdInt=ULSLoad*XSpan*YSpan*0.60*2;
        VdCorn=ULSLoad*XSpan*YSpan*0.5*0.46;
        VdSide=ULSLoad*XSpan*YSpan*0.46;

        dm=(dx+dy)/2;
        k1=1+(200/dm)^0.5;
        k2=2;
        k=min(k1,k2);
        Srmax=dm*0.75;
        Stmaxin=dm*1.5;
        Stmaxout=dm*2;
        fyw1=250+0.25*dm;
        fyw2=0.87*fsy;
        fyw=min(fyw1,fyw2);
        vrdmax=0.2*(1-(fck/250))*fck;

        %Design of Int Column

        beta=1.15;
        ui=2*(ColX+ColY);
        ved0=beta*VdInt*1000/(ui*dm);
        if ved0>vrdmax
            disp('Shear perimiter failure');
            status=0;
            StatShearF=0;
            condition3='Shear perimiter failure';
        end

        u1=2*(ColX+ColY)+2*3.142*2*dm;
        ved=beta*VdInt*1000/(u1*dm);

        rox=AsColIntSup2*RFincrease/(1000*dx);
        roy=AsColIntSupB*RFincrease/(1000*dy);
        ro1=(rox*roy)^0.5;
        ro2=0.02;
        ro=min(ro1,ro2);
        vrdc1=0.12*k*(100*ro*fck)^(1/3);
        vrdc2=0.035*(k^1.5)*(fck^0.5);
        vrdc=max(vrdc1,vrdc2);

        VratioInt=ved/vrdc;

        uout=VdInt*beta*1000/(dm*vrdc);
        Ruout=(uout-ui)/(2*3.142);

        Asw=(ved-0.75*vrdc)*Srmax*u1/(1.5*fyw);
        Aswmin=max(0.08*(fck^0.5)*(Srmax*Stmaxin)/(1.5*fsy),50); %at least use H8
        Stin=min(Aswmin*u1/Asw, Stmaxin);
        Nperiin=ceil((2*dm-0.4*dm)/Srmax);

        Aswmin=max(0.08*(fck^0.5)*(Srmax*Stmaxin)/(1.5*fsy),50); %at least use H8
        Stout=min(Aswmin*u1/Asw, Stmaxout);
        %Nperiout=ceil((Ruout-2*dm)/Srmax);
        Ntot=ceil((Ruout-1.5*dm-0.4*dm)/Srmax)+1;
        %Ntot=Nperiin+Nperiout;
        links=0;
        for i=1:Ntot
            %ColAvg=(ColX+ColY)/2;
            %Peri=(ColAvg/2+0.4*dm+Srmax*(i-1))*8;
            Peri=2*(ColX+ColY)+2*3.142*(0.4*dm+Srmax*(i-1));
            if i>Nperiin
                N=ceil(Peri/Stout);
            else
                N=ceil(Peri/Stin);
            end
            Ans=[num2str(i) ' perimiter with ' num2str(Peri) ' length needs ' num2str(N) ' links'];
            %disp(Ans);
            links=links+N;
        end

        if ved<vrdc
            disp('Punching shear RF not required');
            links=0;
        else if ved>2*vrdc
                disp ('Punching redesign, ved>2vrdc');
            end
        end

        SLinksInt=links;
        SLinksSteelInt=SLinksInt*Aswmin*h/(10^9);

        ShearDesign{2,1}='Internal';
        ShearDesign{2,2}=VdInt;
        ShearDesign{2,3}=Aswmin;
        ShearDesign{2,4}=u1;
        ShearDesign{2,5}=ved;
        ShearDesign{2,6}=uout;
        ShearDesign{2,7}=Ruout;
        ShearDesign{2,8}=Ntot;
        ShearDesign{2,9}=SLinksInt;
        ShearDesign{2,10}=SLinksSteelInt;

            %Design of Corn Column

            beta=1.5;
            ui1=3*dm;
            ui2=ColX+ColY;
            ui=min(ui1,ui2);
            ved0=beta*VdCorn*1000/(ui*dm);
            if ved0>vrdmax
                disp('Shear perimiter failure');
                status=0;
                StatShearF=0;
                condition3='Shear perimiter failure';
            end

            u1=ui+3.142*2*dm;
            ved=beta*VdCorn*1000/(u1*dm);

            rox=AsColEndSup1*RFincrease/(1000*dx);
            roy=AsColEndSupA*RFincrease/(1000*dy);
            ro1=(rox*roy)^0.5;
            ro2=0.02;
            ro=min(ro1,ro2);
            vrdc1=0.12*k*(100*ro*fck)^(1/3);
            vrdc2=0.035*(k^1.5)*(fck^0.5);
            vrdc=max(vrdc1,vrdc2);

            VratioCorn=ved/vrdc;

            uout=VdCorn*beta*1000/(dm*vrdc);
            Ruout=(uout-ui2)/(3.142*0.5);

            Asw=(ved-0.75*vrdc)*Srmax*u1/(1.5*fyw);
            Aswmin=max(0.08*(fck^0.5)*(Srmax*Stmaxin)/(1.5*fsy),50); %at least use H8
            Stin=min(Aswmin*u1/Asw, Stmaxin);
            Nperiin=ceil((2*dm-0.4*dm)/Srmax);

            Aswmin=max(0.08*(fck^0.5)*(Srmax*Stmaxin)/(1.5*fsy),50); %at least use H8
            Stout=min(Aswmin*u1/Asw, Stmaxout);
            %Nperiout=ceil((Ruout-2*dm)/Srmax);
            Ntot=ceil((Ruout-(1.5+0.4)*dm)/Srmax)+1;
            %Ntot=Nperiin+Nperiout;
            links=0;
            for i=1:Ntot
                ColAvg=(ColX+ColY)/2;
                %Peri=(ColAvg/2+0.4*dm+Srmax*(i-1))*2;
                Peri=(ColX+ColY)+0.5*3.142*(0.4*dm+Srmax*(i-1));
                if i>Nperiin
                    N=ceil(Peri/Stout);
                else
                    N=ceil(Peri/Stin);
                end
                Ans=[num2str(i) ' perimiter with ' num2str(Peri) ' length needs ' num2str(N) ' links'];
                %disp(Ans);
                links=links+N;
            end

            if ved<vrdc
                disp('Punching shear RF not required');
                links=0;
            else if ved>2*vrdc
                    disp ('Punching redesign, ved>2vrdc');
                end
            end

            SLinksCorn=links;
            SLinksSteelCorn=SLinksCorn*Aswmin*h/(10^9);

            ShearDesign{3,1}='Corner';
            ShearDesign{3,2}=VdCorn;
            ShearDesign{3,3}=Aswmin;
            ShearDesign{3,4}=u1;
            ShearDesign{3,5}=ved;
            ShearDesign{3,6}=uout;
            ShearDesign{3,7}=Ruout;
            ShearDesign{3,8}=Ntot;
            ShearDesign{3,9}=SLinksCorn;
            ShearDesign{3,10}=SLinksSteelCorn;

        %Design of Side Column

        beta=1.4;
        ui1=3*dm+ColX;
        ui2=ColX+2*ColY;
        ui=min(ui1,ui2);
        ved0=beta*VdSide*1000/(ui*dm);
        if ved0>vrdmax
            disp('Shear perimiter failure');
            status=0;
            StatShearF=0;
        end

        u1=(ColX+2*ColY)+3.142*2*dm;
        ved=beta*VdSide*1000/(u1*dm);

        rox=AsColEndSup2*RFincrease/(1000*dx);
        roy=AsColEndSupB*RFincrease/(1000*dy);
        ro1=(rox*roy)^0.5;
        ro2=0.02;
        ro=min(ro1,ro2);
        vrdc1=0.12*k*(100*ro*fck)^(1/3);
        vrdc2=0.035*(k^1.5)*(fck^0.5);
        vrdc=max(vrdc1,vrdc2);

        VratioSide=ved/vrdc;

        uout=VdSide*beta*1000/(dm*vrdc);
        Ruout=(uout-ui2)/(3.142);

        Asw=(ved-0.75*vrdc)*Srmax*u1/(1.5*fyw);
        Aswmin=max(0.08*(fck^0.5)*(Srmax*Stmaxin)/(1.5*fsy),50); %at least use H8
        Stin=min(Aswmin*u1/Asw, Stmaxin);
        Nperiin=ceil((2*dm-0.4*dm)/Srmax);

        Aswmin=max(0.08*(fck^0.5)*(Srmax*Stmaxin)/(1.5*fsy),50); %at least use H8
        Stout=min(Aswmin*u1/Asw, Stmaxout);
        %Nperiout=ceil((Ruout-2*dm)/Srmax);
        Ntot=ceil((Ruout-(1.5+0.4)*dm)/Srmax)+1;
        %Ntot=Nperiin+Nperiout;
        links=0;
        for i=1:Ntot
            ColAvg=(ColX+ColY)/2;
            %Peri=(ColAvg/2+0.4*dm+Srmax*(i-1))*4;
            Peri=2*ColX+ColY+3.142*(0.4*dm+Srmax*(i-1));
            if i>Nperiin
                N=ceil(Peri/Stout);
            else
                N=ceil(Peri/Stin);
            end
            Ans=[num2str(i) ' perimiter with ' num2str(Peri) ' length needs ' num2str(N) ' links'];
            %disp(Ans);
            links=links+N;
        end

        if ved<vrdc
            disp('Punching shear RF not required');
            links=0;
        else if ved>2*vrdc
                disp ('Punching redesign, ved>2vrdc');
            end
        end

        SLinksSide=links;
        SLinksSteelSide=SLinksSide*Aswmin*h/(10^9);

        ShearDesign{4,1}='Side';
        ShearDesign{4,2}=VdSide;
        ShearDesign{4,3}=Aswmin;
        ShearDesign{4,4}=u1;
        ShearDesign{4,5}=ved;
        ShearDesign{4,6}=uout;
        ShearDesign{4,7}=Ruout;
        ShearDesign{4,8}=Ntot;
        ShearDesign{4,9}=SLinksSide;
        ShearDesign{4,10}=SLinksSteelSide;

        Vratios=[VratioInt VratioCorn VratioSide];
        condition4=[num2str(Vratios) ' ved/vrdc of Internal, Corner and Side'];
        StatShearC=max(Vratios);


        %COLUMN DESIGN

        hc=ColX-cover-10-10;
        bc=ColY;

        ColumnDesign={};
        ColTypeM={'Internal' 'Corner' 'Side'};
        ColumnDesign{1,1}='Column Type';
        ColumnDesign{1,2}='Width';
        ColumnDesign{1,3}='V kN';
        ColumnDesign{1,4}='M kNm';
        ColumnDesign{1,5}='RF';
        ColumnDesign{1,6}='T20 Nos';
        ColumnDesign{1,7}='T10 link spacing';
        ColumnDesign{1,8}='RF m3';
        ColumnDesign{1,9}='Links m3';
        ColumnDesign{1,10}='Concrete m3';


        %Design Forces
        %Vall=[VdInt VdCorn VdSide]*4;   %Consider load of 3 storeys
        Vall=[ULSLoad*XSpan*YSpan*3 ULSLoad*XSpan*YSpan*0.25*3 ULSLoad*XSpan*YSpan*0.5*3];
        Mall=[0 0.25*MColEndSpanA*0.5*XSpan 0.25*MColEndSpanB*0.5*XSpan];

        for c=1:3

            Asc=NaN;

            V=Vall(c);
            M=Mall(c);
            ColType=ColTypeM(c);

            Md0=M+V*CL/400;
            Mdmin=V*max(20,hc/30)/1000;
            Md=max(Md0, Mdmin)*2; %*2 for biaxial bending

            Mratio=Md*1000000/(fckcol*bc*hc*hc);
            Nratio=V*1000/(bc*hc*fckcol);

            Asratio=interp2(MMCH,NNCH,Ach,Mratio,Nratio);
            %Asratio=interp2(MMCH,NNCH,Ach,0.225,0.3);

            Ask=Asratio*bc*hc*fckcol/fsy;

            Asmin=max((0.10*V*1000/(fsy*0.87)),0.002*ColX*ColY);
            Asmax=0.04*ColX*ColY;

            if (Asratio==1) || (Ask>0.04*ColX*ColY) || (Nratio>1.3)
             disp('Redesign Column size');
             status=0;
             StatColumn=0;
            end

            Asc=min(max(Ask,Asmin),Asmax);

            %provide links with T10

            sp=min([12*20, 0.6*ColX, 240]);
            Nlinks=ceil(CL*1000/sp);
            Vonelink=4*(ColX-cover*2)*3.142*5*5/(10^9);
            Vlink=Vonelink*Nlinks;

            VRF=Asc*CL/(10^6);

            VCon=(ColX*ColY*CL/(10^6))-VRF-Vlink;

            ECol=VCon*2400*0.132 + (VRF+Vlink)*7850*1.2;

            ColumnDesign{c+1,1}=ColType;
            ColumnDesign{c+1,2}=ColX;
            ColumnDesign{c+1,3}=V;
            ColumnDesign{c+1,4}=Md;
            ColumnDesign{c+1,5}=Asc;
            ColumnDesign{c+1,6}=Asc/(3.142*100);
            ColumnDesign{c+1,7}=sp;
            ColumnDesign{c+1,8}=VRF;
            ColumnDesign{c+1,9}=Vlink;
            ColumnDesign{c+1,10}=VCon;

        end




        %Embodeid Carbon Calculation

        if status==1

            VSteelShear = ShearDesign{2,10}*4 + ShearDesign{3,10}*4 + ShearDesign{4,10}*8;

            DataVSteel=cell2mat(RFDesign(2:35,8));
            VSteelT1=sum(DataVSteel(17:26));
            VSteelT2=sum(DataVSteel(27:34));
            VSteelB1=sum(DataVSteel(1:8));
            VSteelB2=sum(DataVSteel(9:16));

            VSteelFlex=VSteelT1+VSteelT2+VSteelB1+VSteelB2;

            VSteelTotal=VSteelFlex+VSteelShear;
            
            %Curtailing
            
            VSteelT1Mid=DataVSteel(18)+DataVSteel(22)+DataVSteel(23);
            VSteelT1Col=VSteelT1-VSteelT1Mid;
            VSteelT2Mid=DataVSteel(27)+DataVSteel(31)+DataVSteel(32);
            VSteelT2Col=VSteelT2-VSteelT2Mid;
            
            VSteelTCol=VSteelT1Col+VSteelT2Col;
            VSteelTMid=VSteelT1Mid+VSteelT2Mid;
            VSteelB=VSteelB1+VSteelB2;
            
            %Minimum for 0.34L for two larers on Top net
            CurtMins=AsminC*0.34*XSpan*3*YSpan*3*0.5*2*(10^(-6));
            VSteelTCurt=CurtMins+VSteelTCol*0.435+VSteelTMid;
            VSteelFlexCurt=VSteelTCurt+VSteelB;
            VSteelTotalCurt=VSteelFlexCurt+VSteelShear;

            VConcrete=(XSpan*YSpan*9*h/1000)-VSteelTotalCurt;

            ESteelFlexCurt=VSteelFlexCurt*7850*ecSteel;
            ESteelShear=VSteelShear*7850*ecSteel;
            EConcrete=VConcrete*2400*ecSlabCon;

            VColRF = ColumnDesign{2,8}*4 + ColumnDesign{3,8}*4 + ColumnDesign{4,8}*8;
            VColLink= ColumnDesign{2,9}*4 + ColumnDesign{3,9}*4 + ColumnDesign{4,9}*8;
            VColConc= ColumnDesign{2,10}*4 + ColumnDesign{3,10}*4 + ColumnDesign{4,10}*8;
            EColumn=VColConc*2400*ecColumnCon + (VColRF+VColLink)*7850*ecSteel;

            ESteelTotal=ESteelFlexCurt+ESteelShear;
            ETotal=EConcrete+ESteelTotal+EColumn;

            Ans=['For Span = ' num2str(SP) 'For h = ' num2str(h) ' Total Embodied carbon is ' num2str(ETotal)];
            disp(Ans);
        else
            ETotal=0;
            EConcrete=0;
            ESteelFlex=0;
            ESteelShear=0;
            EColumn=0;
            VSteelB1=0;
            VSteelB2=0;
            VSteelT1=0;
            VSteelT2=0;

            Ans=['For Span = ' num2str(SP) 'For h = ' num2str(h) ' Total Embodied carbon is ' num2str(ETotal)];
            disp(Ans);
        end

        ECM{beam,1}=h;
        ECM{beam,2}=EConcrete/(XSpan*YSpan*9);
        ECM{beam,3}=ESteelFlexCurt/(XSpan*YSpan*9);
        ECM{beam,4}=ESteelShear/(XSpan*YSpan*9);
        ECM{beam,5}=EColumn/(XSpan*YSpan*9);
        ECM{beam,6}=ETotal/(XSpan*YSpan*9);
        ECM{beam,7}=condition1;
        ECM{beam,8}=condition2;
        ECM{beam,9}=condition3;
        ECM{beam,10}=condition4;

        if XSpan>8.5
            LDM=LDM*(8.5/XSpan);
        end

        hdef=XSpan*1000/min(LDM);
        if h>=hdef
            condition5='viable';
        else
            condition5=[num2str(hdef) ' mm is needed at least'];
            StatDef=0;
        end
        ECM{beam,11}=condition5;

        ECM{beam,12}=max(cell2mat(RFDesign(2:17,3)));
        ECM{beam,13}=max(cell2mat(RFDesign(18:35,3)));
        ECM{beam,14}=(VSteelB1+VSteelB2)*(10^6)/(XSpan*YSpan*9*2);
        ECM{beam,15}=(VSteelT1+VSteelT2)*(10^6)/(XSpan*YSpan*9*2);
        ECM{beam,16}=min(cell2mat(RFDesign(2:17,3)));
        ECM{beam,17}=min(cell2mat(RFDesign(18:35,3)));

        ECMT{beam,1}=h;
        ECMT{beam,2}=EConcrete;
        ECMT{beam,3}=ESteelFlexCurt;
        ECMT{beam,4}=ESteelShear;
        ECMT{beam,5}=EColumn;
        ECMT{beam,6}=ETotal;
        ECMT{beam,7}=condition1;
        ECMT{beam,8}=condition2;
        ECMT{beam,9}=condition3;
        ECMT{beam,10}=condition4;
        ECMT{beam,11}=condition5;

        ECSlab(beam-1)=ECM{beam,6};
        
        StatCompRFMat(beam-1)=StatCompRF;
        StatDefMat(beam-1)=StatDef;
        StatColumnMat(beam-1)=StatColumn;
        StatShearFMat(beam-1)=StatShearF;
        StatShearCMat(beam-1)=StatShearC;

        %RFDesign2=[RFDesign num2cell(zeros(35,2))];
        %ThisRFDetails=[RFDesign2; ShearDesign];
        ThisRFDetails=[RFDesign; ShearDesign; ColumnDesign];
        ThisBeam={'For h=' h 0 0 0 0 0 0 0 0};

        AllRFDetails=[AllRFDetails; ThisBeam; ThisRFDetails];
        
        SlabNum(beam-1)=h;

    end

    %PlotEC graphs

    H=cell2mat(ECM(2:end,1));
    ECon=cell2mat(ECM(2:end,2));
    ECon2=ECon;
    ECon2(ECon2==0)=nan;
    EFlex=cell2mat(ECM(2:end,3));
    EFlex2=EFlex;
    EFlex2(EFlex2==0)=nan;
    EShear=cell2mat(ECM(2:end,4));
    EShear2=EShear;
    EShear2(EShear2==0)=nan;
    ECol=cell2mat(ECM(2:end,5));
    ECol2=ECol;
    ECol2(ECol2==0)=nan;
    ETot=cell2mat(ECM(2:end,6));
    ETot2=ETot;
    ETot2(ETot2==0)=nan;

    %q=find(~isnan(C1Pris),1);
    %H=H(q:end);
    %C1Pris=C1Pris(q:end);

    figure(5);
    h(1)=plot(H,ECon2,'-r','LineWidth',2);
    hold on;
    h(2)=plot(H,EFlex2,'--b','LineWidth',2);
    hold on;
    h(3)=plot(H,EShear2,'-g','LineWidth',2);
    hold on;
    h(4)=plot(H,ECol2,':y','LineWidth',2);
    hold on;
    h(5)=plot(H,ETot2,':k','LineWidth',2);
    hold on;
    [m,I]=min(ETot2);
    h(6)=plot(H(I),m,'vr','MarkerFaceColor','y');
    hold on;
    h(7)=plot([H(I) H(I)],[m 0],'y');
    hold on;
    h(8)=plot(H(I),0,'vr','MarkerFaceColor','y');
    hold on;

    xlabel('Slab depth(mm)');
    ylabel('Embodied Carbon (kgCO2/m2)');
    xlim([hmin hmax]);
    ylim([0 200]);
    legend(h([1 2 3 4 5]),'Concrete','Flex RF','Shear RF','Column','Total');
    legend('Location','southoutside');
    set(legend, 'NumColumns' ,3);
    hold on;
    
    saveas(gcf,'C30 7m 400mm Col Curt EC.jpg');
    saveas(gcf,'C30 7m 400mm Col Curt EC.fig');
    saveas(gcf,'C30 7m 400mm Col Curt EC.pdf');

    
    AllECMat=[AllECMat;ECSlab];
    
    AllCompRFMat=[AllCompRFMat; StatCompRFMat];
    AllDefMat=[AllDefMat; StatDefMat];
    AllColumnMat=[AllColumnMat; StatColumnMat];
    AllShearFMat=[AllShearFMat; StatShearFMat];
    AllShearCMat=[AllShearCMat; StatShearCMat];
    
    SpanNum(SpanN)=SP;
    SpanN=SpanN+1;

end

%{

AllECMat(AllECMat==0)=nan;

[THICKNESS,SPAN]=meshgrid(SlabNum,SpanNum);
figure(6);
%contourf(SPAN,THICKNESS,AllECMat,30);
contourf(SPAN,THICKNESS,AllECMat,25,'LineColor','none');
hold on;
ylabel('Slab Depth (mm)');
xlabel('Column Spacing (m)');
zlabel('Embodied Carbon (kgCO2/m2)');
colorbar;
colormap('jet');
ylim([hmin hmax]);
xlim([SpanMin SpanMax]);
view(2);

%Compressive RF line
TComRFMat=[];
for n=1:size(AllCompRFMat,1)
    if (sum(AllCompRFMat(n,:))==0)
        TComRFMat(n)=hmax+100;
        ECComRFMat(n)=500;
    else
        TComRF=find(AllCompRFMat(n,:), 1, 'first');
        TComRFMat(n)=SlabNum(TComRF);
        ECComRFMat(n)=500;
    end
end
h(1)=plot3(SpanNum,TComRFMat,ECComRFMat,'-om','LineWidth',1,'MarkerFaceColor','m','MarkerEdgeColor','m');
hold on;

%Deflection line
TDefMat=[];
for n=1:size(AllDefMat,1)
    if (sum(AllDefMat(n,:))==0)
        TDefMat(n)=hmax+100;
        ECDefMat(n)=500;
        TDefMat2(n)=NaN;
        ECDefMat2(n)=NaN;
    else
        TDef=find(AllDefMat(n,:), 1, 'first');
        TDefMat(n)=SlabNum(TDef);
        ECDefMat(n)=500;
        TDefMat2(n)=SlabNum(TDef);
        ECDefMat2(n)=AllECMat(n,TDef);
    end
end
h(2)=plot3(SpanNum,TDefMat,ECDefMat,'-vw','LineWidth',2,'MarkerFaceColor','w','MarkerEdgeColor','y');
hold on;
    
%Column Failure line
TColumnMat=[];
for n=1:size(AllColumnMat,1)
    if (sum(AllColumnMat(n,:))==0)
        TColumnMat(n)=hmin-100;
        ECColumnMat(n)=500;
    else
        TColumn=find(AllColumnMat(n,:), 1, 'last');
        TColumnMat(n)=SlabNum(TColumn);
        ECColumnMat(n)=500;
    end
end
h(3)=plot3(SpanNum,TColumnMat,ECColumnMat,'-sk','LineWidth',1,'MarkerFaceColor','k','MarkerEdgeColor','k');
hold on;

%Shear Perimeter failure line
TShearFMat=[];
for n=1:size(AllShearFMat,1)
    if (sum(AllShearFMat(n,:))==0)
        TShearFMat(n)=hmax+100;
        ECShearFMat(n)=500;
    else
        TShearF=find(AllShearFMat(n,:), 1, 'first');
        TShearFMat(n)=SlabNum(TShearF);
        ECShearFMat(n)=500;
    end
end
h(4)=plot3(SpanNum,TShearFMat,ECShearFMat,'--sk','LineWidth',1,'MarkerEdgeColor','k');
hold on;

SpansC=[4 5 6 7 8 9 10];
Economic=[200 200 206 227 250 286 343];
SpantoDepth28=[200 200 214 250 285 321 357];
ECDummy=[500 500 500 500 500 500 500];

%h(5)=plot3(SpansC,Economic,ECDummy,'-.','LineWidth',1,'color','#A2142F');
%hold on;
h(6)=plot3(SpansC,SpantoDepth28,ECDummy,'--','LineWidth',1,'color','#D95319');
hold on;


legend(h([1 2 3 4 6]),'Ductile failure','Deflection','Column failure','Shrear failure','S/D=28');
legend('Location','southoutside');
set(legend, 'NumColumns' ,2);
set(gcf,'Position',[200 200 550 550]);
caxis([90 210]);
set(gcf,'Renderer','Painter');
%saveas(gcf,'C30 500mmcol EC-2.png');
%saveas(gcf,'C30 500mmcol EC-2.fig');
%saveas(gcf,'C30 500mmcol EC-2.pdf');


%Plot Graphs for Min EC designs with Col Spacing

%Minimum Possible EC
[MinECs,MinECDepthIndex]=min(AllECMat,[],2);
MinECDepths=SlabNum(MinECDepthIndex);

%Eurocode Deflection span/depth limit
%Within the loop

%Span/Depth=28
RecDepths=max(200,SpanNum*1000/28);
RecECs=interp2(SlabNum,SpanNum,AllECMat,RecDepths,SpanNum);

AllMinECsMat=[AllMinECsMat;transpose(MinECs);ECDefMat2;RecECs];
AllMinECDepthsMat=[AllMinECDepthsMat;MinECDepths;TDefMat2;RecDepths];
%}

 end
 
 %{
 
 AllMinECsMatDup=AllMinECsMat;
 AllMinECsMatDup(isnan(AllMinECsMatDup))=0;
 AllMinECsMatDup(~AllMinECsMatDup==0)=1;
 AllMinECDepthsMatF=AllMinECDepthsMat.*AllMinECsMatDup;
 AllMinECDepthsMatF(AllMinECDepthsMatF==0)=NaN;

%Min EC Graphs
figure(7);
h(1)=plot(SpanNum,AllMinECsMat(1,:),'-g','LineWidth',1.5);
hold on;
h(2)=plot(SpanNum,AllMinECsMat(2,:),'--g','LineWidth',1.5);
hold on;
%h(3)=plot(SpanNum,AllMinECsMat(3,:),':g','LineWidth',1.5);
hold on;
h(4)=plot(SpanNum,AllMinECsMat(4,:),'-b','LineWidth',1.5);
hold on;
h(5)=plot(SpanNum,AllMinECsMat(5,:),'--b','LineWidth',1.5);
hold on;
h(6)=plot(SpanNum,AllMinECsMat(6,:),':k','LineWidth',1.5);
hold on;
h(7)=plot(SpanNum,AllMinECsMat(7,:),'-r','LineWidth',1.5);
hold on;
h(8)=plot(SpanNum,AllMinECsMat(8,:),'--r','LineWidth',1.5);
hold on;
%h(9)=plot(SpanNum,AllMinECsMat(9,:),':r','LineWidth',1.5);
hold on;
ylabel('Embodied Carbon (kgCO2/m2)');
xlabel('Column Spacing (m)');
ylim([0 210]);
xlim([SpanMin SpanMax]);
legend(h([1 2 4 5 7 8 6]),'C20 Optimum for ULS','C20 Calculate Deflection','C30 Optimum for ULS','C30 Calculate Deflection','C40 Optimum for ULS','C40 Calculate Deflection','C30 Span/Depth=28');
legend('Location','southoutside');
set(legend, 'NumColumns' ,2);
set(gcf,'Position',[200 200 550 550]);
supersizeme(1.1);
saveas(gcf,'Min EC with 600mmm Col RF 1.00.png');
saveas(gcf,'Min EC with 600mmm Col RF 1.00.fig');
saveas(gcf,'Min EC with 600mmm Col RF 1.00.pdf');


%Min EC depths graph
figure(8);
h(1)=plot(SpanNum,AllMinECDepthsMatF(1,:),'-g','LineWidth',1.5);
hold on;
h(2)=plot(SpanNum,AllMinECDepthsMatF(2,:),'--g','LineWidth',1.5);
hold on;
%h(3)=plot(SpanNum,AllMinECDepthsMatF(3,:),':g','LineWidth',1.5);
hold on;
h(4)=plot(SpanNum,AllMinECDepthsMatF(4,:),'-b','LineWidth',1.5);
hold on;
h(5)=plot(SpanNum,AllMinECDepthsMatF(5,:),'--b','LineWidth',1.5);
hold on;
%h(6)=plot(SpanNum,AllMinECDepthsMatF(6,:),':b','LineWidth',1.5);
hold on;
h(7)=plot(SpanNum,AllMinECDepthsMatF(7,:),'-r','LineWidth',1.5);
hold on;
h(8)=plot(SpanNum,AllMinECDepthsMatF(8,:),'--r','LineWidth',1.5);
hold on;
h(9)=plot(SpanNum,AllMinECDepthsMatF(9,:),':k','LineWidth',1.5);
hold on;
ylabel('Slab Depth (mm)');
xlabel('Column Spacing (m)');
ylim([0 hmax]);
xlim([SpanMin SpanMax]);
legend(h([1 2 4 5 7 8 9]),'C20 Optimum for ULS','C20 Calculate Deflection','C30 Optimum for ULS','C30 Calculate Deflection','C40 Optimum for ULS','C40 Calculate Deflection','Span/Depth=28');
legend('Location','southoutside');
set(legend, 'NumColumns' ,2);
set(gcf,'Position',[200 200 550 550]);
supersizeme(1.1);
saveas(gcf,'Min EC Depths with 600mmm Col RF 1.00.png');
saveas(gcf,'Min EC Depths with 600mmm Col RF 1.00.fig');
saveas(gcf,'Min EC Depths with 600mmm Col RF 1.00.pdf');

 %}
 