function [fr_pw]=find_lmax_iono(sxx,fr,fr_int,sx_mm,num_bands)
% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% non-working draft 2017/09/21
% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% Finds local maxima and selectrs those
% meeting criteria
% sxx is function, fr is argument
% fr_int is frequency interval where maxinum is selected
% sx_mm is interval (max/min) for PSD in maximum
% fr_hw is halfwidth to compare P_max and P_av
% psd_rat is threshold for P_max/P_av
% n_mx is an integer array with maxima positions
% p_mx is a real array with PSD values in maxima
% i_c is number of curve
% fr_prev and df_prev are previous values of x_max position and diff(x_max)
fr_mx=[];
p_mx=[];
fr_pw=zeros(num_bands,3)-1e11;

L_sxx=length(sxx);

sxx_l=sxx(1:L_sxx-2);
sxx_c=sxx(2:L_sxx-1);
sxx_r=sxx(3:L_sxx);

% positions of local maxima
f_max=find(sxx_c>sxx_l & sxx_c>sxx_r)+1;
L_max=length(f_max);
%pause
if L_max>0
fr_max=fr(f_max);

p_max=sxx(f_max);

%pause
find_ok_0=find(fr_max>=fr_int(1) & fr_max<=fr_int(2) & p_max<sx_mm(1) & p_max>=sx_mm(2));
fr_ok_0=fr_max(find_ok_0);
p_ok_0=p_max(find_ok_0);
n_ok_0=f_max(find_ok_0);
L_mx=length(find_ok_0);
if L_mx>0
fr_pw_u=[fr_ok_0',p_ok_0',n_ok_0'];
fr_pw_s=sortrows(fr_pw_u,2);
if L_mx>num_bands
    fr_pw_h=fr_pw_s(L_mx-num_bands+1:L_mx,:);
else
    fr_pw_h=fr_pw_s;
end
fr_pw_hs=sortrows(fr_pw_h,1);
if L_mx>=num_bands
    fr_pw=fr_pw_hs;
else
    fr_pw(1:L_mx,:)=fr_pw_hs;
end
end
end
    
  