% unisson bi ou trichorde----------------
clear all; format short e;
T = 5;

% la table-------------------------------
fc = 50;
mc = 300;
kc = (2*pi*fc)^2*mc;
critik = 2*sqrt(kc*mc);
c = 5*critik;

% les cordes (masse m1 = 1) -------------
N = 2;
%f = [500 500.00];
f = [500 500.1];
%N = 3;
%f = [500 500.1 500.07];
k = (2*pi*f).^2;
ks = sum(k);

% assemblage-----------------------------
A = zeros(2*N+2,2*N+2);
for i = 1:N+1
    A(i,N+1+i) = 1;
end

for j = 1:N
    A(N+1+j,j) = -k(j);
end
A(N+2:2*N+1,N+1) = k'/2;
%masse table a la fin
A(2*N+2,1:N) = k/2;
A(2*N+2,N+1) = -(kc+ks/2);
A(2*N+2,2*N+2) = -c;
A(2*N+2,:) = A(2*N+2,:)/mc;

% resolution-----------------------------
tau = 1/f(1)/10;
M = expm(tau*A);
FS = 1/tau;
t = 0:tau:T;
nt = length(t);
son = zeros(1,nt);

nt = length(t);
X = [zeros(1,N+1) 1 .9 0]';
son(1) = X(N+1);

for n = 2:nt
    X = M*X;
    son(n) = X(N+1);
end

% resultats-----------------------------
son = .99*son/max(abs(son));
plot(t,20*log10(abs(son)+1e-3));
sound(son,FS);
