%--------------------------------------------------------------------------
%
% Kuramoto Sivashinsky equation
%
% Response to a spatially localised initial condition.
% 
% The time response of the plant to an Gaussian shaped initial condition is
% calculated: the generated wavepacket travels downstream while growing and
% is detected by the two outputs y and z.
% The spatio-temporal time evolution of v(x,t) is plotted, as the the two
% output signals.
% 
% Reference figure(s): Fig. 4 
%
% 
% Nicolo' Fabbiane, March 2014
% KTH Mechanics
% nicolo@mech.kth.se
%
%--------------------------------------------------------------------------
clc, clear all
fprintf('\nResponse to a spatially localised initial condition\n\n')




%% Initialization of the plant --------------------------------------------

% LTI system matrix A and x grid
nq = 400; % number of grid points

[A,x,I] = KS_init(nq);


% Inputs matrix B

% disturbance d (Gaussian shape at x_d with sigma_d variance)
x_d = 35; sigma_d = 4; 
Bd  = exp(-((x-x_d)/sigma_d).^2)/sigma_d;

% actuator u (Gaussian shape at x_u with sigma_u variance)
x_u = 400; sigma_u = 4; 
Bu  = exp(-((x-x_u)/sigma_u).^2)/sigma_u;


% Outputs matrix C

% measurement y (Gaussian shape at x_y with sigma_y variance)
x_y = 300; sigma_y = 4; 
Cy  = (exp(-((x-x_y)/sigma_y).^2)/sigma_y).' * I;

% output z (Gaussian shape at x_z with sigma_z variance)
x_z = 700; sigma_z = 4; 
Cz  = (exp(-((x-x_z)/sigma_z).^2)/sigma_z).' * I;

        

        
%% Time integration -------------------------------------------------------

% Parameters
dt   = 1;
Tend = 2000;

t  = 0:dt:Tend;
nt = length(t);


% Time-stepper matrices (Crank-Nicholson)
CNi = (speye(nq) - dt/2*A);
CNe = (speye(nq) + dt/2*A);


% Variables initialization
v = zeros(nq,nt); % velocity
y = zeros(1 ,nt); % measurement signal
z = zeros(1 ,nt); % output signal


% Initial condition (Gaussian shape at x_d with sigma_d variance)
x_d = 35; sigma_d = 4; 
v(:,1) = exp(-((x-x_d)/sigma_d).^2)/sigma_d;


% Time-loop
for i = 1:nt-1
    
    % outputs
    y(:,i) = Cy*v(:,i);
    z(:,i) = Cz*v(:,i);
    
    % timestep
    v(:,i+1) = CNi\(CNe*v(:,i));
    
    % runtime output
    runtime_output

end




%% Figure 4 ---------------------------------------------------------------
figure(4); clf

imagesc(x,t,v')
axis xy; axis([0 800 0 Tend]);
caxis([-1 1]*2.5)
xlabel('x'); ylabel('t')