%% Annnnnnnna CEM Inlämning 2
format compact
EM_constants
close all

%% Assignment parameters
omega = 3e9;   % Given in §
k0 = omega/c0; % Wavenumber for given frequency (in vacuum).
sigm = 0.02;   % [S/m]
er = 2.5;           % Relative permittivity for glass.
eps_glass = e0*er;  % Abs. permittivity for glass.
E0i = 1;            % Magnitude of incident wave. Arbitrary?

grid_type = "G1"
% N >= 4 and even.
N = 10;
assert(mod(N,2) == 0, "Assert failed: N must be even!")
a = 2e-2; % Half-thickness of glass. 2cm

for grid_type = ["G1" "G2"]
switch grid_type
    case "G1"
    %% G1: Glass boundaries between grid-points.   
    n = -N-1:N
    dx = 2*a/N
    xn = (n + 1/2)*dx % Grid-points.

    N_gp = length(xn); % Number of grid points.
    case "G2"
    %% G2: Glass boundaries on grid points!
    n = -N:N
    dx = 2*a/N
    xn = n*dx

    N_gp = length(xn); % Number of grid points.
    end
%% Finite differences
h = dx; % Step length.

A = zeros(N_gp, N_gp);

% Helmholtz (ie. non-boundary cases)
for i = 1:N_gp-1 % Initialize main diagonal.
    if (abs( xn(i) ) <= a) % Inside glass!!
        d = 2/h^2 + mu0*(j*omega*sigm - omega^2*eps_glass);
    else % Else in air.
        d = 2/h^2 + mu0*(-omega^2*e0);
    end
    A(i,i) = d;
end
s = -1/h^2;% Upper and lower diagonal entries.
for i = 1:N_gp-1
    A(i,i+1) = s; % Upper diagonal.
    A(i+1,i) = s; % Lower diagonal.
end 

% lambda = eig(A)% dunno lol delet this

A(1,1:2) = [-1/h-j*k0/2, 1/h-j*k0/2];                % Left boundary condition.
A(N_gp, end-1:end) = [-1/h+j*k0/2, 1/h+j*k0/2];      % Right boundary condition.

b = [-2*j*k0*(E0i*exp(-j*k0*(xn(2)+xn(1))/2)), zeros(1,N_gp-1)]';

z = A\b; % Solve system.



subplot(211)
plot(xn, angle(z))
grid on
ylabel('Phase')
xlabel('x')
hold on
plot([-a -a],[min(angle(z)) max(angle(z))],'r--','LineWidth',2.0) % Draw (-a,+a) on plot.
plot([+a +a],[min(angle(z)) max(angle(z))],'r--','LineWidth',2.0) % Draw (-a,+a) on plot.

subplot(212) 
plot(xn, abs(z))
grid on
hold on
legend(grid_type)
ylabel('Magnitude')
xlabel('x')
plot([-a -a],[min(abs(z)) max(abs(z))],'r--','LineWidth',2.0) % Draw (-a,+a) on plot.
plot([+a +a],[min(abs(z)) max(abs(z))],'r--','LineWidth',2.0) % Draw (-a,+a) on plot.


end

legend(["G1" "G2"])

%% 2.4
% Analytical values for reflection and transmission
k1 = sqrt(er*k0^2 - j*omega*mu0*sigm)
Delta =  (k0 + k1)^2*exp(j*4*a*k1) - (k0 - k1)^2

R = (k0^2 - k1^2)/Delta * exp(j*2*a*k0) * (exp(j*4*a*k1)-1)
T = k0*k1/Delta * 4*exp(j*2*a*(k0+k1))



% Empirical values 
Tnum = z'./E0i./exp(-j*k0.*xn)

% plot(abs(Tnum))
% hold on
% plot(angle(Tnum))

abs(T)
abs(Tnum(end))

% Compute magnitude of reflected wave E0r (at -a)
% for i = 1:100
%     E0r = (z(i) - E0i*exp(-j*k0*xn(i))) / exp(j*k0*xn(i));
%     abs(E0r)
% end