Matlab - Función 3D encontrar su segundo valor máximo.

clc; clear; close all; format short

data=importdata('Sph3D.txt');                            % Data import

Sph3D=data.data(:,3);                                    % Extract column

Mm=reshape(Sph3D,181,360);                               % Convert column to matrix

Ma=Mm';

theta=linspace(0.01,pi,181);

phi=linspace(0,2*pi,360);

[THETA,PHI]=meshgrid(theta,phi);

dtheta=pi/181;

dphi=pi/360;

minx = 0.50; maxx = 2.23; dx = (maxx-minx).*rand(1,1) + minx;

miny = 1.83; maxy = 3.45; dy = (maxy-miny).*rand(1,1) + miny;

deltax=0;  deltay=0;

k=2*pi;                                                                 % Wave number

psix=(k*dx*sin(THETA).*cos(PHI))+deltay;

psiy=(k*dy*sin(THETA).*sin(PHI))+deltay;

AFy=exp(j*(0)*psiy)+exp(j*(1)*psiy)+exp(j*(2)*psiy);                                                    

AFx=exp(j*(0)*psix)+exp(j*(1)*psix)+exp(j*(2)*psix);

AF=exp(j*(0)*psiy).*exp(j*(0)*psix)+exp(j*(0)*psiy).*exp(j*(1)*psix)+exp(j*(0)*psiy).*exp(j*(2)*psix)+exp(j*(1)*psiy).*exp(j*(0)*psix)+exp(j*(0)*psiy).*exp(j*(1)*psix)+exp(j*(0)*psiy).*exp(j*(2)*psix)+exp(j*(2)*psiy).*exp(j*(0)*psix)+exp(j*(0)*psiy).*exp(j*(1)*psix)+exp(j*(0)*psiy).*exp(j*(2)*psix);

AF=AFx.*AFy;

AFmag=abs(AF);

% Directivity Calculation _____________________________________________

Utheta=AFmag.^2;

Prad=sum(sum(Utheta.*sin(THETA)*dtheta*dphi));

F=4*pi*Utheta/Prad;

D=(F).*abs(Ma);

DdB=10.*log10(D);

figure

 mesh(D);shading interp; colormap('default');

title('3D Directivity as function of \theta and \phi','FontSize',15);

xlabel('\theta [deg]','FontSize',10);

ylabel('\phi [deg]','FontSize',10);

zlabel('Directivity','FontSize',10);

[A1,B1]=max(max(D));

[c1,f1]=find(D==A1)

A1

hold on

plot3(f1(1),c1(1),A1,'*')

D(c1(1),f1(1))=-inf;

[A2,B2]=max(max(D));

[c2,f2]=find(D==A2)

A2

plot3(f2(1),c2(1),A2,'*')

 view(0,0)

grid on

hold off

>> practica

c1 =

   180

f1 =

   146

A1 =

  163.8984

c2 =

   181

f2 =

   146

A2 =

  163.5324