current.m

来自「method of moment to solve Current distri」· M 代码 · 共 72 行

clear all;
format long;
lambda  = 1;

delka=0.5;
polomer=0.01;

l=delka;
a=polomer;
l       = l*lambda;          % length of the dipole
a       = a*lambda;    		  % radius of the antenna wire

mi      = 4*pi*1e-7;         % permeability of vacuum
epsilon = 8.85e-12;          % permittivity of vacuum

N       = 33;                % number of segments
delta   = l/(N+1);  		     % length of one segment

k	     = 2*pi/lambda;       % wavenumber
omega   = k*3e+8;            % angular frequency
alpha   = 0.5*delta;         % one half of the segment's length

psi = zeros( 1, N+1);        % numerical integration
for m=1:(N+1)
  x = (m-1)*delta;
  psi(m) = quad( 'green', -alpha, alpha, 1e-5, [], x, a, k);
end
psi = psi/delta;

Z = zeros( N, N);            % computing impedance matrix
multi = j*omega*mi;
divid = j*omega*epsilon;
for m=1:N
    for n=m:N
        dist  = abs(m-n);
        hlp   = 2*psi(1+dist) - psi(1+abs(dist-1)) - psi(1+abs(dist+1));
        Z(m,n)= multi*(delta^2)*psi(1+dist) + hlp/divid;
        Z(n,m)=Z(m,n);
  end
end

Y = inv( Z);               % computing admitance matrix
I = Y( :, 17);             % current distribution on the antenna
clear psi Z

module = abs( I) * 1e+3;
phase  = angle( I);
                           % plotting current distribution
% plot( 1:N, module, 'b-', 1:N, phase, 'g.');
figure(1)
plot( 1:N, module);
%clear module phase

figure(2)
plot(1:N,phase)
                           % directivity pattern
                           
 figure(3);                          
theta = zeros( 1, 361);
E     = zeros( 1, 361);
r     = zeros( 1, N);
for m=1:361
  theta(m) = m*pi/180;
  for n=1:N
    r(n) = (n-1)*delta*cos( theta(m));
  end
  E(m) = abs( conj(I')*exp( -j*k*r')*sin( theta(m)));
end

polar( theta, E/max(E));

out = 1/Y( 17, 17);     % input impedance of the antenna