theta0 = pi/2;
elevation0 = pi/2 - theta0;
[Azimuth,R] = cart2pol(X,Y);
Azimuth = wrapTo2Pi(Azimuth);
phi0 = pi/4 + 9*pi/56;
arrAz = linspace(0, 2*pi, 2*Rez)';
arrMaxR = nan(numel(arrAz),1);
for iAz = 1:numel(arrAz)
idx = find(abs(Azimuth-arrAz(iAz)) < 2*pi/Rez);
if (~isempty(idx))
tmpR = R(idx);
arrMaxR(iAz) = max(tmpR(:));
end
end
figure()
plot(X,Y);
hold on; axis equal; axis tight;
plot(arrMaxR.*cos(arrAz),arrMaxR.*sin(arrAz),'.k');
xlabel('X'); ylabel('Y');
maxRInDirectionPhi0 = interp1(arrAz,arrMaxR,phi0,'pchip');
p1 = [0 0];
p2 = [ maxRInDirectionPhi0.*cos(phi0) maxRInDirectionPhi0.*sin(phi0)];
dp = p2-p1;
quiver(p1(1),p1(2),dp(1),dp(2),0,'MaxHeadSize',0.3,'color','k','linewidth',3)
title({'Radiation Pattern In X-Y Plane'; ...
['R(\theta=\pi/2,\phi_0=\pi/4 + 9\pi/56)=' num2str(maxRInDirectionPhi0)]});