close all; clear; clc;

%% Simulate electric wave in a thin dielectric
fs = 1:0.01:5; % Driver frequencies
Us = cell(1,numel(fs));

N = 100+1;
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 = L./(h^2)/4;
L = kron(L,speye(N)) + kron(speye(N),L);
lap = @(e) L*e;

v = 1;
t = linspace(0,20,2e2); % Time axis
for fi = 1:numel(fs)
    E = zeros(N,N); e = E(:); e = [e; lap(e)];
    dedt = @(t,e) wave_2d_odefun_4(t,e,fs(fi),v,X,Y,N,lap);
    
    tic
    [~,U] = ode113(dedt,t,e);
    toc
    
    U2 = zeros(N,N,numel(t));
    for ti = 1:numel(t)
        U2(:,:,ti) = reshape(U(ti,1:N^2),N,N);
    end
    Us{fi} = U2;
    
    %{{
    figure(1); clf; colormap(plasma);
    for ti = 1:5:numel(t)
        u = U2(:,:,ti);
        imagesc(X(1,:),Y(:,1),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);
        axis equal tight;
        ylim([-1 1]*0.65);
        drawnow;
    end
    %}
    fprintf('f=%.2f\n',fs(fi));
end
fprintf('Saving...\n');
tic
save('wave_2d_Us.mat','Us','fs','t','X','Y','N','-v7.3');
toc
return

%% Load the saved simulation
close all; clear; clc;
load('wave_2d_Us.mat');

%% Animate surface
x = X(1,:); y = Y(:,1);
figure(1); clf; colormap(plasma);
Ui = 27; % 36
u = Us{Ui}(:,:,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);
xlabel('x'); ylabel('y');
title(sprintf('Driver frequency = %.2f',fs(Ui)));
axis equal tight;
ylim([-1 1]*0.65);
drawnow; pause(0.01);

for ti = 2:numel(t)
    u = Us{Ui}(:,:,ti);
    ih.CData = u;
    drawnow; pause(0.01);
end

%% Animate slice through the middle
x = X(1,:);
I = c2(x,-0.75);
Ui = 32; % 27; % 36;
u = Us{32}(round(N/2),:,1);
figure(1); clf;
ph(1) = plot(x,u); hold on;
ph(2) = plot(x(I),u(I),'o'); % Driver
plot([1 1]*0.2,[-2 2],'k--'); plot([1 1]*-0.2,[-2 2],'k--');
ylim([-1 1]*2);
xlabel('x'); ylabel('Electric field');
grid on;
drawnow; pause(0.05);

for ti = 2:numel(t)
    u = Us{32}(round(N/2),:,ti);
    ph(1).YData = u;
    ph(2).YData = u(I);
    drawnow; pause(0.05);
end

%% Average intensity in 3 regions
input = NaN(1,numel(fs));
middle = NaN(1,numel(fs));
output = NaN(1,numel(fs));
for fi = 1:numel(fs)
    x = X(1,:); y = Y(:,1);
    U = double(Us{fi});
    
    region = y>-0.5&y<0.5; % Y region
    u = squeeze(mean(U(region,:,:),1));
    
    % X regions
    rms1 = rms(u(x>-0.65&x<-0.2,:))*sqrt(2)*2;
    rms2 = rms(u(x>-0.2&x<+0.2,:))*sqrt(2)*2;
    rms3 = rms(u(x>+0.2&x<+0.65,:))*sqrt(2)*2;
    
    region = t>15&t<19.5; % Time region
    input(fi) = mean(rms1(region));
    middle(fi) = mean(rms2(region));
    output(fi) = mean(rms3(region));
    
    if ~mod(fi,round(numel(fs)/20))
        fprintf('%.1f%% done\n',fi/numel(fs)*100)
    end
end
%% 
figure(2); clf;
yyaxis left;
plot(fs,input);
ylabel('Input intensity');
yyaxis right;
plot(fs,(output)./(middle));
ylabel('Output/Middle intensity');
ylim([0 2]);
grid on; box on; axis square;
xlabel('Driving frequency');
