% ==========
% Apartado I
% ==========
%
% Ejemplo de la instrucción case
%
method = 'Gauss';
switch lower(method)
   case {'gauss','gauss-jordan'}
      disp('Método directo por sustitución')
   case {'doolittle','crout','cholesky'}
      disp(' Método directo por factorización')
   case {'jacobi','gauss-seidel','relajacion'}
      disp(' Método iterativo de primer orden')
   otherwise
      disp(' Método desconocido')
end
%
% * * * EJERCICIO 1 * * *
% ##################### Inicio del fichero ####################
function R=opera(A,operacion,B)
%
%function function R=opera(A,operacion,B)
%Realiza operaciones con una y otra matriz dependiendo de la operacion
%Ejemplo de uso de la instrucción case
%
%Entrada
%	A,B         matrices
%	operacion   suma('+'),resta('-'),multiplicación('*') o inversión(#')
%Salida
%	R           resultado
%
% Verificacion de los argumentos de entrada
switch nargin
    case 2 % La unica operacion admitida es la inversion
        if operacion~='#';error('Dos argumentos y la operacion no es inversion');end
        [fA,cA]=size(A);
        if fA~=cA;error('La matriz a invertir no es cuadrada');end
        if fA~=rank(A);error('La matriz a invertir no es regular');end
        R=inv(A);
    case 3
        [fA,cA]=size(A);[fB,cB]=size(B);
        switch operacion
            case {'+','-'}
                if fA~=fB|cA~=cB;error('Las matrices no son congruentes');end
                if operacion='-';signo=-1;else;signo=1;end
                R=A+signo*B;
            case '*'
                if cA~=fB;error('Las matrices no son congruentes');end
                R=A*B;
            otherwise
                error('Operación no valida');
        end
    otherwise
        error('Número de argumentos no válido');
end
% ##################### fin del fichero ####################
%
% ============
% Apartado I.I
% ============
%
% * * * EJERCICIO 2 * * *
A=[ 112,84,48,28; 84,81,77,-7; 48,77,146,-83; 28,-7,-83,85];b=[-204;-186;-193;38];
x=A\b
%
% * * * EJERCICIO 3 * * *
clear all
ang=pi*[1/4,1/6];
ne=12;                  % numero de ecuaciones (nudosx2)
ni=13;                  % numero de incognitas (10 Fuerzas+3reacciones: R1X,R1Y,R6Y)
AA=zeros(ne,ni);
% Equilibrio en nudo 1 
AA(1,[1,2,11])=[-cos(ang(1)),-1,1];
AA(2,[1,12])=[-sin(ang(1)),1];
% Equilibrio en nudo 2
AA(3,[1,4,6])=[cos(ang(1)),-1,-cos(ang(2))];
AA(4,[1,3,6])=[sin(ang(1)),1,sin(ang(2))];
% Equilibrio en nudo 3 
AA(5,[2,5,7])=[1,-cos(ang(2)),-1];
AA(6,[3,5])=[-1,-sin(ang(2))];
% Equilibrio en nudo 4 
AA(7,[4,5,9])=[1,cos(ang(2)),-cos(ang(1))];
AA(8,[5,8,9])=[sin(ang(2)),1,sin(ang(1))];
% Equilibrio en nudo 5 
AA(9,[6,7,10])=[cos(ang(2)),1,-1];
AA(10,[6,8])=[-sin(ang(2)),-1];
% Equilibrio en nudo 6 
AA(11,[9,10])=[cos(ang(1)),1];
AA(12,[9,13])=[-sin(ang(1)),1];
% Seleccionamos las correspondientes a los nudos con grados de libertad
A=AA(3:11,1:10);
[faa,caa]=size(AA);
[fa,ca]=size(A);
% 
% Resistencia de las barras
S=diag([4255,6000,6000,3670,3000,3670,6000,6000,4255,3000]);
%S=100000*eye(10);
%
%
L=chol(A*S*A')
% Cargas en los nudos
PP=[ 0,0,0,-1000,0,0,0,-1000,0,0,0,0;       % carga vertical en los nudos superiores
     0,0,0,0,100,-700,0,0,100,-700,100,0;   % carga por vehiculos en el suelo
     0,0,0,0,0,0,-500,0,0,0,-500,0]';       % carga lateral derecha debido al viento
% Seleccionamos las correspondientes a los nudos con grados de libertad
P=PP(3:11,:);
[fpp,cp]=size(PP);
[fp,cp]=size(P);
D=zeros(fp,cp);
F=zeros(fpp,cp);
%
% Solución
for k=1:cp
    X=L'\(-P(:,k));X=L\X; % Desplazamientos
    D(:,k)=X; % Desplazamientos
    F(1:ca,k)=S*A'*X;        % Fuerzas
    R=-AA([1,2,12],1:10)*F(1:ca,k)+PP([1,2,12],k); %Reacciones
    F(ca+1:caa,k)=R;        % Fuerzas
end

    
%
% Representamos las tensiones y desplazamientos máximos
%
n=8;t=linspace(0,1,n);
[X,Y]=meshgrid(t,t);ZT=zeros(n,n);ZD=ZT;
for i=1:n;for j=1:n;ZT(i,j)=max(abs(F*[X(i,j);Y(i,j);1-X(i,j)-Y(i,j)]));end;end;
figure(1);surf(X,Y,ZT);xlabel('\alpha');ylabel('\beta');zlabel('Tensiones');
for i=1:n;for j=1:n;ZD(i,j)=max(abs(D*[X(i,j);Y(i,j);1-X(i,j)-Y(i,j)]));end;end;
figure(2);surf(X,Y,ZD);xlabel('\alpha');ylabel('\beta');zlabel('Desplazamientos');
%figure(1);plot([0,5],[I1(1:7),I2(1:7)]);grid on;legend('i1','i2','i3','i4','i5','i6','i7');
%figure(2);plot([0,5],[I1(8:end),I2(8:end)]);grid on;legend('v1','v2','v3');