%Zbrojenie belki
clear, clc
% Instrukcje wprowadzania danych do obliczeñ
Dane=inputdlg({'Œrednica prêtów [mm]','Œrednica strzemion [mm]','Gruboœæ otuliny [mm]','Wytrzyma³oœæ betonu [MPa]','Granica plastycznoœci stali [MPa]','Modu³ sprê¿ystoœci stali [GPa]'});
fi=str2double(Dane{1}); fis=str2double(Dane{2}); c=str2double(Dane{3}); fck=str2double(Dane{4}); fyk=str2double(Dane{5}); E=str2double(Dane{6});
Dane1=inputdlg({'Moment zginaj¹cy [kNm]','Si³a tn¹ca [m]'});
M=str2double(Dane1{1}); V=str2double(Dane1{2});
% Pozosta³e sta³e: b i h - odpowiednio wektor szerokoœci i wektor wysokoœci
% przekroju zgodnie z PN-B... 
% ss - minimalne odstêpy pomiêdzy prêtami zbrojenia g³ównego
% dane materia³owe, d³ugosci wektorów szer. i wys.
% ksl - graniczna wzglêdna wysokoœæ strefy œciskanej "ksi"
b=[0.15 0.18 0.2:0.05:1]; h=b; ss=20; lh=length(h); lb=lh; fcd=fck/1.4; fyd=fyk/1.15; Es=E*1000; ksl=0.0035/(0.0035+(fyd/Es));
for j=1:lh
    % petla w której dobierane sa kolejne wysokoœci przekroju
    for i=1:lb
        % petla w której dobierane sa kolejne szerokoœci przekroju
        a1=(c+fis+0.5*fi)/1000; d(i,j)=h(i)-a1; Sc(i,j)=M/(fcd*1000*b(j)*d(i,j)^2);
        if Sc(i,j)<0.5
            ks(i,j)=1-sqrt(1-2*Sc(i,j));
        else
            ks(i,j)=1;
        end
        if ks(i,j)<=ksl
           x(i,j)=ks(i,j)*d(i,j); As(i,j)=(fcd*b(j)*x(i,j)/fyd)*10000; A1=((pi*fi^2)/4)/100; n1(i,j)=ceil(As(i,j)/A1);
           % je¿eli zbrojenie wyjdzie poni¿ej 1 prêta z automatu w
           % przekroju pojawi¹ siê 2 prêty zbrojenia g³ównego zgodnie z
           % normami
           if n1(i,j)==1
               n1(i,j)=2;
           end
           A(i,j)=n1(i,j)*A1; od(i,j)=(b(j)*1000-2*c-2*fis-n1(i,j)*fi)/(n1(i,j)-1);
           if od(i,j)>=ss
               x1(i,j)=(A(i,j)/10000*fyd)/(fcd*b(j)); Mrd(i,j)=fcd*1000*b(j)*x1(i,j)*(d(i,j)-0.5*x1(i,j));
               % podstawowy warunek uzyskania poprawnych wyników 
               %(moment mniejszy od nosnoœci przekroju)
               if M<Mrd(i,j)
                    H(i,j)=h(i); B(i,j)=b(j); N(i,j)=M/Mrd(i,j); AS1(i,j)=A(i,j); LP1(i,j)=n1(i,j); LP2(i,j)=0; LP3(i,j)=LP1(i,j)+LP2(i,j); D(i,j)=d(i,j);
               else
                   H(i,j)=0; B(i,j)=0; N(i,j)=0; AS1(i,j)=0; LP1(i,j)=0; LP2(i,j)=0; LP3(i,j)=LP1(i,j)+LP2(i,j);
               end
           else
               % w tym miejscu nastepuje powtórne obliczenie przekroju z
               % rozmieszczeniem prêtów w dwóch warstwach 
               a1=(c+fis+fi+0.5*ss)/1000; d(i,j)=h(i)-a1; Sc(i,j)=M/(fcd*1000*b(j)*d(i,j)^2);
                if Sc(i,j)<0.5
                    ks(i,j)=1-sqrt(1-2*Sc(i,j));
                else
                    ks(i,j)=1;
                end
                if ks(i,j)<ksl
                    x(i,j)=ks(i,j)*d(i,j); As(i,j)=(fcd*b(j)*x(i,j)/fyd)*10000; A1=((pi*fi^2)/4)/100; n(i,j)=ceil(As(i,j)/A1); n1(i,j)=n(i,j);
                    if n1(i,j)==1
                        n1(i,j)=2;
                    end
                    A(i,j)=n1(i,j)*A1; od(i,j)=(b(j)*1000-2*c-2*fis-n1(i,j)*fi)/(n1(i,j)-1);
                    while od(i,j)<ss
                        n1(i,j)=n1(i,j)-1; od(i,j)=(b(j)*1000-2*c-2*fis-n1(i,j)*fi)/(n1(i,j)-1); n2(i,j)=n(i,j)-n1(i,j);
                    end
                    if n1(i,j)>=n2(i,j)
                        a11=c+fis+0.5*fi; a12=a11+fi+ss; a1=(a11*n1(i,j)+a12*n2(i,j))/n(i,j); d(i,j)=h(i)-a1/1000;
                        x1(i,j)=(A(i,j)/10000*fyd)/(fcd*b(j)); Mrd(i,j)=fcd*1000*b(j)*x1(i,j)*(d(i,j)-0.5*x1(i,j));
               % podstawowy warunek uzyskania poprawnych wyników 
               %(moment mniejszy od nosnoœci przekroju)
                        if M<Mrd(i,j)
                            H(i,j)=h(i); B(i,j)=b(j); N(i,j)=M/Mrd(i,j); AS1(i,j)=A(i,j); LP1(i,j)=n1(i,j); LP2(i,j)=n2(i,j); LP3(i,j)=LP1(i,j)+LP2(i,j); D(i,j)=d(i,j);
                        else
                            H(i,j)=0;
                        end
                    end
                end
           end
        end
    end
end
W1=H(find(H(:))); W2=B(find(H(:))); W3=N(find(H(:))); W4=AS1(find(H(:))); W5=LP1(find(H(:))); W6=LP2(find(H(:))); W7=LP3(find(H(:))); W8=D(find(H(:))); lw=length(W1);
for i=1:lw
    % kolejny warunek dopasowania wyników
    % przekroje z zapasem noœnoœci mniejszym ni¿ 15% 
    %(mo¿na dowolnie modyfikowaæ)
    if W3(i)>0.85
        if W1(i)>1.5*W2(i) && W1(i)<2.5*W2(i)
            H1(i)=W1(i); B1(i)=W2(i); N1(i)=W3(i); AS11(i)=W4(i); LP11(i)=W5(i); LP21(i)=W6(i); LP31(i)=W7(i); D1(i)=W8(i);
        end
    end
end
% pierwsza partia wyników dotycz¹ca prêtów zbrojenia g³ównego
disp('K1 - Indeks; K2 - Wysokoœæ; K3 - Szerokoœæ; K4 - Wysokoœæ uzyteczna; K5 - Moment / Noœnoœæ;') 
disp('K6 - Zbrojenie; K7 - liczba prêtów w 1 warstwie; K8 - liczba prêtów w 2 warstwie; K9 - wszystkie prêty')
Wyniki(:,2)=H1(find(H1)); Wyniki(:,3)=B1(find(H1)); Wyniki(:,5)=N1(find(H1)); Wyniki(:,6)=AS11(find(H1)); Wyniki(:,7)=LP11(find(H1)); Wyniki(:,8)=LP21(find(H1)); Wyniki(:,9)=LP31(find(H1)); Wyniki(:,4)=D1(find(H1));
LW=length(B1(find(H1))); Wyniki(:,1)=1:1:LW;
Wyniki
% algorytm obliczaj¹cy rozstaw strzemion
% algorytm posiada zabezpieczenia w postaci minimalnych i maksymalnych
% rozstawów strzemion, oblicza strzemiona 2 lub 4 ciête
for i=1:LW
    k=1+sqrt(200/(Wyniki(i,4)*1000)); l=ceil(0.3333*Wyniki(i,9)); vmin=0.035*(k^(3/2))*(fck^(0.5));
    if l<2
        l==2;
    end
    Asl=l*A1; ro=(Asl/10000)/(Wyniki(i,3)*Wyniki(i,4)); Vrdc=((0.18/1.4)*k*(100*ro*fck)^(1/3))*Wyniki(i,3)*Wyniki(i,4); Vrds=vmin*Wyniki(i,3)*Wyniki(i,4); Vr(i)=1000*max(Vrdc, Vrds); smax=min([40 0.75*Wyniki(i,4)*100]);
    if Vr(i)>=V
        sk(i)=1; S(i)=smax; nn(i)=2;
    else
        sk(i)=0; nn(i)=2; Asw=nn(i)*fis^2*pi*0.25; s(i)=(Asw*0.9*Wyniki(i,4)*1000*fyd*2/(V*1000)); S(i)=floor((s(i)-10)/10); lw(i)=(V-Vr(i))/(2*V);
        if S(i)<8
            nn(i)=nn(i)*2; Asw=Asw*2; s(i)=(Asw*0.9*Wyniki(i,4)*1000*fyd*2/(V*1000)); S(i)=floor((s(i)-10)/10);
        end
        if S(i)>smax
            S(i)=smax;
        end 
    end
    cott=2; tt=1/cott; v1=0.6*(1-(fck/250)); Vrmax(i)=(1*Wyniki(i,3)*0.9*Wyniki(i,4)*fcd*1000*v1)/(cott+tt); NN(i)=V/Vrmax(i);
end
disp('K1 - Indeks; K2 - Rozstaw; K3 - Si³a / Noœnoœæ; K4 - Iloœæ ramion strzemienia')
Wyniki1(:,1)=1:1:LW; Wyniki1(:,2)=S; Wyniki1(:,3)=NN; Wyniki1(:,4)=nn;
Wyniki1