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');

  1 %%
  2 %dipole antenna calc
  3 %author:Li Yang 
  4 % email: yangli0534@yahoo.com
  5 
  6 function [beamwidth,D0,Rr,Rin,Xin]=dipole_pattern(len,flag );
  7 
  8 format long
  9 
 10 if nargin<2
 11    if nargin == 0
 12        error('too few argument')
 13    else
 14    flag = 0;
 15    end
 16 end
 17 
 18 N = 1000;
 19 theta = linspace(0,pi,N)+pi/5/N;
 20 phi = linspace(0,2*pi,N);
 21 
 22 I0 = 1 ;% constant current
 23 freq = 1e9 ;% frequency 
 24 e0 = 8.854187817e-12;
 25 u0 = 4*pi*1e-7;
 26 c = 1.0/sqrt(e0*u0);
 27 lambda = c/freq;
 28 l = len*lambda;%length
 29 rad = 0.001*lambda;%radius of the dipole
 30 imp = sqrt(u0/e0);
 31 k = 2*pi/lambda;
 32 r = 10;
 33 if k*r <1
 34     error('kr < 1');
 35 end   
 36 %%
 37 %electric field intensity and magnetic field intensity in far-field region
 38 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
 39 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
 40 %plot(theta,E0)
 41 %%
 42 %power density
 43 Wav = 0.5*real(E0.*conj(H0));
 44 Wavx= (Wav.*sin(theta))'*cos(phi);
 45 Wavy= (Wav.*sin(theta))'*sin(phi);
 46 Wavz= repmat((Wav.*cos(theta))',1,N);
 47 Wavc= repmat((Wav)',1,N);
 48 if flag == 1
 49 figure;
 50 surf(Wavx,Wavy,Wavz,Wavc);
 51 shading interp;
 52 light;
 53 lighting gouraud;
 54 colorbar;
 55 title('average power density');
 56 end
 57 %%
 58 %radiation intensity
 59 U =r^2*Wav;
 60 Ux= (U.*sin(theta))'*cos(phi);
 61 Uy= (U.*sin(theta))'*sin(phi);
 62 Uz= repmat((U.*cos(theta))',1,N);
 63 Uc= repmat((U)',1,N);
 64 if flag == 1
 65 figure;
 66 surf(Ux,Uy,Uz,Uc);
 67 shading interp;
 68 light;
 69 lighting gouraud;
 70 colorbar;
 71 title('radiation intensity');
 72 end
 73 %Prad = imp*pi/3*(I0*l/lambda)^2;
 74 Prad = 0;
 75 for i = 1:1:N
 76     Prad  = Prad +U(i)*sin(theta(i))*(theta(2)-theta(1))*2*pi;
 77 end
 78 D = 4*pi*U/Prad;
 79 D0 = max(D);
 80 Dx= (D.*sin(theta))'*cos(phi);
 81 Dy= (D.*sin(theta))'*sin(phi);
 82 Dz= repmat((D.*cos(theta))',1,N);
 83 Dc= repmat((D)',1,N);
 84 if flag == 1
 85 figure; 
 86 surf(Dx,Dy,Dz,Dc);
 87 shading interp;
 88 light;
 89 lighting gouraud;
 90 colorbar;
 91 title('directivity');
 92 end
 93 
 94 
 95 % ND = 10*log10(D(1:round(N/2))/max(D(1:round(N/2))));
 96 ND = 10*log10(D/max(D));
 97 peak = find(ND==max(ND));
 98 NDD= abs(ND+3);
 99 left = find(NDD(1:peak)==min(NDD(1:peak)));
100 find((NDD)==min(NDD));
101 right = find(NDD(peak:end)==min(NDD(peak:end)))+ peak-1;
102 beamwidth = (theta(right)-theta(left))/pi*180;
103 if flag == 1
104 figure;
105 plot(theta/pi*180,ND);
106 title('E plane directivity pattern');
107 end
108 str = strcat('3dB beam width is  ',num2str(round(beamwidth*100)/100), '° ');
109 
110 E_dB=10*log10(abs(E0)/max(abs(E0)));
111 E_dB(E_dB<=-60)=-60;
112 if flag == 1
113 figure;
114 polar_dB(theta*180/pi,E_dB,-60,0,4);
115 title('Elevation plane normalized amplitude pattern (dB)','fontsize',16);
116 end
117 D_dB=10*log10(D);
118 D_dB(D_dB<=-60)=-60;
119 if flag == 1
120 figure;
121 polar_dB(theta*180/pi,D_dB,-60,max(D_dB),4);
122 title(strcat('Elevation plane directivity pattern (dB):',str),'fontsize',12);
123 end
124 %text(60,0,str)
125 
126 Rr = 2*Prad;
127 Rin = Rr/(sin(k*l/2)^2);
128 Xm=30*(2*si(k*l)+cos(k*l)*(2*si(k*l)-si(2*k*l))- ...
129           sin(k*l)*(2*ci(k*l)-ci(2*k*l)-ci(2*k*rad^2/l)));
130 Xin=Xm/(sin(k*l/2))^2;
131  
132 function [y]=si(x);
133 v=linspace(0,x/pi,500);
134 dv=v(2)-v(1);
135 y=pi*sum(sinc(v)*dv);
136 
137 function [y]=ci(x);
138 v=linspace(0,x/(2*pi),500);
139 dv=v(2)-v(1);
140 y1=2*pi*sum(sinc(v).*sin(pi*v)*dv);
141 y=.5772+log(x)-y1;
142

OPTIMISM, PASSION & HARDWORK