%% solve d^2/dt^2 u = c^2*d^2/dx^2 u
% use modified forward Euler

%% inputs
Nx=22;           % number of x-nodes
Lx=1;            % size of domain
dt=.002;         % time step
TT=10.24;        % total simulation time
c=1;             % speed of wave

%% ploting
plotskip=10;  % number of steps to skip before plot

%% calc
dx=Lx/(Nx-1);           % number of x intervals
x=(0:dx:Lx)';           % x-node positions
k=(1:Nx-2);             % eigenmodes numbers
Nt=fix(TT/dt);          % number of time steps
t=(0:Nt-1)*dt;          % time positions
f=(-Nt/2:Nt/2-1)/Nt/dt; % frequencies for fft

%% save u,v
us=zeros(Nx-2,Nt);
vs=zeros(Nx-2,Nt);

%% spatial derivatives: K=-(d/dx)^2
K=ktbc(Nx-2);    % Ku=u''
[Q,L]=eig(K);    % eigensystem for K w^2=L*(c/dx)^2

%% initial conditions
%u=Q*ones(Nx-2,1)/10; % all mode equal amplitude
u=Q(:,6);
v=zeros(Nx-2,1);     % zero 'velocity'

us(:,1)=u;
vs(:,1)=v;

% used for plotting 
umax=max(abs(u));

%% main loop
for nt=2:Nt; 
  % modified Euler integration
  v=v+(3*D*u/dx-c^2*K*u/dx^2)*dt;  % 'a=f/m'
  u=u+v*dt;             % v=du/dt
  
  % save for later
  us(:,nt)=u;
  vs(:,nt)=v;
  
  % plot every plotskip
  if (mod(nt,plotskip)==0)
    plot(x,[0;u;0]);
    axis([0 Lx 1.1*umax*[-1 1]]);
    drawnow;
  end
end

%% plots

% space-time plot
imagesc([0 TT],[1 Nx-2]*Lx/Nx,us); 
xlabel('time'); 
ylabel('position');
pause;

% mode-time plot
imagesc([0 TT],[1 Nx-2],Q'*us); 
xlabel('time'); 
ylabel('mode');
pause;

% mode-frequency (dispersion relation)
w=sqrt(diag(L*(c/dx)^2)); % eigenfrequencies
fmax=max(w)/2/pi;         % max freq w=2*pi*f

fr=[-Nt/2 Nt/2-1]/Nt/dt;  % frequency range
imagesc(fr,[1 Nx-2],abs(fftshift(fft(Q'*us,[],2),2))); 
xlabel('frequency'); 
ylabel('mode'); 
axis('xy'); 
axis([1.1*fmax*[-1 1] -inf inf]); 

hold on; 
plot(w/2/pi,k,'yo'); % dispersion relation k(w)
hold off;
pause;

% frequency-mode (dispersion relation)
w=sqrt(diag(L*(c/dx)^2)); % eigenfrequencies
fmax=max(w)/2/pi;         % max freq w=2*pi*f

fr=[-Nt/2 Nt/2-1]/Nt/dt;  % frequency range
imagesc([1 Nx-2],fr,abs(fftshift(fft(Q'*us,[],2),2))'); 
ylabel('frequency'); 
xlabel('mode'); 
axis('xy'); 
axis([-inf inf 1.1*fmax*[-1 1]]); 

hold on; 
plot(k,w/2/pi,'yo'); % dispersion relation w(k)
hold off;


