Dipole Antenna : 2

Characteristics of dipole antenna.

%%
% characteristics of dipole antenna
% author : Leon
% email:yangli0534@yahoo.com
% url : www.cnblogs.com/hiramlee0534

%%
close all;
clear all;
clc;
global k;
global l;

M = 50;% numbers pattern of different length
N = 1000; % numbers of theta
theta = linspace(-pi, pi, N);
lamda = 1 ; % wavelength
k = 2*pi/lamda; % wave number

%%
% radiation pattern
fig = figure('color','white');
title('radiation pattern of dipole antenna');

for i = 1:1:M
    l = i/M*2*lamda;
f =(cos(k*l*cos(theta))-cos(k*l))./sin(theta);
f = abs(f)/max(abs(f));
polar(theta,f);
xlabel(strcat('antenna length is  ', num2str(l*2),' lamda'));
pause(0.1);
MF(i) = getframe(gcf);
end

%%
% radiation resistance
for i = 1:1:M
    l = i/M*2*lamda;
    r(i) = quad('Fr',0, pi);
end
figure;
plot((1:1:M)/M*2*lamda,r);
title('radiation resistance of dipole antenna');
xlabel('antenna length');
ylabel('radiation resistance');

%%
%dipole antenna calc
%author:Li Yang 
% email: yangli0534@yahoo.com

function [beamwidth,D0,Rr,Rin,Xin]=dipole_pattern(len,flag );

format long

if nargin<2
   if nargin == 0
       error('too few argument')
   else
   flag = 0;
   end
end

N = 1000;
theta = linspace(0,pi,N)+pi/5/N;
phi = linspace(0,2*pi,N);

I0 = 1 ;% constant current
freq = 1e9 ;% frequency 
e0 = 8.854187817e-12;
u0 = 4*pi*1e-7;
c = 1.0/sqrt(e0*u0);
lambda = c/freq;
l = len*lambda;%length
rad = 0.001*lambda;%radius of the dipole
imp = sqrt(u0/e0);
k = 2*pi/lambda;
r = 10;
if k*r <1
    error('kr < 1');
end   
%%
%electric field intensity and magnetic field intensity in far-field region
E0 = 1i*imp*I0/(2*pi*r)*exp(-1i*k*r)*(cos(k*l/2*cos(theta))-cos(k*l/2))./sin(theta);%the theta component of e field
H0 = 1i*I0/(2*pi*r)*exp(-1i*k*r)*(cos(k*l/2*cos(theta))-cos(k*l/2))./sin(theta);% the phi component of h compoent
%plot(theta,E0)
%%
%power density
Wav = 0.5*real(E0.*conj(H0));
Wavx= (Wav.*sin(theta))'*cos(phi);
Wavy= (Wav.*sin(theta))'*sin(phi);
Wavz= repmat((Wav.*cos(theta))',1,N);
Wavc= repmat((Wav)',1,N);
if flag == 1
figure;
surf(Wavx,Wavy,Wavz,Wavc);
shading interp;
light;
lighting gouraud;
colorbar;
title('average power density');
end
%%
%radiation intensity
U =r^2*Wav;
Ux= (U.*sin(theta))'*cos(phi);
Uy= (U.*sin(theta))'*sin(phi);
Uz= repmat((U.*cos(theta))',1,N);
Uc= repmat((U)',1,N);
if flag == 1
figure;
surf(Ux,Uy,Uz,Uc);
shading interp;
light;
lighting gouraud;
colorbar;
title('radiation intensity');
end
%Prad = imp*pi/3*(I0*l/lambda)^2;
Prad = 0;
for i = 1:1:N
    Prad  = Prad +U(i)*sin(theta(i))*(theta(2)-theta(1))*2*pi;
end
D = 4*pi*U/Prad;
D0 = max(D);
Dx= (D.*sin(theta))'*cos(phi);
Dy= (D.*sin(theta))'*sin(phi);
Dz= repmat((D.*cos(theta))',1,N);
Dc= repmat((D)',1,N);
if flag == 1
figure; 
surf(Dx,Dy,Dz,Dc);
shading interp;
light;
lighting gouraud;
colorbar;
title('directivity');
end


% ND = 10*log10(D(1:round(N/2))/max(D(1:round(N/2))));
ND = 10*log10(D/max(D));
peak = find(ND==max(ND));
NDD= abs(ND+3);
left = find(NDD(1:peak)==min(NDD(1:peak)));
find((NDD)==min(NDD));
right = find(NDD(peak:end)==min(NDD(peak:end)))+ peak-1;
beamwidth = (theta(right)-theta(left))/pi*180;
if flag == 1
figure;
plot(theta/pi*180,ND);
title('E plane directivity pattern');
end
str = strcat('3dB beam width is  ',num2str(round(beamwidth*100)/100), '° ');

E_dB=10*log10(abs(E0)/max(abs(E0)));
E_dB(E_dB<=-60)=-60;
if flag == 1
figure;
polar_dB(theta*180/pi,E_dB,-60,0,4);
title('Elevation plane normalized amplitude pattern (dB)','fontsize',16);
end
D_dB=10*log10(D);
D_dB(D_dB<=-60)=-60;
if flag == 1
figure;
polar_dB(theta*180/pi,D_dB,-60,max(D_dB),4);
title(strcat('Elevation plane directivity pattern (dB):',str),'fontsize',12);
end
%text(60,0,str)

Rr = 2*Prad;
Rin = Rr/(sin(k*l/2)^2);
Xm=30*(2*si(k*l)+cos(k*l)*(2*si(k*l)-si(2*k*l))- ...
          sin(k*l)*(2*ci(k*l)-ci(2*k*l)-ci(2*k*rad^2/l)));
Xin=Xm/(sin(k*l/2))^2;
 
function [y]=si(x);
v=linspace(0,x/pi,500);
dv=v(2)-v(1);
y=pi*sum(sinc(v)*dv);

function [y]=ci(x);
v=linspace(0,x/(2*pi),500);
dv=v(2)-v(1);
y1=2*pi*sum(sinc(v).*sin(pi*v)*dv);
y=.5772+log(x)-y1;