format long %formato longo das variáveis para poder comparar resultados utilizando mais casas decimais
%Configuração da viga
comprimento= input('Digite o comprimento da viga : ');
E= input('Digite o módulo de elasticidade da viga : ');
I= input('Digite o momento de inércia z da viga : ');
EI=E*I;
%Quantificações
quant_apoios=input('Digite a quantidade de apoios : ');
quant_cargas_concen=input('Digite a quantidade de carregamentos concentrados : ');
quant_cargas_distri=input('Digite a quantidade de carregamentos distribuídos : ');
quant_rotulas=input('Digite a quantidade de rótulas : ');
quant_momentos=input('Digite a quantidade de momentos concentrados : ');
%Localizações
matriz_loc_apoio=zeros(1,quant_apoios);
for i=1:quant_apoios
    matriz_loc_apoio(i)= input(['Digite o local do apoio ' num2str(i) ': ']);
end
matriz_loc_rot=zeros(1,quant_rotulas);
for i=1:quant_rotulas
    matriz_loc_rot(i)= input(['Digite o local da rótula ' num2str(i) ': ']);
end
matriz_loc_cargas_conc=zeros(1,quant_cargas_concen);
for i=1:quant_cargas_concen
    matriz_loc_cargas_conc(i)= input(['Digite o local da carga concentrada ' num2str(i) ': ']);
end
matriz_loc_momento=zeros(1,quant_momentos);
for i=1:quant_momentos
    matriz_loc_momento(i)= input(['Digite o local da aplicação do momento ' num2str(i) ': ']);
end
matriz_loc_cargas_distri=zeros(1,2,quant_cargas_distri);
for i=1:quant_cargas_distri
    matriz_loc_cargas_distri(1,1,i)= input(['Digite o local de início da carga distribuída ' num2str(i) ': ']);
    matriz_loc_cargas_distri(1,2,i)= input(['Digite o local de fim da carga distribuída ' num2str(i) ': ']);
end
%adicionando nós no carregamento distribuído:
nos_distri=zeros(1, 100,quant_cargas_distri);
for i=1: quant_cargas_distri
    nos_distri(1,:,i)=linspace(matriz_loc_cargas_distri(1,1,i), matriz_loc_cargas_distri(1,2,i), 100);
end
nos_distrif=reshape(nos_distri, 1, 100*quant_cargas_distri);
%ordem crescente dos nós
%a cada par de nós contidos no carregamento distribuído será colocado um nó intermediário
nos=horzcat(matriz_loc_apoio, matriz_loc_cargas_conc, matriz_loc_momento, matriz_loc_rot, nos_distrif);
ordenada_nos = sort(nos);
matriz_loc_nos= unique(ordenada_nos);
quant_eleme=numel(matriz_loc_nos)-1;
%Tipos de apoios:Apoio 1 liberdade em x, apoio 2 liberdade em y, apoio 3
%liberdade de rotação, apoio 4 engaste
matriz_tipo_apoio=zeros(1,quant_apoios);
for i=1:quant_apoios
    matriz_tipo_apoio(i)=input(['Digite o tipo do apoio ' num2str(i) ':']);
end
%Inputs de valores das cargas aplicadas
matriz_vlr_carga_concen=zeros(1,quant_cargas_concen);
for i=1:quant_cargas_concen
    matriz_vlr_carga_concen(i)=input(['Digite o valor da carga concentrada ' num2str(i) ':']);
end
matriz_vlr_momento=zeros(1,quant_momentos);
for i=1:quant_momentos
    matriz_vlr_momento(i)=input(['Digite o valor do momento ' num2str(i) ':']);
end
matriz_vlr_carga_distri=zeros(1,2,quant_cargas_distri);
for i=1:quant_cargas_distri
    matriz_vlr_carga_distri(1,1,i)=input(['Digite o valor do início da carga distribuída ' num2str(i) ':']);
    matriz_vlr_carga_distri(1,2,i)=input(['Digite o valor do fim da carga distribuída ' num2str(i) ':']);
end
%A matriz abaixo foi criada com o intuito de permitir que existam nós
%contidos na extensão do carregamento distribuído, ela representa a equação
%da reta do carregamento distribuído.
func_carg_distri=zeros(1,2,quant_cargas_distri);
for i=1:quant_cargas_distri
    if matriz_vlr_carga_distri(1,1,i)==matriz_vlr_carga_distri(1,2,i)
        func_carg_distri(1,1,i)=0;
        func_carg_distri(1,2,i)=matriz_vlr_carga_distri(1,2,i);
    else
        func_carg_distri(1,1,i)=(matriz_vlr_carga_distri(1,2,i)-matriz_vlr_carga_distri(1,1,i))/(matriz_loc_cargas_distri(1,2,i)-matriz_loc_cargas_distri(1,1,i));
        func_carg_distri(1,2,i)=matriz_vlr_carga_distri(1,1,i);
    end
end
%Definição de elementos
matriz_comprim_elem=zeros(1,quant_eleme);
for i=1:quant_eleme
    matriz_comprim_elem(i)=matriz_loc_nos(i+1)-matriz_loc_nos(i);
end
matriz_loc_elem=zeros(1,2,quant_eleme);
for i=1:numel(matriz_loc_nos)-1
    matriz_loc_elem(1,1,i)=matriz_loc_nos(i);
    matriz_loc_elem(1,2,i)=matriz_loc_nos(i+1);
end
%Matriz dos esforços equivalentes
%Para as cargas distribuídas
no_rot=zeros(1,quant_rotulas);
for i=1:quant_rotulas
    for j=1:numel(matriz_loc_nos)
        if matriz_loc_rot(i)==matriz_loc_nos(j)
            no_rot(i)=j;
        end
    end
end
MEE_distri_1=zeros(4,1,quant_eleme);
for i=1:quant_cargas_distri
    for j=1:quant_eleme
        if matriz_loc_cargas_distri(1,1,i)==matriz_loc_elem(1,1,j)
            contador=j;
            while matriz_loc_cargas_distri(1,2,i)~=matriz_loc_elem(1,2,contador)
                xum=matriz_loc_elem(1,1,contador)-matriz_loc_cargas_distri(1,1,i);
                xdois=matriz_loc_elem(1,2,contador)-matriz_loc_cargas_distri(1,1,i);
                qum=func_carg_distri(1,1,i)*xum+func_carg_distri(1,2,i);
                qdois=func_carg_distri(1,1,i)*xdois+func_carg_distri(1,2,i);
                MEE_distri_1(1,1,contador)=(((7*matriz_comprim_elem(contador)*qum)/20)+((3*matriz_comprim_elem(contador)*qdois)/20));
                MEE_distri_1(2,1,contador)=((((matriz_comprim_elem(contador)^2)*qum)/20)+(((matriz_comprim_elem(contador)^2)*qdois)/30));
                MEE_distri_1(3,1,contador)=(((3*matriz_comprim_elem(contador)*qum)/20)+((7*matriz_comprim_elem(contador)*qdois)/20));
                MEE_distri_1(4,1,contador)=-((((matriz_comprim_elem(contador)^2)*qum)/30)+(((matriz_comprim_elem(contador)^2)*qdois)/20));
                contador=contador+1;
            end
            if matriz_loc_cargas_distri(1,2,i)==matriz_loc_elem(1,2,contador)
                xum=matriz_loc_elem(1,1,contador)-matriz_loc_cargas_distri(1,1,i);
                xdois=matriz_loc_elem(1,2,contador)-matriz_loc_cargas_distri(1,1,i);
                qum=func_carg_distri(1,1,i)*xum+func_carg_distri(1,2,i);
                qdois=func_carg_distri(1,1,i)*xdois+func_carg_distri(1,2,i);
                MEE_distri_1(1,1,contador)=(((7*matriz_comprim_elem(contador)*qum)/20)+((3*matriz_comprim_elem(contador)*qdois)/20));
                MEE_distri_1(2,1,contador)=((((matriz_comprim_elem(contador)^2)*qum)/20)+(((matriz_comprim_elem(contador)^2)*qdois)/30));
                MEE_distri_1(3,1,contador)=(((3*matriz_comprim_elem(contador)*qum)/20)+((7*matriz_comprim_elem(contador)*qdois)/20));
                MEE_distri_1(4,1,contador)=-((((matriz_comprim_elem(contador)^2)*qum)/30)+(((matriz_comprim_elem(contador)^2)*qdois)/20));
            end
        end
    end
end
MEE_distri_2=MEE_distri_1;
%aloação das condições especiais de engastamento perfeito devido à rotula
%no caso 1
for i=1:quant_cargas_distri
    for j=1:quant_eleme
        if matriz_loc_cargas_distri(1,1,i)==matriz_loc_elem(1,1,j)
            contador=j;
            while matriz_loc_cargas_distri(1,2,i)~=matriz_loc_elem(1,2,contador)
                for k=1:numel(matriz_loc_rot)
                    if matriz_loc_elem(1,2,contador)==matriz_loc_rot(k)
                        xum=matriz_loc_elem(1,1,contador)-matriz_loc_cargas_distri(1,1,i);
                        xdois=matriz_loc_elem(1,2,contador)-matriz_loc_cargas_distri(1,1,i);
                        qum=func_carg_distri(1,1,i)*xum+func_carg_distri(1,2,i);
                        qdois=func_carg_distri(1,1,i)*xdois+func_carg_distri(1,2,i);
                        if qum==qdois
                            MEE_distri_1(1,1,contador)=((5*matriz_comprim_elem(contador)*qum)/8);
                            MEE_distri_1(2,1,contador)=(((matriz_comprim_elem(contador)^2)*qum)/8);
                            MEE_distri_1(3,1,contador)=((3*matriz_comprim_elem(contador)*qum)/8);
                            MEE_distri_1(4,1,contador)=0;
                        elseif abs(qum)>abs(qdois)
                            MEE_distri_1(1,1,contador)=((5*matriz_comprim_elem(contador)*qdois)/8);
                            MEE_distri_1(2,1,contador)=(((matriz_comprim_elem(contador)^2)*qdois)/8);
                            MEE_distri_1(3,1,contador)=((3*matriz_comprim_elem(contador)*qdois)/8);
                            MEE_distri_1(4,1,contador)=0;
                            qum=qum-qdois;
                            qdois=0;
                            MEE_distri_1(1,1,contador)=MEE_distri_1(1,1,contador)+((2*matriz_comprim_elem(contador)*qum)/5);
                            MEE_distri_1(2,1,contador)=MEE_distri_1(2,1,contador)+(((matriz_comprim_elem(contador)^2)*qum)/15);
                            MEE_distri_1(3,1,contador)=MEE_distri_1(3,1,contador)+((matriz_comprim_elem(contador)*qum)/10);
                            MEE_distri_1(4,1,contador)=0;
                        elseif abs(qdois)>abs(qum)
                            MEE_distri_1(1,1,contador)=((5*matriz_comprim_elem(contador)*qum)/8);
                            MEE_distri_1(2,1,contador)=(((matriz_comprim_elem(contador)^2)*qum)/8);
                            MEE_distri_1(3,1,contador)=((3*matriz_comprim_elem(contador)*qum)/8);
                            MEE_distri_1(4,1,contador)=0;
                            qdois=qdois-qum;
                            qum=0;
                            MEE_distri_1(1,1,contador)=MEE_distri_1(1,1,contador)+((9*matriz_comprim_elem(contador)*qdois)/40);
                            MEE_distri_1(2,1,contador)=MEE_distri_1(2,1,contador)+(((7*matriz_comprim_elem(contador)^2)*qdois)/120);
                            MEE_distri_1(3,1,contador)=MEE_distri_1(3,1,contador)+((11*matriz_comprim_elem(contador)*qdois)/40);
                            MEE_distri_1(4,1,contador)=0;
                        end
                    end
                end
                contador=contador+1;
            end
            if matriz_loc_cargas_distri(1,2,i)==matriz_loc_elem(1,2,contador)
                for k=1:numel(matriz_loc_rot)
                    if matriz_loc_elem(1,2,contador)==matriz_loc_rot(k)
                        xum=matriz_loc_elem(1,1,contador)-matriz_loc_cargas_distri(1,1,i);
                        xdois=matriz_loc_elem(1,2,contador)-matriz_loc_cargas_distri(1,1,i);
                        qum=func_carg_distri(1,1,i)*xum+func_carg_distri(1,2,i);
                        qdois=func_carg_distri(1,1,i)*xdois+func_carg_distri(1,2,i);
                            if qum==qdois
                                MEE_distri_1(1,1,contador)=((5*matriz_comprim_elem(contador)*qum)/8);
                                MEE_distri_1(2,1,contador)=(((matriz_comprim_elem(contador)^2)*qum)/8);
                                MEE_distri_1(3,1,contador)=((3*matriz_comprim_elem(contador)*qum)/8);
                                MEE_distri_1(4,1,contador)=0;
                            elseif abs(qum)>abs(qdois)
                                MEE_distri_1(1,1,contador)=((5*matriz_comprim_elem(contador)*qdois)/8);
                                MEE_distri_1(2,1,contador)=(((matriz_comprim_elem(contador)^2)*qdois)/8);
                                MEE_distri_1(3,1,contador)=((3*matriz_comprim_elem(contador)*qdois)/8);
                                MEE_distri_1(4,1,contador)=0;
                                qum=qum-qdois;
                                qdois=0;
                                MEE_distri_1(1,1,contador)=MEE_distri_1(1,1,contador)+((2*matriz_comprim_elem(contador)*qum)/5);
                                MEE_distri_1(2,1,contador)=MEE_distri_1(2,1,contador)+(((matriz_comprim_elem(contador)^2)*qum)/15);
                                MEE_distri_1(3,1,contador)=MEE_distri_1(3,1,contador)+((matriz_comprim_elem(contador)*qum)/10);
                                MEE_distri_1(4,1,contador)=0;
                            elseif abs(qdois)>abs(qum)
                                MEE_distri_1(1,1,contador)=((5*matriz_comprim_elem(contador)*qum)/8);
                                MEE_distri_1(2,1,contador)=(((matriz_comprim_elem(contador)^2)*qum)/8);
                                MEE_distri_1(3,1,contador)=((3*matriz_comprim_elem(contador)*qum)/8);
                                MEE_distri_1(4,1,contador)=0;
                                qdois=qdois-qum;
                                qum=0;
                                MEE_distri_1(1,1,contador)=MEE_distri_1(1,1,contador)+((9*matriz_comprim_elem(contador)*qdois)/40);
                                MEE_distri_1(2,1,contador)=MEE_distri_1(2,1,contador)+(((7*matriz_comprim_elem(contador)^2)*qdois)/120);
                                MEE_distri_1(3,1,contador)=MEE_distri_1(3,1,contador)+((11*matriz_comprim_elem(contador)*qdois)/40);
                                MEE_distri_1(4,1,contador)=0;
                            end
                    end
                end
            end
        end
    end
end
%aloação das condições especiais de engastamento perfeito devido à rotula
%no caso 2
for i=1:quant_cargas_distri
    for j=1:quant_eleme
        if matriz_loc_cargas_distri(1,1,i)==matriz_loc_elem(1,1,j)
            contador=j;
            while matriz_loc_cargas_distri(1,2,i)~=matriz_loc_elem(1,2,contador)
                for k=1:numel(matriz_loc_rot)
                    if matriz_loc_elem(1,1,contador+1)==matriz_loc_rot(k)
                        xum=matriz_loc_elem(1,1,contador+1)-matriz_loc_cargas_distri(1,1,i);
                        xdois=matriz_loc_elem(1,2,contador+1)-matriz_loc_cargas_distri(1,1,i);
                        qum=func_carg_distri(1,1,i)*xum+func_carg_distri(1,2,i);
                        qdois=func_carg_distri(1,1,i)*xdois+func_carg_distri(1,2,i);
                        disp(['O valor de qdois é: ' num2str(qdois)]);
                        disp(['O valor de qum é: ' num2str(qum)]);                    
                        if qum==qdois
                            MEE_distri_2(1,1,contador+1)=((3*matriz_comprim_elem(contador+1)*qum)/8);
                            MEE_distri_2(2,1,contador+1)=0;
                            MEE_distri_2(3,1,contador+1)=((5*matriz_comprim_elem(contador+1)*qum)/8);
                            MEE_distri_2(4,1,contador+1)=((-(matriz_comprim_elem(contador+1)^2)*qum)/8);
                        elseif abs(qdois)>abs(qum)
                            disp('entrou aqui');
                            MEE_distri_2(1,1,contador+1)=((3*matriz_comprim_elem(contador+1)*qum)/8);
                            MEE_distri_2(2,1,contador+1)=0;
                            MEE_distri_2(3,1,contador+1)=((5*matriz_comprim_elem(contador+1)*qum)/8);
                            MEE_distri_2(4,1,contador+1)=((-(matriz_comprim_elem(contador+1)^2)*qum)/8);
                            qdois=qdois-qum;
                            qum=0;
                            MEE_distri_2(1,1,contador+1)=MEE_distri_2(1,1,contador+1)+((matriz_comprim_elem(contador+1)*qdois)/10);
                            MEE_distri_2(2,1,contador+1)=0;
                            MEE_distri_2(3,1,contador+1)=MEE_distri_2(3,1,contador+1)+((2*matriz_comprim_elem(contador+1)*qdois)/5);
                            MEE_distri_2(4,1,contador+1)=MEE_distri_2(4,1,contador+1)+(((-matriz_comprim_elem(contador+1)^2)*qdois)/15);                            
                        elseif abs(qum)>abs(qdois)
                            MEE_distri_2(1,1,contador+1)=((3*matriz_comprim_elem(contador+1)*qdois)/8);
                            MEE_distri_2(2,1,contador+1)=0;
                            MEE_distri_2(3,1,contador+1)=((5*matriz_comprim_elem(contador+1)*qdois)/8);
                            MEE_distri_2(4,1,contador+1)=((-(matriz_comprim_elem(contador+1)^2)*qdois)/8);
                            qum=qum-qdois;
                            qdois=0;
                            MEE_distri_2(1,1,contador+1)=MEE_distri_2(1,1,contador+1)+((11*matriz_comprim_elem(contador+1)*qum)/40);
                            MEE_distri_2(2,1,contador+1)=0;
                            MEE_distri_2(3,1,contador+1)=MEE_distri_2(3,1,contador+1)+((9*matriz_comprim_elem(contador+1)*qum)/40);
                            MEE_distri_2(4,1,contador+1)=MEE_distri_2(4,1,contador+1)+(((-7*matriz_comprim_elem(contador+1)^2)*qum)/120);
                        end
                    end
                end
                contador=contador+1;
            end
            if matriz_loc_cargas_distri(1,2,i)==matriz_loc_elem(1,2,contador)
                for k=1:numel(matriz_loc_rot)
                    if matriz_loc_elem(1,1,contador)==matriz_loc_rot(k)
                        xum=matriz_loc_elem(1,1,contador)-matriz_loc_cargas_distri(1,1,i);
                        xdois=matriz_loc_elem(1,2,contador)-matriz_loc_cargas_distri(1,1,i);
                        qum=func_carg_distri(1,1,i)*xum+func_carg_distri(1,2,i);
                        qdois=func_carg_distri(1,1,i)*xdois+func_carg_distri(1,2,i);
                        if qum==qdois
                            MEE_distri_2(1,1,contador)=((3*matriz_comprim_elem(contador)*qum)/8);
                            MEE_distri_2(2,1,contador)=0;
                            MEE_distri_2(3,1,contador)=((5*matriz_comprim_elem(contador)*qum)/8);
                            MEE_distri_2(4,1,contador)=((-(matriz_comprim_elem(contador)^2)*qum)/8);
                        elseif abs(qdois)>abs(qum)
                            MEE_distri_2(1,1,contador)=((3*matriz_comprim_elem(contador)*qum)/8);
                            MEE_distri_2(2,1,contador)=0;
                            MEE_distri_2(3,1,contador)=((5*matriz_comprim_elem(contador)*qum)/8);
                            MEE_distri_2(4,1,contador)=((-(matriz_comprim_elem(contador)^2)*qum)/8);
                            qdois=qdois-qum;
                            qum=0;
                            MEE_distri_2(1,1,contador)=MEE_distri_2(1,1,contador)+((matriz_comprim_elem(contador)*qdois)/10);
                            MEE_distri_2(2,1,contador)=0;
                            MEE_distri_2(3,1,contador)=MEE_distri_2(3,1,contador)+((2*matriz_comprim_elem(contador)*qdois)/5);
                            MEE_distri_2(4,1,contador)=MEE_distri_2(4,1,contador)+(((-matriz_comprim_elem(contador)^2)*qdois)/15);
                        elseif abs(qum)>abs(qdois)
                            MEE_distri_2(1,1,contador)=((3*matriz_comprim_elem(contador)*qdois)/8);
                            MEE_distri_2(2,1,contador)=0;
                            MEE_distri_2(3,1,contador)=((5*matriz_comprim_elem(contador)*qdois)/8);
                            MEE_distri_2(4,1,contador)=((-(matriz_comprim_elem(contador)^2)*qdois)/8);
                            qum=qum-qdois;
                            qdois=0;
                            MEE_distri_2(1,1,contador)=MEE_distri_2(1,1,contador)+((11*matriz_comprim_elem(contador)*qum)/40);
                            MEE_distri_2(2,1,contador)=0;
                            MEE_distri_2(3,1,contador)=MEE_distri_2(3,1,contador)+((9*matriz_comprim_elem(contador)*qum)/40);
                            MEE_distri_2(4,1,contador)=MEE_distri_2(4,1,contador)+(((-7*matriz_comprim_elem(contador)^2)*qum)/120);
                        end
                    end
                end
            end
        end
    end
end                           
%Para cargas concentradas e momentos:
MEE_concen=zeros(4,1,quant_eleme);
for i=1:quant_cargas_concen
    if matriz_loc_cargas_conc(i)==matriz_loc_elem(1,2,quant_eleme)
        MEE_concen(3,1,quant_eleme)=matriz_vlr_carga_concen(i);
    end
    for j=1:quant_eleme
        if matriz_loc_cargas_conc(i)==matriz_loc_elem(1,1,j)
            MEE_concen(1,1,j)=matriz_vlr_carga_concen(i);
        end
    end
end
for i=1:quant_momentos
    if matriz_loc_momento(i)==matriz_loc_elem(1,2,quant_eleme)
        MEE_concen(4,1,quant_eleme)=matriz_vlr_momento(i);
    end
    for j=1:quant_eleme
        if matriz_loc_momento(i)==matriz_loc_elem(1,1,j)
            MEE_concen(2,1,j)=matriz_vlr_momento(i);
        end
    end
end
%somando os vetores de carregamentos equivalentes para os casos 1 e 2:
MEE_soma_1=MEE_concen+MEE_distri_1;
MEE_soma_2=MEE_concen+MEE_distri_2;
MEE_glob_1=zeros(numel(matriz_loc_nos)*2,1);
MEE_glob_2=zeros(numel(matriz_loc_nos)*2,1);
for i=2:quant_eleme
    MEE_glob_1(2*i-1)=MEE_soma_1(3,1,i-1)+MEE_soma_1(1,1,i);
    MEE_glob_1(2*i)=MEE_soma_1(4,1,i-1)+MEE_soma_1(2,1,i);
end
MEE_glob_1(1)=MEE_soma_1(1,1,1);
MEE_glob_1(2)=MEE_soma_1(2,1,1);
MEE_glob_1(numel(matriz_loc_nos)*2-1)=MEE_soma_1(3,1,quant_eleme);
MEE_glob_1(numel(matriz_loc_nos)*2)=MEE_soma_1(4,1,quant_eleme);
for i=2:quant_eleme
    MEE_glob_2(2*i-1)=MEE_soma_2(3,1,i-1)+MEE_soma_2(1,1,i);
    MEE_glob_2(2*i)=MEE_soma_2(4,1,i-1)+MEE_soma_2(2,1,i);
end
MEE_glob_2(1)=MEE_soma_2(1,1,1);
MEE_glob_2(2)=MEE_soma_2(2,1,1);
MEE_glob_2(numel(matriz_loc_nos)*2-1)=MEE_soma_2(3,1,quant_eleme);
MEE_glob_2(numel(matriz_loc_nos)*2)=MEE_soma_2(4,1,quant_eleme);
%geração das matrizes de rigidez local para os casos 1 e 2
matriz_rig_loc_1=zeros(4,4,quant_eleme);
for i=1:quant_eleme
    matriz_rig_loc_1(:,:,i)=(EI/(matriz_comprim_elem(i)^3))*[12 6*matriz_comprim_elem(i) -12 6*matriz_comprim_elem(i); 6*matriz_comprim_elem(i) 4*matriz_comprim_elem(i)^2 -6*matriz_comprim_elem(i) 2*matriz_comprim_elem(i)^2;-12 -6*matriz_comprim_elem(i) 12 -6*matriz_comprim_elem(i); 6*matriz_comprim_elem(i) 2*matriz_comprim_elem(i)^2 -6*matriz_comprim_elem(i) 4*matriz_comprim_elem(i)^2];
end
matriz_rig_loc_2=matriz_rig_loc_1;
%identificação e alocação das rótulas para o caso 1
if quant_rotulas ~= 0
    for i =1:quant_rotulas
        for j=1:numel(matriz_loc_nos)
            if matriz_loc_rot(i)==matriz_loc_nos(j)
                matriz_rig_loc_1(:,:,j-1)=(3*EI/(matriz_comprim_elem(j-1)^3))*[1 matriz_comprim_elem(j-1) -1 0; matriz_comprim_elem(j-1) matriz_comprim_elem(j-1)^2 -matriz_comprim_elem(j-1) 0; -1 -matriz_comprim_elem(j-1) 1 0; 0 0 0 0];
            end
        end
    end
end
%identificação e alocação das rótulas para o caso 2
if quant_rotulas ~= 0
    for i =1:quant_rotulas
        for j=1:numel(matriz_loc_nos)
            if matriz_loc_rot(i)==matriz_loc_nos(j)
                matriz_rig_loc_2(:,:,j)=(3*EI/(matriz_comprim_elem(j)^3))*[1 0 -1 matriz_comprim_elem(j); 0 0 0 0; -1 0 1 -matriz_comprim_elem(j); matriz_comprim_elem(j) 0 -matriz_comprim_elem(j) matriz_comprim_elem(j)^2];
            end
        end
    end
end      
%Criação das matrizes de rigidez global para os casos 1 e 2
MK_1=zeros(numel(matriz_loc_nos)*2,numel(matriz_loc_nos)*2);
for n = 1:quant_eleme   
    for i=1:4    
        for j=1:4   
            if MK_1(i+2*n-2,j+2*n-2)==0    
                MK_1(i+2*n-2,j+2*n-2)= matriz_rig_loc_1(i,j,n);
            else
                MK_1(i+2*n-2,j+2*n-2)=matriz_rig_loc_1(i,j,n)+MK_1(i+2*n-2,j+2*n-2);
            end
        end
    end
end
MK_2=zeros(numel(matriz_loc_nos)*2,numel(matriz_loc_nos)*2);
for n = 1:quant_eleme   
    for i=1:4    
        for j=1:4   
            if MK_2(i+2*n-2,j+2*n-2)==0    
                MK_2(i+2*n-2,j+2*n-2)= matriz_rig_loc_2(i,j,n);
            else
                MK_2(i+2*n-2,j+2*n-2)=MK_2(i+2*n-2,j+2*n-2)+matriz_rig_loc_2(i,j,n);
            end
        end
    end
end
%definição das condições de contorno
%vetor de giros e deslocamentos nodais
VGD = zeros (numel(matriz_loc_nos)*2,1);
for i= 1:numel(matriz_loc_nos)*2
    if mod(i,2)==0
        nome_variavel = ['d', num2str(i)];
        % Criar a variável avaliando a string
        eval([nome_variavel, ' = 1;']);
        VGD(i,1)=eval(nome_variavel);
    else
        nome_variavel = ['o', num2str(i)];
        % Criar a variável avaliando a string
        eval([nome_variavel, ' = 1;']);
        VGD(i,1)=eval(nome_variavel);
    end
end
%condições de cotorno dos apoios
contador=0;
for i=1:quant_apoios
    if matriz_tipo_apoio(i)==1
        contador=contador+1;
    elseif matriz_tipo_apoio(i)==2
        contador=contador+1;
    elseif matriz_tipo_apoio(i)==3
        contador=contador+1;
    elseif matriz_tipo_apoio(i)==4
        contador=contador+2;
    end
end
%alocação das condições de contorno de apoio em VGD
for i=1:numel(matriz_loc_apoio)
    for j=1:numel(matriz_loc_nos)
        if matriz_loc_apoio(i)==matriz_loc_nos(j)
            if matriz_tipo_apoio(i)==1
               VGD(2*j-1)=0;
            elseif matriz_tipo_apoio(i)==2
                VGD(2*j-1)=0;
            elseif matriz_tipo_apoio(i)==3
                VGD(2*j-1)=0;
            elseif matriz_tipo_apoio(i)==4
                VGD(2*j-1)=0;
                VGD(2*j)=0;
            end
        end
    end
end
%Geração dos vetores de coeficientes para montagem do sistema nos casos 1 e 2
%vetor dos coeficientes 1
vet_coef_1=MK_1;
for i=numel(VGD):-1:1 %apagar colunas
    if VGD(i)==0
        vet_coef_1(:,i)=[];
    end
end
for i=numel(VGD):-1:1 %apagar linhas
    if VGD(i)==0
        vet_coef_1(i,:)=[];
    end
end
%vetor dos coeficientes 2
vet_coef_2=MK_2;
for i=numel(VGD):-1:1 %apagar colunas
    if VGD(i)==0
        vet_coef_2(:,i)=[];
    end
end
for i=numel(VGD):-1:1 %apagar linhas
    if VGD(i)==0
        vet_coef_2(i,:)=[];
    end
end
%Geração dos vetores de constantes para montagem do sistema nos casos 1 e 2
%vetor de constantes 1
vet_const_1=MEE_glob_1;
for i=numel(vet_const_1):-1:1 %apagar linhas
    if VGD(i)==0
        vet_const_1(i,:)=[];
    end
end
%vetor de constantes 2
vet_const_2=MEE_glob_2;
for i=numel(vet_const_2):-1:1 %apagar linhas
    if VGD(i)==0
        vet_const_2(i,:)=[];
    end
end
%Resolução do sistema linear de equações para os casos 1 e 2
solu_loc_1 = vet_coef_1 \ vet_const_1;
solu_loc_2 = vet_coef_2 \ vet_const_2;
%alocação das soluções no vetor de giros e deslocamentos
VGD_1=VGD;
VGD_2=VGD;
VGD_1(VGD_1 == 1) = solu_loc_1;%onde tiver 1 em VGD será alterado pelo valor de solu_loc_1
VGD_2(VGD_2 == 1) = solu_loc_2;%onde tiver 1 em VGD será alterado pelo valor de solu_loc_2
%Vetor dos giros e deslocamentos dos elementos, 
VGDE_1=zeros(4,1,quant_eleme);
for i=1:quant_eleme 
    VGDE_1(:,:,i)=VGD_1(2*i-1:(2*i+2),:);
end
VGDE_2=zeros(4,1,quant_eleme);
for i=1:quant_eleme
    VGDE_2(:,:,i)=VGD_2(2*i-1:(2*i+2),:);
end
%Vetor de reações de apoio
solu_geral_1=(MK_1)*VGD_1;
Rap_1=solu_geral_1-MEE_glob_1;
%solução de esforços internos para os casos 1 e 2
solu_elementos_1=zeros(4,1,quant_eleme);
for i=1:quant_eleme
        solu_elementos_1(:,:,i)=matriz_rig_loc_1(:,:,i)*VGDE_1(:,:,i)-MEE_soma_1(:,:,i);
end
%Vetor de giros e deslocamentos global, concatenando os casos 1 e 2
if quant_rotulas~=0
    for i=1:quant_rotulas
        for j=1:numel(matriz_loc_nos)
            if matriz_loc_rot(i)==matriz_loc_nos(j)
                VGDE_1(4,1,j-1)=-VGD_2(2*j);
            end
        end
    end
end
VGDG=zeros(numel(matriz_loc_nos*2+quant_rotulas,1));

%criação do gráfico de esforço cortante
x_cort=zeros((quant_eleme+1)*2,1);
x_cort(1)=0;
x_cort(2)=0.00001;
x_cort(numel(x_cort)-1)=comprimento-0.00001;
x_cort(numel(x_cort))=comprimento;
n=3;
while n<numel(x_cort)-2
    for i=2:numel(matriz_loc_nos)-1
        x_cort(n)=matriz_loc_nos(i)-0.00001;
        x_cort(n+1)=matriz_loc_nos(i)+0.00001;
        n=n+2;
    end
end
y_cort=zeros(numel(x_cort),1);
y_cort(1)=0;
y_cort(numel(y_cort))=0;
for i=2:quant_eleme+1
    y_cort(2*i-2)=round(solu_elementos_1(1,1,i-1), 8);
    y_cort(2*i-1)=-round(solu_elementos_1(3,1,i-1),8);
end
if quant_cargas_concen~=0
    for i=1:quant_cargas_concen
        for j=1:numel(matriz_loc_nos)
            if matriz_loc_cargas_conc(i)==matriz_loc_nos(j)
                y_cort(j*2)=y_cort(j*2)+matriz_vlr_carga_concen(i);
            end
        end
    end
end
%Adicionando o Gráfico de esforço cortante
subplot(1,2,1)
plot(x_cort,y_cort,'LineWidth',2)
title('Diagrama de Esforço Cortante (kN)')
xlabel('X (m)');
ylabel('V (kN)');
xlim([-1 comprimento+1]);
% Para destacar os eixos x e y com linhas adicionais
hold on;
plot([min(x_cort), max(x_cort)], [0, 0], 'k--'); 
plot([0, 0], [min(y_cort), max(y_cort)], 'k--'); 
x_momen=zeros((quant_eleme+1)*2,1);
x_momen(1)=0;
x_momen(2)=0.00001;
x_momen(numel(x_momen)-1)=comprimento-0.00001;
x_momen(numel(x_momen))=comprimento;
n=3;
while n<numel(x_momen)-2
    for i=2:numel(matriz_loc_nos)-1
        x_momen(n)=matriz_loc_nos(i)-0.00001;
        x_momen(n+1)=matriz_loc_nos(i)+0.00001;
        n=n+2;
    end
end
y_momen=zeros(numel(x_momen),1);
y_momen(1)=0;
y_momen(numel(y_momen))=0;
for i=2:quant_eleme+1
    y_momen(2*i-2)=round(solu_elementos_1(2,1,i-1),8);
    y_momen(2*i-1)=-round(solu_elementos_1(4,1,i-1),8);
end
if quant_momentos~=0
    for i=1:quant_momentos
        for j=1:numel(matriz_loc_nos)
            if matriz_loc_momento(i)==matriz_loc_nos(j)
                y_momen(j*2)=y_momen(j*2)+matriz_vlr_momento(i);
            end
        end
    end
end
y_momen=y_momen*-1;
%Adicionando o gráfico de momento fletor
subplot(1,2,2)
plot(x_momen,y_momen,'LineWidth',2,'Color',[1 0 0])
set(gca, 'YDir', 'reverse');
title('Diagrama de Momento Fletor (kN.m)')
xlabel('X (m)');
ylabel('M (kN.m)');
xlim([-1 comprimento+1]);
% Para destacar os eixos x e y com linhas adicionais
hold on;
plot([min(x_momen), max(x_momen)], [0, 0], 'k--'); 
plot([0, 0], [min(y_momen), max(y_momen)], 'k--');
%Exportando os dados para o excel
%Esforços internos:
tabelacort=table(x_cort,y_cort, 'VariableNames',{'x(m)', 'V (kN)'});
writetable(tabelacort,'Resultados.xlsx', 'Sheet','Esforços internos');
tabelamoment=table(x_momen,y_momen, 'VariableNames',{'x(m)', 'M (kN.m)'});
writetable(tabelamoment,'Resultados.xlsx', 'Sheet','Esforços internos','Range','D1');
%Reação de Apoio
xRap=repelem(matriz_loc_nos',2);
tabelaRap=table(xRap,Rap_1, 'VariableNames',{'x(m)', 'Força e momento intercalados'});
writetable(tabelaRap,'Resultados.xlsx', 'Sheet','Reações de apoio');
%Giros e deslocamentos
x_VGDE = repmat({'deslocamento'; 'giro'}, 2, 1);
tabelaVGD=table(x_VGDE,VGDE_1, 'VariableNames',{ 'Deslocamentos e giros','Elemento'});
writetable(tabelaVGD,'Resultados.xlsx', 'Sheet','Deslocamentos');