helmholz.m

来自「ok &ouml nemli bir kaybnkao lasdhlas sdo」· M 代码 · 共 88 行

% FINITE DIFFERENCE SOLUTION OF HELMHOLZ EQUATION 
%     FOR TM MODES IN RECTANGULAR WAVEGUIDE 
%		EIGENVALUES AND PLOT OF THE AXIAL ELECTRIC FIELD : Ez 
%Note that true eigenvalue is given for the dominant TM(1,1) mode and can be used for convergence test
clear all
K = input ('Number of mesh points along i=x   ');
L = input ('Number of mesh points along j=y   ');
a = input ('Length along x = a   =  ');
h = a/(K+1);
b = h * (L+1)
% Set coefficient matrix D 
Order=K*L;
%initialize
for i=1:Order
    for j=1:Order
    D(i,j)=0;
    end
end
k=1;
for i=1:Order
   D(i,i)=4;
   if i>k*L
      k=k+1;
   end
   kL=k*L;
      if i > 1+(k-1)*L
         D(i,i-1)=-1;
      end
      if i < kL
         D(i,i+1)=-1;
      end
      if i > L
         D(i,i-L)=-1;
      end
      if i <= Order - L 
         D(i,i+L)=-1;
      end
   end
%   Diagonal elements of c are eigenvalues of D
%   to ith eigenvalue corresponds the ith column vector of V   
[V,C]=eig(D);
% gamman is the column vector of eigenvalues.
gamman=sqrt(diag(C))/h;
% E is gamman obtained by organizing the elements of gamman in ascending order
% o(i) is an index which denotes the position of i-th element of E in gamman 
[E,o]=sort(gamman);
% Calculate correct eigenvalues of the modes 
kont1=0;
for i=1:2*K
   for j=1:2*L
      kont1=kont1+1;
      kapareal(kont1)=sqrt((i/a)^2+(j/b)^2)*pi;
   end
end
R=sort(kapareal);
for i=1:Order
   R1(i)=R(i);
end
kont = input ('enter 1 if plot is desired ' );
while kont == 1 
   modeno=input (' enter mode number   ')
   colno=o(modeno);
   result=E(modeno)
   correct_value=R1(modeno)
   percentage_error=(correct_value-result)/correct_value*100
   B=V(:,colno);
% plot Ez 
   for i=1:K
   for j=1:L
      x(i,j)=i;
      y(i,j)=j;
      z(i,j)=B(j+(i-1)*L);
   end
end
% select plot type
p1=input ('enter 1 for surface_plot and any number for contour_plot');
if p1==1
surf(x,y,z)
else 
   level = input ('enter number of desired equal amplitude contours');
   contour(x,y,z,level)
end
kont=input('enter 1 if another plot is desired');
end