% Define theta and phi vectors
phi = [0 pi];
theta = linspace(-pi,pi,501).';
% Compute U
U(theta,phi)=eta/2((-1j*k*I*L*exp(-1j*k*r))/4*pi).^2 *(cos(theta)*cos(phi) + sin(phi)).^2
% Normalize U
U = U./max(U(:));
% Plot
figure;
polarplot(theta,U(:,1),'linewidth',2);
hold on;
polarplot(theta,U(:,2),'linewidth',2);
legend("\Phi = 0","\Phi = 180")