%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Math 3520: Differential Equations and Dynamical Systems
% Compute local maximums and local maximum values using
% quadratic interpolation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%% Find local arg-maxima and maximum values:
function [Xs,Ys] = lmax(t,z)
% Input
%  t = Row vector of independent variable samples
%  z = row vector of dependent variable samples
% Output
%  X = row vector of local maxima
%  Y = row vector local maximum values
 dz = diff(z);  % differences z(k+1)-z(k), k=1,2,...
 n = 0;         % initial local max index
 for k=1:(length(dz)-1)
   if ( dz(k)>0 && dz(k+1)<0 )  % found local max at k+1
      n = n+1;      % increment the count of local maxes
      x1=t(k); x2=t(k+1); x3=t(k+2);
      y1=z(k); y2=z(k+1); y3=z(k+2);
      xnum=x1*x1*(y3-y2)+x2*x2*(y1-y3)+x3*x3*(y2-y1);
      xden=x1*(y3-y2)+x2*(y1-y3)+x3*(y2-y1);
      Xs(n) = 0.5*xnum/xden;  % append nth arg max
      d12=x1-x2; d31=x3-x1; d23=x2-x3;
      ynum1= 2*(y1*y1*d23^4+y2*y2*d31^4+y3*y3*d12^4);
      ynum2= (y1*d23^2+y2*d31^2+y3*d12^2)^2;
      yden = 4*d12*d23*d31*(y1*d23+y2*d31+y3*d12);
      Ys(n) = (ynum1-ynum2)/yden; % append nth max val
   end  % endif
 end  % endfor
 return
end % endfunction