%notched rectangular waveguide cutoff frequency eigenvalue solver

%enter dimensions (only built to enter even numbered values)
a=input('enter width of waveguide ');    %width
b=input('enter height of waveguide ');    %height

an=input('enter width of notch ');
if an >= a
    an = input('error! enter notch width smaller than waveguide width ');
end

bn=input('enter height of notch');
if bn >= b
    bn = input('error! enter notch height smaller than waveguide height ');
end

%enter mesh settings
h=input('enter mesh size ');

jmax = a/h+1; %x axis
imax = b/h+1; %y axis
jn = an/h+1;
in = bn/h+1;

%total number of nodes
K = imax*jmax;
Kn = in*jn;

Ko = K-Kn+2*in+jn-2; 
%calculate boundary location
nbc = 2*imax+2*jmax-4;
[X Y] = meshgrid(1:jmax,1:imax);
mdl = (Y-1)*jmax + X;
mdl(2:end-1,2:end-1) = 0;

%top boundary
for n = 1:jmax
    top(n,1) = n;
end

%bottom
for n = 1:jmax
    bot(n,1) = (imax-1)*jmax+n;
end

%left
for n = 1:imax
    left(n,1) = jmax*n-(jmax-1);
end

%right
for n = 1:imax
    right(n,1) = jmax*n;
end

%Build A Matrix
A = eye(K)*4;
for k = 1:K
    if k-1 > 0
        A(k,k-1) = 1;
    end
    if k+1 < K
        A(k,k+1) =1;
    end
    if k-jmax > 0
        A(k,k-jmax) = 1;
    end
    if k+jmax < K
        A(k,k+jmax) = 1;
    end
end

%Enforce Boundary Conditions
%top
for n = 1:jmax
    A(n,n+jmax) = 2;
end

%bottom
for n = 1:jmax
    p = (imax-1)*jmax+n; %index of x value
    A(p,p-jmax) = 2;
end

%left
for n = 1:imax
    q = (n*jmax-(jmax-1));
    if q-1 > 0
        A(q,q-1) = 0;
    end
    if q+1 < K
        A(q,q+1) = 2;
    end
end

%right
for n = 1:imax
    q = n*jmax;
    if q+1 < K
        A(q,q+1) = 0;
    end
    if q-1 > 0
        A(q,q-1) = 2;
    end
end

yc=eig(A)

%equations
u = 4*pi*10^-7;
e = 8.85*10^-13;
Y = 1;
N = 0;
more = 1;
while more ~= N
yc
m = input('enter x value of node 1-jmax ' );
n = input('enter y value of node 1-imax ');

fc = 1/2/pi/sqrt(u*e)*sqrt((m*pi/a)^2+(n*pi/b)^2)
more = input('Calculate another node? Type Y or N (CAPS) ');
end