function [] = show_eigmode(kl,km,kr)

% function show_eigmode
% 
% Purpose: Simulates the first two eigenmodes
%          of the bridge in Exercise 1 of Exercise sheet 4
%          of the lecture Geometric Approximation Theory
% Input:  kl,km,kr
%         Stiffness of the springs on the left part the middle spring and
%         the right part of the bridge, respectively 
%
% Output: none
%

%Simulation time
sim_time=3;

%number of timesteps
nz=100;

%stable configuration
A=[[0 1 2 3 4;0 0 0 0 0] [0.8 1.6 2.4 3.2;-.5 -.5 -.5 -.5] [1.3 2 2.7;-1 -1 -1] [1.75 2.25;-1.5 -1.5] [2;-1.9]];
A=[A(1,:) A(1,:)+5;A(2,:) A(2,:)];

%Edges in one part of the bridge
s=zeros(15);
for i=1:14
    if i~=5 && i~=9 && i~=12 && i~=14
        s(i,i+1)=1;
    end
end
for i=1:4
    s(i,i+5)=1;
end
for i=2:5
    s(i,i+4)=1;
end
for i=6:8
    s(i,i+4)=1;
end
for i=7:9
    s(i,i+3)=1;
end
s(10,13)=1;
s(11,13)=1;
s(11,14)=1;
s(12,14)=1;
s(13,15)=1;
s(14,15)=1;
s=s+s';

%Create Stiffness matrix
K=[kl*s zeros(15);zeros(15) kr*s];
K(5,16)=km;
K(16,5)=km;

%compute eigenmode
V=geteigenmode(A,K);

for eigm=1:2   %1st and 2nd eigenmode
    vi=V(:,eigm);
%scale vi
vi=vi/2/max(abs(vi));

%add fixed points
vn=zeros(2,30);
tild=[2:14 16:19 21:29];
vn(:,tild)=reshape(vi,[2,26]);

figure
for k=1:nz
    clf;
    An=A+cos(4*pi*k/nz)*vn;
    for i=1:30
        for j=1:30
            if K(i,j)>0
                plot([An(1,i),An(1,j)],[An(2,i),An(2,j)],'k');
                hold on;
            end
        end
    end
    for i=[1 15 20 30]
        plot(A(1,i),A(2,i),'o','Color','k','MarkerFaceColor','k');
        hold on;
    end
    axis equal
    axis([0 9 -2 1]);
    title(['Eigenmode ' num2str(eigm)]);
    pause(sim_time/nz);
end

end
end

function [v] = geteigenmode(A,K)
%computes eigenmode of a springsystem corresponding to the smallest
%eigenvalue
%A:dxn position of masses
%K:nxn spring coefficients

tild=[2:14 16:19 21:29];
[d,n]=size(A);
S=zeros(d*n);
for i=1:n
    ai=A(:,i);
    indi=1+d*(i-1):d*i;
    for j=1:n
        if i~=j
        aj=A(:,j);
        indj=1+d*(j-1):d*j;
        aimaj=ai-aj;
        temp=K(i,j)*(aimaj*aimaj')/(aimaj'*aimaj);
        S(indi,indj)=-temp;
        S(indi,indi)=S(indi,indi)+temp;
        end
    end
end
tildn=(1:d)'*ones(size(tild))+ones(d,1)*((tild-1)*d);
tildn=tildn(:);
S=S(tildn,tildn);
[V,~,~]=svd(S);
v=V(:,end:-1:end-1);
end

