clc
clear all
close all

%% Entrada de dados
L=2;%m
b=0.1;%m
h=0.3;%m
E=210e9;
alpha=10^(-5);%oC-1
M0=5e3;
theta=10;%oC
dx=0.1; %Discretizacao ao longo da barra
%%
%%Calculos simples
I=b*h^3/12;
EI=E*I;
curvterm_modulo=2*alpha*theta/h;
%%
%% discretizando o dominio e expressoes para esforcos solicitantes
vetorx=[0:dx:2*L];

for cont=1:length(vetorx)
  vetorm(cont)=M0;
  if vetorx(cont)<=L
  vetormbarra(cont)=L-vetorx(cont);
  else
  vetormbarra(cont)=0;
  end
end
Fig=figure('units','normalized','position',[.1 .1 .8 .8],'color','w')
subplot(2,1,1)
plot(vetorx,vetorm,'k')
xlabel('x[m]','fontsize',18)
ylabel('m','fontsize',18)
set(gca,'fontsize',16,'ydir','reverse','ylim',[0 1.1*M0])

subplot(2,1,2)
plot(vetorx,vetormbarra,'k')
xlabel('x[m]','fontsize',18)
ylabel('\delta M','fontsize',18)
set(gca,'fontsize',16,'ydir','reverse','ylim',[0 1.1*L])

%% Integrais
I1=trapz(vetorx,1/EI*vetormbarra.*vetorm);
I2=trapz(vetorx,1/EI*vetormbarra.^2);
I3=-trapz(vetorx,vetormbarra*curvterm_modulo);

% Calculo da incognita hiperestatica
X=-(I1+I3)/I2
X_analitico=3*EI/L^3*(alpha*theta*L^2/h-M0*L^2/(2*EI)) % Solucao Analitica

%% momento fletor real
vetorM_real=vetorm+X*vetormbarra;

Fig=figure('units','normalized','position',[.1 .1 .8 .8],'color','w')
plot(vetorx,vetorM_real,'k')
xlabel('x[m]','fontsize',18)
ylabel('M','fontsize',18)
set(gca,'fontsize',16,'ydir','reverse','ylim',[0 1.1*max(vetorM_real)])