% Definimos el paso y los dos vectores de posibles valores de la función.

delta = (2*pi)/50;
theta = 0.001:delta/2:pi+0.001;
phi = 0:delta:(2*pi);
n = [1, 2, 3, 4, 5, 10, 25, 50:93];
n(50<n) = 0;
[Theta, Phi] = meshgrid(theta, phi);
[M, N] = meshgrid(theta, n);

% Definimos la función a representar.

i = 1;

if i == 1
    funcion = (sin(2.5*pi*sin(Theta))./(2.5*pi*sin(Theta)))^2;
elseif i == 2
    funcion = (sin(7*pi*sin(Theta))./(7*pi*sin(Theta)))^2;
else
    funcion = (cos(Theta)).^N;
end

% Dividimos la función entre su valor máximo para normalizar entre 0 y 1.

funcion = funcion/max(max(funcion));

% Obtenemos las componentes cartesianas de la función.

x = funcion.*sin(Theta).*cos(Phi);
y = funcion.*sin(Theta).*sin(Phi);
z = funcion.*cos(Theta);

% Y representamos la función con los ejes normalizados entre -30 y 0 dB.

pantalla = get(0, 'ScreenSize');
ventana = [pantalla(3)/2-300, pantalla(4)/2-300, 600, 600];
figure('Name', 'Diagrama de radiación en 3D', 'Position', ventana);

surf(x, y, z);

axis equal;
xlim([-1, 1]);
ylim([-1, 1]);
zlim([-1, 1]);

xlabel('dB');
ylabel('dB');
zlabel('dB');

L1 = get(gca, 'XLim');
set(gca, 'XTick', linspace(L1(1), L1(2), 6));
L2 = get(gca, 'YLim');
set(gca, 'YTick', linspace(L2(1), L2(2), 6));
L3 = get(gca, 'ZLim');
set(gca, 'ZTick', linspace(L3(1), L3(2), 11));

set(gca, 'XTickLabel', {0:-12:-30; -24:12:0});
set(gca, 'YTickLabel', {0:-12:-30; -24:12:0});
set(gca, 'ZTickLabel', {0:-6:-30; -24:6:0});