% Mika Juntunen / 12.2.2007
%
% This script computes the dynamics of simple rod
% 

% Define coefficienst of the rod
E = 10;
A = 1;
density = 1;
L = 1;
c = sqrt(density/E);

% Define intial conditions
f = @(x)(zeros(size(x)));   % deformation
g = @(x)(10*(x.*(x-L)));  % velocity

% Compute coefficients
N = 10;
f_fun = @(x,k)(2/L*sin((2*k-1)*pi*x/2/L));
g_fun = @(x,k)(4*c/(2*k-1)/pi*sin((2*k-1)*pi*x/2/L));

Ck = zeros(N,1);
Dk = zeros(N,1);

for i = 1:N
  apu_f = @(x)(f(x).*f_fun(x,i));
  apu_g = @(x)(g(x).*g_fun(x,i));
  
  Ck(i) = quad(apu_g,0,L);
  Dk(i) = quad(apu_f,0,L);
end

% plot solution
x = (0:0.05:1)';
fx = @(x,k)(sin((2*k-1)*pi*x/L/2));
ft_c = @(t,k)(sin((2*k-1)*pi*t/c/L/2));
ft_d = @(t,k)(cos((2*k-1)*pi*t/c/L/2));

i = 0;
for t=0:0.01:4
  y = ratkaisu(x,t,Ck,Dk,fx,ft_c,ft_d);
  figure(1)
  plot(x+y,zeros(size(x)),'-o')
  axis([-0.1,2*L+0.1,-1,1])
  drawnow
  pause(0.1)
  
%  M(i+1) = getframe(gcf);
end

%movie(M)
