%% Propagation constant k_z(omega) for TE10 mode. 
% Shows cutoff frequency and regions where mode is propagating or
% evanescent.
close all, clear all
c0=299792458;
f = 0:1e7:1.2e10;
w = 2 * pi * f; 
a = 0.04; %Lx
b = 22.5e-3; % Ly
Lz = 160e-3;

kz = sqrt((w/c0).^2 - (pi/a)^2)

S21 = exp(-j*kz*Lz)

plot(f, abs(S21))
prettify_axes(gca, 3.9e9, 6.5e9)
xlim([3.9e9, 6.5e9])
ylim([0, 2])

%%

% Propagation constants for modes
gamma_mn = @(m,n,w) sqrt((m*pi/a)^2 + (n*pi/b)^2 - (w/c0).^2)

fc_mn = @(m,n) sqrt((m*pi/a)^2 + (n*pi/b)^2)*c0/(2*pi)

% Plots
figure; ax = gca;
hold on
plot(f, imag(gamma_mn(1,0,w)))
plot(f, imag(gamma_mn(0,1,w)))
plot(f, imag(gamma_mn(1,1,w)))
plot(f, imag(gamma_mn(2,1,w)))

plot(f, w/c0, "k--")
fc10 = fc_mn(1,0)
fc01 = fc_mn(0,1)

titletext = sprintf('$TE_{10}$, $f_c$ = %.4e Hz', fc10);
title(titletext, 'interpreter','latex')
xlabel('$f$','interpreter','latex')
ylabel('$k_z(\omega)$ ','Interpreter','latex')


x1 = 3.9e9; x2  = 6.5e9; % shade x region
prettify_axes(ax, x1, x2)


legend(ax,'','$TE_{10}$', '$TE_{01}$', '$TE_{11}$', '$TE_{21}$', 'interpreter', 'latex', 'FontSize', 14, 'Location', 'best');


%% Modes in narrow section 
a = 30e-3;
b = 22.5e-3; % Ly
fc_mn = @(m,n) sqrt((m*pi/a)^2 + (n*pi/b)^2)*c0/(2*pi)

fc_mn(1,0)
fc_mn(1,1)
fc_mn(0,1)