clear all;format long;lambda  = 1;delka=0.5;polomer=0.01;l=delka;a=polomer;l       = l*lambda;          % length of the dipolea       = a*lambda;    		  % radius of the antenna wiremi      = 4*pi*1e-7;         % permeability of vacuumepsilon = 8.85e-12;          % permittivity of vacuumN       = 33;                % number of segmentsdelta   = l/(N+1);  		     % length of one segmentk	     = 2*pi/lambda;       % wavenumberomega   = k*3e+8;            % angular frequencyalpha   = 0.5*delta;         % one half of the segment's lengthpsi = zeros( 1, N+1);        % numerical integrationfor m=1:(N+1)  x = (m-1)*delta;  psi(m) = quad( 'green', -alpha, alpha, 1e-5, [], x, a, k);endpsi = psi/delta;Z = zeros( N, N);            % computing impedance matrixmulti = 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);  endendY = inv( Z);               % computing admitance matrixI = Y( :, 17);             % current distribution on the antennaclear psi Zmodule = 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 phasefigure(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)));endpolar( theta, E/max(E));out = 1/Y( 17, 17);     % input impedance of the antenna