clc, clear,
format short

% ---- Pow³oka walcowa ----
% ---- Wprowadzenie danych ----
Dane1=inputdlg({'Wysokoœæ pow³oki [m]','Œrednica pow³oki [m]','Iloœæ pierœcieni p³aszcza'});
L=str2double(Dane1{1}); D=str2double(Dane1{2}); n=str2double(Dane1{3});
a=D/2; l=L/n; u1=0:0.001:L; LL=length(u1); uu1=0:0.001:l; ll=length(uu1)-1; hw=L; gs=78.5; 
gc=10; v=0.3; E=210000000; m1=-sqrt(1/(2*a)); n1=m1; g=a^2; g11=1; g22=a^2; beta22=1/a; 

for i=1:n
    Dane=inputdlg({'WprowadŸ gruboœci pierœcieni [mm] (od najgrubszego)'});
    tpl(i)=str2double(Dane{1})/1000; h(i)=0.5*tpl(i);
    
    % ---- Stan b³onowy ----
    for j=i*ll-ll+1:(i+1)*ll-ll
            N11B(j)=2*h(i)*gs*(u1(j)-L); N22B(j)=a*gc*(hw-u1(j));
            M12B(j)=-2*h(i)^3*gs*(u1(j)-L)/(3*a); M21B(j)=-h(i)^2*gc*(hw-u1(j))/3;
            Q1B(j)=2*gs*h(i)^3/(3*a); W1B(j)=(h(i)*gs+0.5*v*a*gc)*(u1(j)^2)-...
                (2*h(i)*gs*L+v*a*gc*hw)*u1(j)/(2*E*h(i));
            W3B(j)=(a^2*gc*(hw-u1(j))-a*v*2*h(i)*gs*(u1(j)-L))/(2*E*h(i));

            % ---- stan zgiêciowy ----
            eps=1/a; om4(i)=3*(1-v^2)/h(i)^2; om2(i)=om4(i)^(1/2); om(i)=om4(i)^(1/4);
            S12(j)=exp(m1*u1(j)*om(i))*sin(om(i)*n1*u1(j)); S22(j)=exp(m1*u1(j)*om(i))*cos(om(i)*n1*u1(j));
            C222(j)=(a^2*gc*hw+L*a*v*2*h(i)*gs)/g; C122(j)=((a^2*gc+a*v*2*h(i)*gs)/(om(i)*g*n1))+C222(j);
            C212(j)=-C122(j); C112(j)=C222(j);
            
            N22Z(j)=-a*(C212(j)*S12(j)+C222(j)*S22(j)); M12Z(j)=(g*eps*g11/om2(i))*(C112(j)*S12(j)+C122(j)*S22(j));
            M21Z(j)=-(a^2)*(((1-v)*beta22-g*eps*g22)*(C112(j)*S12(j)+C122(j)*S22(j)))/om2(i);
            Q1Z(j)=(eps*g*g11*n1/om(i))/((C212(j)+C112(j))*S12(j)+(C222(j)+C122(j))*S22(j));
            W1Z(j)=v*Q1Z(j)/(E*2*h(i)*eps); W3Z(j)=(g/(E*2*h(i)))*(C212(j)*S12(j)+C222(j)*S22(j));
    end
end
i=n; j=LL;

% ---- Stan b³onowy CD ----
N11B(j)=2*0.5*tpl(i)*gs*(u1(j)-L); N22B(j)=a*gc*(u1(j)-hw);  M12B(j)=-(2*gs*(u1(j)-L)*(0.5*tpl(i))^3)/(3*a);
M21B(j)=-(a^2)*(((0.5*tpl(i))^2)/(3*a))*gc*(u1(j)-hw); Q1B(j)=(2*gs*(0.5*tpl(i))^3)/(3*a); W3B(j)=(-1/(2*E*(0.5*tpl(i))))*(a^2*gc*(u1(j)-hw)-a*v*2*(0.5*tpl(i))*gs*(u1(j)-L));
W1B(j)=(((0.5*tpl(i))*gs+0.5*v*a*gc)*(u1(j)^2)-(2*(0.5*tpl(i))*gs*L+v*a*gc*hw)*u1(j))/(2*E*(0.5*tpl(i)));

% ---- Stan zgiêciowy CD ----
N11B(j)=2*h(i)*gs*(u1(j)-L); N22B(j)=a*gc*(hw-u1(j)); M12B(j)=-2*h(i)^3*gs*(u1(j)-L)/(3*a);
M21B(j)=-h(i)^2*gc*(hw-u1(j))/3; Q1B(j)=2*gs*h(i)^3/(3*a); W1B(j)=(h(i)*gs+0.5*v*a*gc)*(u1(j)^2)-...
    (2*h(i)*gs*L+v*a*gc*hw)*u1(j)/(2*E*h(i)); W3B(j)=(a^2*gc*(hw-u1(j))-a*v*2*h(i)*gs*(u1(j)-L))/(2*E*h(i));

% ---- stan zgiêciowy CD----
eps=1/a; om4(i)=3*(1-v^2)/h(i)^2; om2(i)=om4(i)^(1/2); om(i)=om4(i)^(1/4);
S12(j)=exp(m1*u1(j)*om(i))*sin(om(i)*n1*u1(j)); S22(j)=exp(m1*u1(j)*om(i))*cos(om(i)*n1*u1(j));
C222(j)=(a^2*gc*hw+L*a*v*2*h(i)*gs)/g; C122(j)=((a^2*gc+a*v*2*h(i)*gs)/(om(i)*g*n1))+C222(j);
C212(j)=-C122(j); C112(j)=C222(j);
            
N22Z(j)=-a*(C212(j)*S12(j)+C222(j)*S22(j)); M12Z(j)=(g*eps*g11/om2(i))*(C112(j)*S12(j)+C122(j)*S22(j));
M21Z(j)=-(a^2)*(((1-v)*beta22-g*eps*g22)*(C112(j)*S12(j)+C122(j)*S22(j)))/om2(i);
Q1Z(j)=(eps*g*g11*n1/om(i))/((C212(j)+C112(j))*S12(j)+(C222(j)+C122(j))*S22(j));
W1Z(j)=v*Q1Z(j)/(E*2*h(i)*eps); W3Z(j)=(g/(E*2*h(i)))*(C212(j)*S12(j)+C222(j)*S22(j));

% ---- Koñcowe wartoœci si³ i przemieszczeñ ----
N11=N11B; N22=N22B+N22Z; M12=M12B+M12Z; M21=M21B+M21Z;
Q1=Q1B+Q1Z; W1=W1B+W1Z; W3=W3B+W3Z;

% ---- Maksymalne wartoœci si³ w obrêbie pojedyñczego pierœcienia ----
for i=1:n
    for j=i*ll-ll+1:(i+1)*ll-ll
        NN11(i,j)=N11B(j); NN22(i,j)=N22B(j); QQ1(i,j)=Q1B(j);
    end
end
i=n; j=LL;
NN11(i,j)=N11B(j); NN22(i,j)=N22B(j); QQ1(i,j)=Q1B(j);
for i=1:n
    N1(i)=abs(min(NN11(i,:))); N2(i)=max(NN22(i,:)); Q(i)=abs(max(QQ1(i,:)));
end

% ---- Wykresy stanu b³onowego ----
figure(1)
subplot(1,2,1)
plot(N11B,u1)
title('Wykres si³ po³udnikowych - N11 - stan b³onowy')
xlabel('Wartoœci si³ w kN/m')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
subplot(1,2,2)
plot(N22B,u1)
title('Wykres si³ obwodowych - N22 - stan b³onowy')
xlabel('Wartoœci si³ w kN/m')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
figure(2)
subplot(1,2,1)
plot(M12B,u1)
title('Wykres momentów zginaj¹cych - M12 - stan b³onowy')
xlabel('Wartoœci momentów w kNm/m')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
subplot(1,2,2)
plot(M21B,u1)
title('Wykres momentów zginaj¹cych - M21 - stan b³onowy')
xlabel('Wartoœci momentów w kNm/m')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
figure(3)
subplot(1,3,1)
plot(Q1B,u1)
title('Wykres si³ tn¹cych - Q - stan b³onowy')
xlabel('Wartoœci si³ w kN/m')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
subplot(1,3,2)
plot(W1B,u1)
title('Wykres przemieszczeñ - w3 - stan b³onowy')
xlabel('Wartoœci przemieszczeñ w metrach')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
subplot(1,3,3)
plot(W3B,u1)
title('Przemieszczeñ - w1 - stan b³onowy')
xlabel('Wartoœci przemieszczeñ w metrach')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
% ---- Wykresy stanu zgiêciowego ----
figure(4)
subplot(1,2,1)
plot(Q1Z,u1)
title('Wykres si³ tn¹cych - Q1 - stan zgiêciowy')
xlabel('Wartoœci si³ w kN/m')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
subplot(1,2,2)
plot(N22Z,u1)
title('Wykres si³ obwodowych - N22 - stan zgiêciowy')
xlabel('Wartoœci si³ w kN/m')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
figure(5)
subplot(1,2,1)
plot(M12Z,u1)
title('Wykres momentów zginaj¹cych - M12 - stan zgiêciowy')
xlabel('Wartoœci momentów w kNm/m')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
subplot(1,2,2)
plot(M21Z,u1)
title('Wykres momentów zginaj¹cych - M21 - stan zgiêciowy')
xlabel('Wartoœci momentów w kNm/m')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
figure(6)
subplot(1,2,1)
plot(W1Z,u1)
title('Wykres przemieszczeñ - w3 - stan zgiêciowy')
xlabel('Wartoœci przemieszczeñ w metrach')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
subplot(1,2,2)
plot(W3Z,u1)
title('Przemieszczeñ - w1 - stan zgiêciowy')
xlabel('Wartoœci przemieszczeñ w metrach')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
% ---- Wykresy koñcowych wartoœci si³ i przemieszczeñ ----
figure(7)
subplot(1,2,1)
plot(N11,u1)
title('Wykres si³ po³udnikowych - N11')
xlabel('Wartoœci si³ w kN/m')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
subplot(1,2,2)
plot(N22,u1)
title('Wykres si³ obwodowych - N22')
xlabel('Wartoœci si³ w kN/m')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
figure(8)
subplot(1,2,1)
plot(M12,u1)
title('Wykres momentów zginaj¹cych - M12')
xlabel('Wartoœci momentów w kNm/m')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
subplot(1,2,2)
plot(M21,u1)
title('Wykres momentów zginaj¹cych - M21')
xlabel('Wartoœci momentów w kNm/m')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
figure(9)
subplot(1,3,1)
plot(Q1,u1)
title('Wykres si³ tn¹cych - Q')
xlabel('Wartoœci si³ w kN/m')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
subplot(1,3,2)
plot(W1,u1)
title('Wykres przemieszczeñ - w3')
xlabel('Wartoœci przemieszczeñ w metrach')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
subplot(1,3,3)
plot(W3,u1)
title('Przemieszczeñ - w1')
xlabel('Wartoœci przemieszczeñ w metrach')
ylabel('Wysokoœæ zbiornika w metrach')
grid on
% ---- Noœnoœæ pow³oki ----

gm0=1.1; lambx0=0.2; lambtet0=0.4; lambt0=lambtet0; Ctet=1.5; beta=0.6; eta=1;
S=menu('Klasa stali konstrukcyjnej','S235','S275','S355','S420','S460');
switch S; 
    case 1; fy=235; fyk=fy/gm0; case 2; fy=275; fyk=fy/gm0;
    case 3; fy=355; fyk=fy/gm0; case 4; fy=460; fyk=fy/gm0; case 5; fy=460; fyk=fy/gm0;
end
KL=menu('Klasa jakoœci wytworzenia','A - Najwy¿sza','B - Wysoka','C - Normalna');
switch KL;
    case 1; QQ=40; alfatet=0.75; alfat=alfatet; case 2; QQ=25; alfatet=0.65; alfat=alfatet;
    case 3; QQ=16; alfatet=0.50; alfat=alfatet; 
end
for i=1:n
    omega(i)=L/sqrt(a*tpl(i));
    if omega(i)<1.7;
        Cx(i)=1.36-(1.83/omega(i))+(2.07/omega(i)^2);
    elseif 1.7<=omega(i)<=0.5*a/tpl(i);
        Cx(i)=1;
    elseif omega(i)>0.5*a/tpl(i);
        Cxb=6; Cx(i)=max(1+(0.2/Cxb)*(1-2*omega(i)*tpl(i)/a),0.6);
    end
    if omega(i)/Ctet<20;
        Cteti(i)=1.5+(10/omega(i)^2)-(5/omega(i)^3);
        sigtetrcr(i)=(0.92*E*Cteti(i)*tpl(i)/(omega(i)*a))/1000;
    elseif 20<=omega(i)/Ctet<=1.63*a/tpl(i);
        Ctet=1.5;  sigtetrcr(i)=(0.92*E*Ctet*tpl(i)/(omega(i)*a))/1000;
    elseif omega(i)/Ctet>1.63*a/tpl(i);
        Ctet=1.5;  sigtetrcr(i)=(E*(tpl(i)/a)^2*(0.275+2.03*(Ctet*a/(omega(i)*tpl(i)))^4))/1000;
    end
    if omega(i)<10;
        Ct(i)=sqrt(1+42/omega(i)^3);
    elseif 10<=omega(i)<=8.7*a/tpl(i);
        Ct(i)=1;
    elseif omega(i)>0.5*a/tpl(i);
        Ct(i)=(1/3)*sqrt(omega(i)*t/a);
    end
    alfax(i)=0.62/(1+1.91*((1/QQ)*sqrt(a/tpl(i)))^1.44);
    Txtetrcr(i)=(0.75*E*Ct(i)*sqrt(1/omega(i))*(tpl(i)/a))/1000;
    sigxrcr(i)=(0.605*E*Cx(i)*tpl(i)/a)/1000;
    lambxp(i)=sqrt(alfax(i)/(1-beta)); lambtetp=sqrt(alfatet/(1-beta));
    lambtp=sqrt(alfat/(1-beta));
    lambx(i)=sqrt(fyk/sigxrcr(i)); lambtet(i)=sqrt(fyk/sigtetrcr(i)); lambt(i)=sqrt(fyk/sqrt(3)*Txtetrcr(i));
    if lambx(i)<=lambx0
        Xx(i)=1;
    elseif (lambx(i)>lambx0 && lambx(i)<lambxp(i))
        Xx(i)=abs(1-(beta*((lambx(i)-lambx0)/(lambxp(i)-lambx0))^eta));
    elseif lambx(i)>=lambxp(i)
        Xx(i)=abs(alfax(i)/lambx(i)^2);
    end
        if lambt(i)<=lambt0
        Xt(i)=1;
    elseif (lambt0<lambt(i) && lambtp>lambt(i))
        Xt(i)=abs(1-beta*((lambt(i)-lambt0)/(lambtp-lambt0))^eta);
    elseif lambt(i)>=lambtp
        Xt(i)=abs(alfat/(lambt(i))^2);
    end
    if lambtet(i)<=lambtet0
        Xtet(i)=1;
    elseif (lambtet0<lambtet(i) && lambtetp>lambtet(i))
        Xtet(i)=(1-beta*((lambtet(i)-lambtet0)/(lambtetp-lambtet0))^eta);
    elseif lambtet(i)>=lambtetp
        Xtet(i)=(alfatet/lambtet(i)^2);
    end
    sigxrd(i)=Xx(i)*fyk/gm0; sigtetetrd(i)=Xtet(i)*fyk/gm0; Txtetrd(i)=Xt(i)*fyk/(3^(1/3)*gm0);
    sigxed(i)=(N1(i)/tpl(i))/1000; sigteted(i)=(N2(i)/tpl(i))/1000; Txteted(i)=(Q(i)/tpl(i))/1000;
    WW1(i,1)=N1(i); WW1(i,2)=sigteted(i); WW1(i,3)=Xx(i); WW1(i,4)=sigxrd(i); WW1(i,5)=sigxed(i)/sigxrd(i);
    WW2(i,1)=N2(i); WW2(i,2)=sigteted(i); WW2(i,3)=Xtet(i); WW2(i,4)=sigtetetrd(i); WW2(i,5)=sigteted(i)/sigtetetrd(i);
    WW3(i,1)=Q(i); WW3(i,2)=Txteted(i); WW3(i,3)=Xt(i); WW3(i,4)=Txtetrd(i); WW3(i,5)=Txteted(i)/Txtetrd(i);
    
    % ---- Z³o¿ony stan naprê¿enia ----
    kx(i)=1.25+0.75*Xx(i); ktet(i)=1.25+0.75*Xtet(i); kt(i)=1.75+0.25*Xt(i); ki(i)=(Xx(i)*Xtet(i))^2;
    Z(i)=WW1(i,5)^kx(i)-ki(i)*WW1(i,5)*WW2(i,5)+WW2(i,5)^ktet(i)+WW3(i,5)^kt(i);
end
disp('K1 - Si³a; K2 - Naprê¿enia; K3 - Wp³. redukcyjny; K4 - Noœnoœæ; K5 - Wytê¿enie')
disp('Œciskanie po³udnikowe'), WW1,
disp('Œciskanie obwodowe'), WW2,
disp('Œcinanie'), WW3
disp('z³o¿ony stan naprê¿enia'), Z,