%rectangular waveguide cutoff frequency eigenvalue solver

%enter dimensions (only built to enter even numbered values)
a=4;    %width
b=2;    %height

%enter mesh settings
h=1;

jmax = a/h+1; %x axis
imax = b/h+1; %y axis

%total number of nodes
K = imax*jmax;

%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

y=eig(A)

%equations
u = 4*pi*10^-7;
e = 8.85*10^-13;
Y = 1;
N = 0;

more = 1;
while more ~= N    
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