[rts,frts,rtsm,frtsm,func]=interact(func,xlim,ylim)

function [rts,frts,rtsm,frtsm,func]=interact(func,xlim,ylim)
% [rts,frts,rtsm,frtsm,func]=interact(func,xlim,ylim)
% Function interact searches for complex function roots in a 
% rectangular region. Surface and contour plots of the function
% -log(abs(func(z))are created to indicate root locations within
% the search region. Roots are indicated by upward pointing 
% spikes, whereas poles produce downward pointing spikes. Root
% estimates selected interactively using a mouse are refined
% with a routine based on Muller's method. Data input can be
% made using the argument list or through interactive prompts
% when no call list parameters are given (when nargin==0).  
% func - the function for which roots are sought. This function
%        is typically expressed in anonymous form typified by
%        @(z) z.*besselh(0,1,z)./(1+z). To facilitate plotting,
%        the function to be analyzed should be vectorized to
%        accept a matrix argument.           
% xlim - a vector specifying the x limits of search. The vector
%        can have forms like [10,20], 10:.5:20, linspace(10,20,100)
% ylim - a vector specifying the y limits of search. The vector
%        forms are the same as those for xlim.

%              Howard Wilson, August, 2009
close all, format long
if nargin==0, disp(' ')
  disp('INTERACTIVE ROOT COMPUTATION USING MULLER''S METHOD')
  func=input('Enter a function handle: ');
  u=input('Enter search limits xmin, xmax, or an x vector: ','s');
  xlim=eval(['[',u,']']); 
  u=input('Enter search limits ymin, ymax, or a  y vector: ','s');
  ylim=eval(['[',u,']']);   
end
if length(xlim)==2, xlim=linspace(xlim(1),xlim(2),80); end
if length(ylim)==2, ylim=linspace(ylim(1),ylim(2),80); end
[x,y]=meshgrid(xlim,ylim); z=x+i*y; 
f=feval(func,z); f=-log(scale(f));
if nargin==0, disp(' ')
  disp('UPWARD spikes on a surface plot indicate roots.')
  disp('DOWNWARD spikes indicate poles and should not be')
  disp('selected. Press ENTER to show the surface. Press')
  disp('ENTER again to show a contour plot. Approximate') 
  disp('roots are selected by left clicking the mouse.')
  disp('Press ENTER once more to end the root selection.')
  pause
end
fmin=min(f(:)); fmax=max(f(:)); fdif=fmax-fmin;
%figure
surf(x,y,f),  caxis([fmin-.25*fdif,fmax]); %colormap([1 1 0])  
xlabel('x axis'), ylabel('y axis'), zlabel('-log(abs(func(z)))')
title(['-log( abs( ',func2str(func),' ) )']) 
shift, shg, pause, close
contour(x,y,f,50), colormap default; 
str1=['Contours of  -log( abs ( ',func2str(func),' ) )'];
str2='LEFT CLICK WITH THE MOUSE TO SELECT APPROXIMATE ROOTS';
title({str1,str2})
xlabel('x axis'), ylabel('y axis')
shift 
[xm,ym]=ginput; 
close
if isempty(xm)
  disp(' '), disp('NO ROOTS WERE COMPUTED'), disp(' ')
  rtsm=[]; frtsm=[]; rts=[]; frts=[]; return
end
rtsm=xm+i*ym; frtsm=feval(func,rtsm);  
[rts,frts]=muller(func,0,500,0,0,0,rtsm);
summary(rts,frts,rtsm,frtsm)


function shift(left,bottom,width,height)
% shift(left,bottom,width,height) moves the graphics window
if nargin==0
  left=.52; bottom=.38; width=.7/1.6; height=.6;
end
set(gcf,'units','normalized','outerposition',...
    [left,bottom,width,height])

function summary(rfinl,ffinl,rinit,finit)
% summary(rfinl,ffinl,rinit,finit)
format long, disp(' ')
disp('Final Roots Corrected by Muller''s Method')
disp(rfinl)
disp('Largest Magnitude of Root Correction')
disp(max(abs(rfinl-rinit)))
disp('Largest Magnitude of Final Function Values')
disp(max(abs(ffinl)))

function ha=scale(h)
% ha=scale(h) replaces infinite and zero values in h to allow 
% surf and contour plots of the function -log(abs(h))
s=size(h); ha=abs(h(:)); 
kfin=find(isfinite(ha)); % find finite values
kinf=find(isinf(ha));    % find infinite values
kzero=find(ha==0);       % find zero values
% eliminate infinite and zero values
hsrt=sort(ha(kfin)); kpos=find(hsrt>0);
small=hsrt(kpos(1))/10;  % smallest nonzero value
big=hsrt(end);           % largest non-infinite value
ha(kinf)=big;            % remove infinite values
ha(kzero)=small;         % remove zero values
ha=reshape(ha,s);