close all; clear; clc;

%% Setup grid
N = 60;
x = linspace(-1,1,N); h = x(2)-x(1);
[X,Y] = meshgrid(x);

% Laplacian operator
L = spdiags(bsxfun(@times,[1 -2 1],ones(N^2,1)),[-1 0 1],N,N);
L(1,end) = 0; L(end,1) = 0;
L = L./(h^2)/4;
L = kron(L,speye(N)) + kron(speye(N),L);

lap = @(e) L*e;

%% Wave in a channel with a loss
% Initial electric field
E = exp(-((X-0.5).^2+(Y).^2)*50);
%E = exp(-((X-0.5).^2)*50);
%E = peaks(N)/40;

% Set up differential eqn.
vph = 2e-2;
del2e = @(e) vph^2.*lap(e);
e = E(:); e = [e; del2e(e)];
dedt = @(t,e) wave_2d_odefun_1(t,e,X,Y,N,del2e);

tic
[t,U] = ode113(dedt,linspace(0,1e3,1e3),e);
toc

%% Animate
figure(1); colormap(plasma);
u = reshape(U(1,1:N^2),N,N);
sh = surf(u);
zlim([-1 1]*0.5); caxis([-1 1]*0.5);
grid on; box on; axis square;
drawnow;

for ti = 2:numel(t)
    sh.ZData = reshape(U(ti,1:N^2),N,N);
    drawnow;
end

%% With a driver
[X,Y] = meshgrid(x);
E = zeros(N,N);

v = 1;
e = E(:); e = [e; lap(e)];
f = 2;
dedt = @(t,e) wave_2d_odefun_2(t,e,f,v,X,Y,N,lap);

tic
[t,U] = ode113(dedt,linspace(0,10,1e2),e);
toc

%% Animate
figure(1); colormap(plasma);
u = reshape(U(1,1:N^2),N,N);
sh = surf(u);
zlim([-1 1]*2); caxis([-1 1]*2);
grid on; box on; axis square;
drawnow;

for ti = 2:numel(t)
    sh.ZData = reshape(U(ti,1:N^2),N,N);
    drawnow; pause(0.05);
end

%% With an etalon
v = 1; % Phase velocity
E = zeros(N,N); e0 = E(:); e0 = [e0; lap(e0)]; % Starting point
f = 3.00; % Driver frequency
dedt = @(t,e) wave_2d_odefun_3(t,e,f,v,X,Y,N,lap);

tic
[t,U] = ode113(dedt,linspace(0,3,2e2),e0);
toc

% Reshape back to 2d coordinates
U_square = zeros(N,N,numel(t));
for ti = 1:numel(t)
    U_square(:,:,ti) = reshape(U(ti,1:N^2),N,N);
end

%% Animate
x = X(1,:); y = Y(:,1);
figure(1); clf; colormap(plasma);
u = U_square(:,:,1);
ih = imagesc(x,y,u); hold on;
plot([1 1]*-0.2,[-1 1],'k--'); plot([1 1]*0.2,[-1 1],'k--');
plot([-1 1],[1 1]*0.5,'k-'); plot([-1 1],[1 1]*-0.5,'k-');
caxis([-1 1]*2);
grid on; axis equal tight;
ylim([-0.65 0.65]); xlim([-0.65 0.65]);
xlabel('x'); ylabel('y');
drawnow; pause(0.01);

for ti = 2:numel(t)
    ih.CData = U_square(:,:,ti);
    drawnow;
end

