clc
clear all
close all

%% Dimensoes da viga [mm]
b=10;
h=30;
L=100;

%% Posicoes ao longo da altura
y=linspace(-h/2,h/2,30);

%% Posicoes ao longo da viga
x=linspace(0,L,1000); 
%% Propriedades
I=b*h^3/12; 
Ms=abs(b*(h/2-y).*(y+1/2*(h/2-y)));

Fig1=figure('units','normalized','position',[.1 .1 .8 .8],'color','w');
plot(Ms,y,'k');
ylabel('y[mm]','fontsize',12);
xlabel('M_s[mm^3]','fontsize',12);
set(gca,'fontsize',12,'ydir','reverse')
grid

%% Carrgamento uniformemente distribuido [N/mm]
q=10;

%% Esforcos solicitantes

V=q*L/2-q*x;
M=-1/2*q*x.^2+q*L*x/2;



%% Plotagem
Fig2=figure('units','normalized','position',[.1 .1 .8 .8],'color','w');
subplot(2,1,1)
plot(x,V,'k');
xlabel('x[mm]','fontsize',12);
ylabel('V(x)[N]','fontsize',12);
set(gca,'fontsize',12)
grid
subplot(2,1,2)
plot(x,M,'k');
xlabel('x[mm]','fontsize',12);
ylabel('M(x)[Nmm]','fontsize',12);
set(gca,'fontsize',12,'ydir','reverse')
grid


MatrizSigma1=zeros(length(x),length(y));
MatrizSigma2=zeros(length(x),length(y));
Matrizv1x=zeros(length(x),length(y));
Matrizv1y=zeros(length(x),length(y));
Matrizv2x=zeros(length(x),length(y));
Matrizv2y=zeros(length(x),length(y));
MatrizSigma=zeros(length(x),length(y));
MatrizTau=zeros(length(x),length(y));
Matrizrho1=zeros(length(x),length(y));
Matrizrho2=zeros(length(x),length(y));


for contx=1:length(x)
    for conty=1:length(y)
        
        Sigma=M(contx)*y(conty)/I;
        Tau=V(contx)*Ms(conty)/(I*b);
        T=[Sigma Tau;Tau 0];
        rho=T*[1;0]; %tensao na seção de normal unitaria +e1
        Matrizrho1(contx,conty)=rho(1);
        Matrizrho2(contx,conty)=rho(2);
        [AutoVetor,Diag]=eig(T);
        [Sigma1,ind1]=max([Diag(1,1) Diag(2,2)]);
        [Sigma2,ind2]=min([Diag(1,1) Diag(2,2)]);
        
        MatrizSigma1(contx,conty)=Sigma1;
        MatrizSigma2(contx,conty)=Sigma2;
        MatrizSigma(contx,conty)=Sigma;
        MatrizTau(contx,conty)=Tau;
        
        Matrizv11(contx,conty)=AutoVetor(1,ind1)/norm(AutoVetor(:,ind1));
        Matrizv12(contx,conty)=AutoVetor(2,ind1)/norm(AutoVetor(:,ind1));
        Matrizv21(contx,conty)=AutoVetor(1,ind2)/norm(AutoVetor(:,ind2));
        Matrizv22(contx,conty)=AutoVetor(2,ind2)/norm(AutoVetor(:,ind2));
        
    end
end

[malhax,malhay]=meshgrid(x,y);


Fig3=figure('units','normalized','position',[.1 .1 .8 .8],'color','w');
subplot(2,1,2)
h=pcolor(x,y,MatrizSigma1');
set(h,'edgecolor','none')
hc=colorbar('eastoutside')
set(gca,'fontsize',12)
set(get(hc,'ylabel'),'string','\sigma_1 [MPa]','fontsize',20)
xlabel('x[mm]','fontsize',12);
ylabel('y[mm]','fontsize',12);
set(gca,'fontsize',12,'ydir','reverse')

subplot(2,1,1)
h=pcolor(x,y,MatrizSigma2');
set(h,'edgecolor','none')
hc=colorbar('eastoutside')
set(gca,'fontsize',12)
set(get(hc,'ylabel'),'string','\sigma_2 [MPa]','fontsize',20)
xlabel('x[mm]','fontsize',12);
ylabel('y[mm]','fontsize',12);
set(gca,'fontsize',12,'ydir','reverse')

%%
Fig4=figure('units','normalized','position',[.1 .1 .8 .8],'color','w');
subplot(2,1,1)
h=pcolor(x,y,MatrizSigma');
set(h,'edgecolor','none')
hc=colorbar('eastoutside')
set(gca,'fontsize',12)
set(get(hc,'ylabel'),'string','\sigma [MPa]','fontsize',20)
xlabel('x[mm]','fontsize',12);
ylabel('y[mm]','fontsize',12);
set(gca,'fontsize',12,'ydir','reverse')

subplot(2,1,2)
h=pcolor(x,y,MatrizTau');
set(h,'edgecolor','none')
hc=colorbar('eastoutside')
set(gca,'fontsize',12)
set(get(hc,'ylabel'),'string','\tau [MPa]','fontsize',20)
xlabel('x[mm]','fontsize',12);
ylabel('y[mm]','fontsize',12);
set(gca,'fontsize',12,'ydir','reverse')


%% Mostrando tensões ao longo da seção transversal
secao1=min(find(x>L/4)) %posicao do vetor x onde tenho x=L/4
xauxsecao=zeros(1,length(y));
rho1_L4=Matrizrho1(secao1,:);
rho2_L4=Matrizrho2(secao1,:);
h11_L4=Matrizv11(secao1,:);
h12_L4=Matrizv12(secao1,:);
Sigma1_L4=MatrizSigma1(secao1,:);


Fig5=figure('units','normalized','position',[.1 .1 .8 .8],'color','w');
quiver(xauxsecao,y,rho1_L4,rho2_L4,'k');
ylabel('y[mm]','fontsize',12);
xlabel('-','fontsize',12);
title('Secao x=L/4; distribuicao da tensao \rho','fontsize',18)
set(gca,'fontsize',12,'ydir','reverse')
grid

Fig7=figure('units','normalized','position',[.1 .1 .8 .8],'color','w');
quiver(xauxsecao,y,Sigma1_L4.*h11_L4,Sigma1_L4.*h12_L4,'k');
ylabel('y[mm]','fontsize',12);
xlabel('-','fontsize',12);
title('Secao x=L/4; \sigma_1\bf{h_1}','fontsize',18)
set(gca,'fontsize',12,'ydir','reverse')
grid

secao2=min(find(x>L/2)) %posicao do vetor x onde tenho x=L/2
xauxsecao=zeros(1,length(y));
rho1_L2=Matrizrho1(secao2,:);
rho2_L2=Matrizrho2(secao2,:);
h11_L2=Matrizv11(secao2,:);
h12_L2=Matrizv12(secao2,:);
Sigma1_L2=MatrizSigma1(secao2,:);

Fig8=figure('units','normalized','position',[.1 .1 .8 .8],'color','w');
quiver(xauxsecao,y,rho1_L2,rho2_L2,'k');
ylabel('y[mm]','fontsize',12);
xlabel('-','fontsize',12);
title('Secao x=L/2; distribuicao da tensao \rho','fontsize',18)
set(gca,'fontsize',12,'ydir','reverse')
grid

Fig9=figure('units','normalized','position',[.1 .1 .8 .8],'color','w');
quiver(xauxsecao,y,Sigma1_L2.*h11_L2,Sigma1_L2.*h12_L2,'k');
ylabel('y[mm]','fontsize',12);
xlabel('-','fontsize',12);
title('Secao x=L/2;  \sigma_1\bf{h_1}','fontsize',18)
set(gca,'fontsize',12,'ydir','reverse')
grid