function [n,l,w,amp,I]=select_modes(n,l,w,params);
% MATLAB HELP:
%
% function [n,l,w,amp,I]=select_modes(n,l,w,params);
%
% select the modes out of n,l,w that maximize some function,
% which is specified by the filter parameters in params
%
% inputs:
%  n, radial orders
%  l, angular degrees
%  w, frequencies
%  params, filter paramters, this structure is complicated, see codes.ps
% output:
%  n, radial orders of the selected modes
%  l, angular degrees of the selected modes
%  w, frequencies of the selected modes
%  amp, amplitude of the mode given the filter in params
%  I, index into the original mode list that tells which modes
%   this code has picked out, this is useful to know for later
%
%
% see codes.ps for details
%
% aaron birch may 21 2002






if (params.EXCLUDE_F)
 I=find(n~=0);
 n=n(I);
 l=l(I);
 w=w(I);
end;


if (params.FILTER_TYPE==0)
g=zeros(size(l));
for (the_n=0:max(n))
	I=find(n==the_n);
	if length(I>0)
         l1=l(I);
         w1=w(I);
         g1=gradient(w1,l1);
         g(I)=g1;
	end;
end;
speed=params.FILTER_PARAM1;
speed_WIDTH=params.FILTER_PARAM2;
amp=exp(-0.5*(g-speed).^2 / (speed_width)^2 );
end;


if (params.FILTER_TYPE==1)  %% use the phase speed filter
    speed=params.FILTER_PARAM1;
    speed_WIDTH=params.FILTER_PARAM2;
    p=w./(l+0.00001);  %% phase speed
    amp=exp(-0.5*(p-speed).^2 / (speed_width)^2 );
end;


if (params.FILTER_TYPE==2)  %% use the l filter
    amp=(l>params.FILTER_PARAM1).*(l<params.FILTER_PARAM2);     
end;

w_zero=params.FILTER_PARAM3;
w_width=params.FILTER_PARAM4;

amp=amp.*exp(-0.5*(w-w_zero).^2 / (w_width^2) );
[tmp,I]=sort(amp);

I=I(length(I)-params.NMODES+1:length(I));

n=n(I);
l=l(I);
w=w(I);
amp=amp(I);