clc
clear all
close all

T=100; %momento torçor na seção transversal, kNm
D=10; % diâmetro da seção transversal, mm
J=pi*D^4/32; % momento de inercia polar da seção transversal, mm4

%% coordenadas polares
vetorr=linspace(0,D/2,10);
vetortheta=linspace(0,2*pi,36);
%[vetorx,vetory]=pol2cart(vetortheta,vetorr);

[thetamalha,rmalha]=meshgrid(vetortheta,vetorr); %cria uma malha no dominio (theta,r)

Tau=T*rmalha/J; %tensão de cisalhamento na torção, na direção tangencial (et)
[Linhas,Colunas]=size(Tau);

%convertendo para coordenadas cartesianas
for cont=1:Linhas
Tyy(cont,:)=Tau(cont,:).*(-sin(vetortheta));
Tzz(cont,:)=Tau(cont,:).*(cos(vetortheta));
end
[ymalha,zmalha]=pol2cart(thetamalha,rmalha); 
%criando o contorno para a seção trasnversal
[ycont,zcont]=pol2cart(vetortheta,D/2);
%Calculando o momento torçor na seção transversal
dMt=Tau.*rmalha; %momento torçor elementar
dMtr=dMt.*rmalha; %incluindo o Jacobiano da transformação de coordenada
for cont=1:Colunas
auxMt(cont)=trapz(vetorr',dMtr(:,cont));
end
Mtcalc=trapz(vetortheta,auxMt)
%plotando a tensão na seçao transversal
Fig=figure('units','normalized','position',[.1 .1 .8 .8],'color','w');
quiver(ymalha,zmalha,Tyy,Tzz);
hold on
plot(ycont,zcont,'k','linewidth',1.2)
axis equal
title('Tensao de cisalhamento na secao transversal')

%plotando o contorno da tensão de cisalhamento na seção transversal
Fig=figure('units','normalized','position',[.1 .1 .8 .8],'color','w');
contour(ymalha,zmalha,Tau);
hold on
plot(ycont,zcont,'k','linewidth',1.2)
axis equal
title('Contorno da tensao de cisalhamento na secao transversal')

%plotando o contorno da tensão de cisalhamento na seção transversal
Fig=figure('units','normalized','position',[.1 .1 .8 .8],'color','w');
surf(ymalha,zmalha,Tau);
axis equal
title('Contorno da tensao de cisalhamento na secao transversal')