% This LaTeX was auto-generated from MATLAB code.
% To make changes, update the MATLAB code and export to LaTeX again.

\documentclass{article}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{lmodern}
\usepackage{graphicx}
\usepackage{color}
\usepackage{listings}
\usepackage{hyperref}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{epstopdf}
\usepackage[table]{xcolor}
\usepackage{matlab}

\sloppy
\epstopdfsetup{outdir=./}
\graphicspath{ {./ssy200_ass2_images/} }

\begin{document}

\matlabheading{SSY200 Assignment 2}


\vspace{1em}
\begin{par}
\begin{flushleft}
Set up environment
\end{flushleft}
\end{par}

\begin{matlabcode}
format compact
EM_constants
close all
\end{matlabcode}

\begin{par}
\begin{flushleft}
Parameters, as given in §2.4
\end{flushleft}
\end{par}


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

\matlabheadingtwo{Generate grid}

\begin{enumerate}
\setlength{\itemsep}{-1ex}
   \item{\begin{flushleft} Glass boundaries (-a,+a) fall between grid points. \end{flushleft}}
   \item{\begin{flushleft} Glass boundaries fall exactly on grid points. \end{flushleft}}
\end{enumerate}


\begin{matlabcode}
grid_type = "G1";
N = 10; % Assert N >= 4 and even.
assert(N >= 4, "Assert failed: N must be >4")
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;
        h = 2*a/N;
        xn = (n + 1/2)*h; % Grid-points.
    
        N_gp = length(xn); % Number of grid points.
    case "G2"
        %% G2: Glass boundaries on grid points!
        n = -N:N;
        h = 2*a/N;
        xn = n*h;
    
        N_gp = length(xn); % Number of grid points.
end
\end{matlabcode}

\matlabheadingtwo{Create system of linear equations for the 1D Helmholtz equation}

\begin{par}
\begin{flushleft}
Create square tri-diagonal A matrix.
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
Each \textbf{row }corresponds to a position along the x-axis.
\end{flushleft}
\end{par}

\begin{par}
\begin{flushleft}
$A=\left\lbrack \begin{array}{cccc}
-\frac{1}{h}-\frac{{\mathrm{jk}}_0 }{2} & \frac{1}{h}-\frac{{\mathrm{jk}}_0 }{2} & 0 & \cdots \\
-\frac{1}{h^2 } & \frac{2}{h^2 }+\mu_0 {\left\lbrack \cdots \right\rbrack }_1  & -\frac{1}{h^2 } & 0\\
0 & -\frac{1}{h^2 } & \frac{2}{h^2 }+\mu_0 {\left\lbrack \cdots \right\rbrack }_2  & -\frac{1}{h^2 }\\
\vdots  & 0 & -\frac{1}{h^2 } & \frac{2}{h^2 }+\mu_0 {\left\lbrack \cdots \right\rbrack }_3 
\end{array}\right\rbrack$, ie. row 1 is the left boundary.
\end{flushleft}
\end{par}

\begin{matlabcode}
A = zeros(N_gp, N_gp); % Initialize
% Helmholtz (non-boundary cases)
% Initialize main diagonal.
for i = 1:N_gp-1 
    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 
\end{matlabcode}

\begin{par}
\begin{flushleft}
Insert boundary conditions in A matrix.
\end{flushleft}
\end{par}

\begin{matlabcode}
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.
\end{matlabcode}

\begin{par}
\begin{flushleft}
Construct b column vector (right hand side at each point in space).
\end{flushleft}
\end{par}

\begin{par}
$$b=\left\lbrack \begin{array}{c}
-2{\mathrm{jk}}_0 E_z^i \frac{\left(x_1 +x_0 \right)}{2}\\
0\\
\vdots \\
0
\end{array}\right\rbrack$$
\end{par}

\begin{matlabcode}
b = [-2*j*k0*(E0i*exp(-j*k0*(xn(2)+xn(1))/2)); zeros(N_gp-1,1)];
\end{matlabcode}

\begin{par}
\begin{flushleft}
Solve
\end{flushleft}
\end{par}

\begin{matlabcode}
z = A\b; % Solve system.
\end{matlabcode}

\matlabheadingtwo{Plot total electric field ${\mathit{\mathbf{E}}}_z \left(x\right)$}

\begin{matlabcode}
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('G1')
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
\end{matlabcode}
\begin{center}
\includegraphics[width=\maxwidth{56.196688409433015em}]{figure_0}
\end{center}


\vspace{1em}
\matlabheadingtwo{Analytical values for $k_n$, R, T}

\begin{matlabcode}
k1 = sqrt(er*k0^2 - j*omega*mu0*sigm)
\end{matlabcode}
\begin{matlaboutput}
k1 = 15.9969 - 2.3567i
\end{matlaboutput}
\begin{matlabcode}
Delta =  (k0 + k1)^2*exp(j*4*a*k1) - (k0 - k1)^2; % Intermediate value (not used elsewhere)
R = (k0^2 - k1^2)/Delta * exp(j*2*a*k0) * (exp(j*4*a*k1)-1)
\end{matlabcode}
\begin{matlaboutput}
R = -0.1716 - 0.2049i
\end{matlaboutput}
\begin{matlabcode}
T = k0*k1/Delta * 4*exp(j*2*a*(k0+k1))
\end{matlabcode}
\begin{matlaboutput}
T = 0.8242 - 0.2129i
\end{matlaboutput}

\end{document}
