[Moved]: near field to far field transformation code error

Code PHP (brief) - [expand]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
lambda=0.1;
%----------- scanning area definition ------------------------------------
velikostX = 20*lambda; % size of the scanning area
velikostY = velikostX;
Nx = (velikostX/lambda)*2+3; % samples at x axis
x = linspace(-(velikostX/2), (velikostX/2), Nx);
dx = abs(x(2)-x(1));
Ny = Nx; % samples at y axis
y = linspace(-(velikostY/2), (velikostY/2), Ny);
r0 = 2*lambda; % distance of scanning surface from the AUT
 
 
%------- testing antenna (AUT) definition (nine-dipole antenna) ----------
Nd = 9; % number of dipoles
dd = lambda/2; % dipole spacing
d = linspace((-(Nd-1)/2)*dd, ((Nd-1)/2)*dd, Nd);
I = ones(1,Nd); % currents for dipoles (1 A)
l = lambda/4; % dipole lenght
Fd = j*(cos(k*l)*cos(theta)-cos(k*l))/(sin(theta)); % radiation function of symetric dipole
 
 
%----------- near-field samples computation -------------------------------
E = zeros(Nx,Ny);
Ehelp = zeros(1,Nd); % help vector
for m=1:Ny
for n=1:Nx
for i=1:Nd
E(n,m) = E(n,m) + 60*(cos(k*l*(x(n)/sqrt(r0^2+(y(m)-d(i))^2+x(n)^2)))-cos(k*l))/(sqrt(r0^2+(y(m)-d(i))^2)/sqrt(r0^2+(y(m)-d(i))^2+x(n)^2))*I(i)*(exp(-j*k*sqrt(r0^2+(y(m)-d(i))^2+x(n)^2))/sqrt(r0^2+(y(m)-d(i))^2+x(n)^2));
end
end
end
Emax = max(E);
maxE = max(Emax);
Eneardb = 20*log(abs(E/maxE));
figure(1); % 3D near-fiel distribution
surf(x/lambda,y/lambda,Eneardb);
title('Rozlozeni normovaneho modulu pole v blizke oblasti (z = 2*lambda)');
xlabel('y[lambda]');
ylabel('x[lambda]');
zlabel('abs(Enear)');
% depiction of E (module and phase) in the near-field for two orthogonal cuts
figure(2);
subplot(2,2,1);
plot(x/lambda,Eneardb(:,((Ny+1)/2)));
title('rez rovinou y = 0');
xlabel('x[lambda]');
ylabel('abs(Enear)');
subplot(2,2,2);
plot(y/lambda,Eneardb(((Nx+1)/2),:));
title('rez rovinou x = 0');
xlabel('y[lambda]');
ylabel('abs(Enear)');
subplot(2,2,3);
plot(x/lambda,angle(E(:,((Ny+1)/2)))/(2*pi)*360);
xlabel('x[lambda]');
ylabel('phase(Enear)');
subplot(2,2,4);
plot(y/lambda,angle(E(((Nx+1)/2),:))/(2*pi)*360);
xlabel('y[lambda]');
ylabel('phase(Enear)');

Code PHP (brief) - [expand]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
%----------- near-field to far-field transformation --------------------
A = zeros(t,t); % spatial spectrum computation
for m=1:t
for n=1:t
for h=1:Ny
for i=1:Nx
A(n,m) = A(n,m) + E(i,h)*exp(j*k*(x(i)*kx(n,m)+y(h)*ky(n,m)));
end
end
end
t-m
end
A = A/(Nx*Ny);
 
 
 
 
 
% ------ far-field computation (theoretical values for comparison) -----------
Evzdal = zeros(t,t);
for m=1:t
for n=2:(t-1)
Evzdal(m,n) = ((cos(k*l*cos(phi(n)))-cos(k*l)) / sin(phi(n)))*(sin(Nd/2*k*dd*cos(theta(m)+pi/2)))/(sin(1/2*k*dd*cos(theta(m)+pi/2)));
end
end
Evzdalmax = max(Evzdal);
maxEvzdal = max(Evzdalmax);
Evzdal = abs(Evzdal/maxEvzdal);
Evzdaldb = 20*log(Evzdal);
figure(3);
surf(phi/(2*pi)*360,theta/(2*pi)*360,Evzdaldb);
title('Vyzarovaci diagram ziskany vypoctem');
xlabel('phi[°]');
ylabel('theta[°]');
zlabel('abs(Efar)');
% -------- multiplying the angular spectrum by a wavenumber ------------------
for i=1:t
A(i,:) = kz.*A(i,:);
A(:,i) = rot90(kz).*A(:,i);
end
 
 
 
 
 
Amax = max(A);
maxA = max(Amax);
A = A/maxA;
F = abs(A);
Fdb = 20*log(F);
figure(4); % depiction of spatial spectrum for two orthogonal cuts
subplot(2,1,1);
stem(ky(:,45),abs(A(:,45)));
title('Rez prostorovym spektrem v rovine kx=0');
xlabel('ky');
ylabel('abs(A)');
subplot(2,1,2);
stem(kx(45,:),abs(A(45,:)));
title('Rez prostorovym spektrem v rovine ky=0');
xlabel('kx');
ylabel('abs(A)');
figure(6); % theoretical and transformed values comparison
subplot(2,1,1);
plot(theta(2:(t-1))/(2*pi)*360,Evzdaldb(2:(t-1),45),'b',theta(2:(t-1))/(2*pi)*360,Fdb(2:(t-1),45),'r');
xlabel('theta[°]');
ylabel('abs(Evzdal/Emax)');
subplot(2,1,2);
plot(theta(2:(t-1))/(2*pi)*360,Evzdaldb(45,2:(t-1)),'b',theta(2:(t-1))/(2*pi)*360,Fdb(45,2:(t-1)),'r');
xlabel('phi[°]');
ylabel('abs(Evzdal/Emax)');