function [x,ssq,cnt] = cxroot(FUN,x0,varargin)
% CXROOT Root of a complex function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The function finds a complex root of a complex function FUN by solving
% a system of two nonlinear real equations
% real(FUN(x)) = 0;
% imag(FUN(x)) = 0;
% of a complex variable x = complex(real(x), imag(x)) by calling the
% function LMFnlsq from MathWorks, File Exchange:
% http://www.mathworks.com/matlabcentral/fileexchange/17534
% ie. as a least squares problem without any need of Optimization Toolbox.
%
% Calls of the function:
% x = cxroot(FUN)
% x = cxroot(FUN,x0)
% x = cxroot(FUN,x0,options)
%
% FUN name or handle of the function the root of which is sought.
% The function itselve is either normal m-function, or a nested
% function. In the later case, constant parameters of the function
% may be transfered from the calling procedure without globals
% (see Example).
% x0 Initial guess of the complex root,
% varargin a set of optional names and values pairs of parameters,
% which control iterations in the function LMFnlsq (see its help)
% If no options are given, default ones are used.
% x complex root of the function FUN,
% ssq sum of squares of residuals (differences from zero),
% cnt number of evaluations of FUN and Jacobian matrix,
% negative if no convergence (see LMFnlsq)
%
% Example:
% function [w,ssq,cnt] = AireAiPrime(k,c1,w0,Options)
% % AIREAIPRIME Function for complex root search
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% [w,ssq,cnt] = cxroot(@funw,w0,Options{:});
% function fw = funw(w)
% Z = i*w*(6i*k)^(-2/3);
% fw = airy(k,Z)-c1*exp(i*pi/6);
% end % fw
% end % AireAiPrime
%
% The function finds a root (with ssq -> 0), which is close to initial
% guess. Sometimes it finishes iterations in points, which are not roots.
% Such a situation is recognized by a nonzero value of ssq.
% Results of a run of the script Airekc1 with w0 = 5+5i,
% Options = {'Display',-2} and parameters
% k = 1.0000 =>
% c1 = 1.0000 =>
% Re w0 = 5.0000 =>
% Im w0 = 5.0000 =>
%**************************************************************************
% itr nfJ SUM(r^2) x dx
%**************************************************************************
% 0 1 7.1377e-001 5.0000e+000 0.0000e+000
% 5.0000e+000 0.0000e+000
% 2 3 8.8108e-002 6.2365e+000 1.2514e+000
% 6.4646e+000 -1.2406e+000
% 4 5 4.9009e-007 6.7102e+000 -6.8413e-002
% 6.2243e+000 -2.0625e-002
%
% 6 7 1.9376e-027 6.7097e+000 -2.0633e-007
% 6.2256e+000 4.6777e-007
% w =
% 6.7097 + 6.2256i % ROOT
% ssq =
% 1.9376e-027 % SUM OF SQUARES
% cnt =
% 6 % NUMBER OF ITERATIONS
%
% It is obvious that w is a root, because of ssq, which dropped from almost
% 1 to almost zero.
% Miroslav Balda
% miroslav AT balda DOT cz
% 2009-01-04 v 1.0
%if nargin==0 && nargout==0, help cxroot, return, end % display help
Id = ''; % Test whether LMFnlsq is on the matlabpath
if ~exist('LMFnlsq.m','file'), Id=[Id ' LMFnlsq (FEX Id=17534)\n']; end
if ~isempty(Id)
error(['Download missing function\n' Id ' from File Exchange']);
end
if nargin<3, varargin = {''};
if nargin<2, x0=0; end
end
x = [real(x0); imag(x0)]; % decompose complex number into components
[x,ssq,cnt] = LMFnlsq(@res,x,varargin{:}); % solve nonlinear equations
x = complex(x(1),x(2));
% Nested function for evaluating residuals:
function r = res(x)
Fx = feval(FUN,complex(x(1),x(2)));
r = [real(Fx);imag(Fx)];
end % res
end % cxroot
