Difficulty plotting the polar plot of the electric field (theta component) radiated by a horn antenna in MATLAB.

I've implemented the following equation in MATLAB after defining all the equations involved in it. I've plotted the polar plot of the equation. But my graph is not coming out to be as expected. Can somebody help me plot the reqired graph?

The expected graph is this (The grey one):

I've written the following piece of code.

    k = 2*pi;
    theta = 0:0.01:2*pi;
    ky = k*sin(theta);
    a1 = 5.5;
    rho2 = 6;
    b = 0.25;
    kx_p = pi/a1;
    kx_pp = -pi/a1;
    t1_p = sqrt(1/(pi*k*rho2))*((-k*a1/2)-(kx_p*rho2));
    t2_p = sqrt(1/(pi*k*rho2))*((k*a1/2)-(kx_p*rho2));
    t1_pp = sqrt(1/(pi*k*rho2))*((-k*a1/2)-(kx_pp*rho2));
    t2_pp = sqrt(1/(pi*k*rho2))*((k*a1/2)-(kx_pp*rho2));
    f1 = ((kx_p^2)*rho2)/(2*k);
    f2 = ((kx_pp^2)*rho2)/(2*k);
    Y = (k*b/2).*sin(theta);
    F = @(t1, t2) (fresnelc(t2) - fresnelc(t1)) - 1j.*(fresnels(t2) - fresnels(t1));
    E_theta = (1j*(b/8)*((sqrt(k*rho2))/pi)*(exp(-1j*k))).*... %assuming E_2 = 1
    ((1 + cos(theta)).*sinc(Y).*((exp(1j*f1).*F(t1_p, t2_p)) + (exp(1j*f1).*F(t1_pp, t2_pp))));
    polarplot(abs(E_theta))

Answer 1

[Vinod Dwarapudi]

Hi Muhammad Hamza Shashid,

The H-plane expression needs to be implemented with angle-dependent spectral terms (kx_p, kx_pp) and corresponding Fresnel integration limits. Also verify that the second phase term uses f2 rather than f1. The code below reproduces the expected pattern from Balanis.

%% H-Plane Sectoral Horn Antenna Radiation Pattern
% E-plane and H-plane patterns in dB (polar plot)
% Based on Balanis, Antenna Theory

% Parameters
k = 2*pi;                   % Wavenumber (lambda = 1)
theta = 0:0.01:2*pi;        % Observation angles
a1 = 5.5;                   % Horn aperture width (H-plane, wavelengths)
b1 = 0.25;                  % Horn aperture height (E-plane, wavelengths)
rho2 = 6;                   % Horn slant length (wavelengths)
dB_floor = -40;

% Fresnel integral combination
F = @(t1, t2) (fresnelc(t2) - fresnelc(t1)) - 1j*(fresnels(t2) - fresnels(t1));

%% E-Plane (phi = 90 deg cut, uniform aperture)
Y = (k*b1/2) .* sin(theta);
E_E = (1 + cos(theta)) .* sinc(Y/pi);

%% H-Plane (phi = 0 deg cut, quadratic phase aperture)
kx_p = k*sin(theta) + pi/a1;
kx_pp = k*sin(theta) - pi/a1;

t1_p = sqrt(1/(pi*k*rho2)) .* ((-k*a1/2) - kx_p*rho2);
t2_p = sqrt(1/(pi*k*rho2)) .* ((k*a1/2) - kx_p*rho2);
t1_pp = sqrt(1/(pi*k*rho2)) .* ((-k*a1/2) - kx_pp*rho2);
t2_pp = sqrt(1/(pi*k*rho2)) .* ((k*a1/2) - kx_pp*rho2);

f1 = (kx_p.^2 * rho2) / (2*k);
f2 = (kx_pp.^2 * rho2) / (2*k);

E_H = (1 + cos(theta)) .* (exp(1j*f1).*F(t1_p, t2_p) + exp(1j*f2).*F(t1_pp, t2_pp));

%% Convert to dB and plot
E_E_dB = max(20*log10(abs(E_E)/max(abs(E_E))), dB_floor) - dB_floor;
E_H_dB = max(20*log10(abs(E_H)/max(abs(E_H))), dB_floor) - dB_floor;

figure;
polarplot(theta, E_E_dB, 'b-', 'LineWidth', 1.5)
hold on
polarplot(theta, E_H_dB, 'r-', 'LineWidth', 1.5)
hold off

ax = gca;
ax.ThetaZeroLocation = 'top';
ax.ThetaTick = 0:30:330;
ax.ThetaTickLabel = {'0','30','60','90','120','150','180','150','120','90','60','30'};
title('Amplitude Pattern (dB)')
legend('E-Plane', 'H-Plane', 'Location', 'southoutside')
rlim([0 abs(dB_floor)])
rticks(0:10:abs(dB_floor))
rticklabels(string((0:10:abs(dB_floor)) + dB_floor))