//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

function [xres, yres] = CL_detectZero(x, y, direct, ytol, xbounds, opts)
// Zero crossing detection
//
// Calling Sequence
// [xres, yres] = CL_detectZero(x, y [, direct, ytol, xbounds, opts])
//
// Description
// <itemizedlist><listitem>
// <p>Computes zero crossings, that is (abscissa,ordinate) pairs such that ordinate = 0. </p> 
// <p>There are 2 possibilities for the arguments <b>x</b> and <b>y</b>: </p>
// <p>1) x and y are 2 row vectors. Lagrange interpolation in then used to compute data at intermediate abscissae. </p>
// <p>2) x is a row vector and y is a list that gives the function to be called. y must be defined as follows: </p>
// <p>y = list(fct) => fct is a function with the interface: z = fct(x) </p>
// <p>y = list(fct, args) => fct is a function with the interface: z = fct(x, args) </p>
// <p>y = list(fct, ind, args) => fct is a function with the interface: z = fct(x, ind, args) [ind: unused] </p>
// <p></p>
// <p>The other arguments are: </p>
// <p>- <b>direct</b>: Specifies the direction of y: "incr"=increasing, "decr"=decreasing, "any"=both. </p>
// <p>- <b>ytol</b>: Desired accuracy on y. </p>
// <p>- <b>xbounds</b>: Abscissae intervals where the abscissae are looked-for. </p>
// <p>- <b>opts</b>: structure containing additional optional parameters. </p>
// <p></p>
// <p>The <b>opts</b> structure may contain the fields: </p>
// <p>- <b>nsub</b>: (integer) Number of intermediate points (>= 0). </p>
// <p>- <b>meth</b>: Method used to finely compute the zeros of y(x): "s": secant, "d": dichotomy, "ds": mix of
// secant and dichotomy (see <link linkend="CL_fsolveb">CL_fsolveb</link>). </p>
// <p>- <b>impr</b>: True (%t) to activate the improved detection mechanism. Default value is %t.</p>
// <p>- <b>ninterp</b>: Number of interpolation points. Default value is 8.</p>
// <p></p>
// <p>The method consists in selecting the abscissae for which the sign of y changes, and then computing a more 
// exact value using CL_fsolveb. If the improved detection mechanism is active, the function first detects when 
// y reaches an extremum, which improves the zero detection capability. </p> 
// <p></p></listitem>
// <listitem>
// <p><b>Notes:</b></p>
// <p>- The abscissae x must be strictly increasing.</p>
// <p>- The data must be continuous and smooth enough so that polynomial interpolation or zero detection works. </p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// x: Abscissae. (1xN)
// y: (1xN or list) Corresponding ordinates or list(function, ...). See above for details. 
// direct: (string, optional) "incr"=increasing, "decr"=decreasing, "any"=both. Default value is "any". 
// ytol: (optional) Accuracy on y (empty <=> default value). Default value <= max(abs(y)) / 1.e8. 
// xbounds: (optional) Abscissae intervals where the results are looked-for. Default value is [-%inf, %inf]. (2xQ)
// opts: (structure, optional) Additional options - see above
// xres: Abscissae of solutions. (1xP)
// yres: Ordinates of solutions. (1xP)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_detectSign
// CL_fsolveb
//
// Examples
// x = 0 : 10;  
// a = 0.99; 
// 
// // Example 1: Data interpolation 
// y = sin(x) - a; 
// xres = CL_detectZero(x,y); 
// sin(xres) - a  // => ~0 (includes effect of interpolation errors)
// 
// // Example 2: User-defined function 
// function [y] = myfunction(x,c)
//   y = sin(x) - c; 
// endfunction
// xres = CL_detectZero(x, list(myfunction, a)); 
// myfunction(xres,a)  // => 0 (more accurate) 


// -------------------------------------------------
// Function for zero detection => interpolation of input data
// args contains: x, y, ninterp
// NB: ind: not used but necessary for CL_fsolveb
// -------------------------------------------------
function [y] = detectZero_fct1(x, ind, args)
  y = CL_interpLagrange(args.x, args.y, x, args.ninterp);
endfunction

// -------------------------------------------------
// Function for zero detection => interface to user defined function
// args contains: fct (interface: fct(x))
// NB: ind: not used but necessary for CL_fsolveb
// -------------------------------------------------
function [y] = detectZero_fct2(x, ind, args)
  if (args.nb == 1)
    y = args.f(x);
  elseif (args.nb == 2)
    y = args.f(x, args.a);
  else // if (args.nb == 3)
    y = args.f(x, [], args.a);
  end
endfunction

// -------------------------------------------------
// Arguments Checking
// (additional checks done in called function)
// -------------------------------------------------
if (argn(2) < 2)
  CL__error("Invalid number of input arguments (at least 2 expected)"); 
end

if (~exists("direct", "local")); direct = "any"; end
if (~exists("ytol", "local")); ytol = []; end  // default value
if (~exists("xbounds", "local")); xbounds = [-%inf; %inf]; end  
if (~exists("opts", "local")); opts = struct(); end

// check that all field names in opts are correct
Names = ["nsub", "meth", "impr", "ninterp"]; 
optsfields = fieldnames(opts); 
if (setdiff(optsfields, Names) <> [])
  CL__error("Invalid fields in opts structure"); 
end

if (~isfield(opts, "nsub")); opts.nsub = 0; end  
if (~isfield(opts, "meth")); opts.meth = "ds"; end
if (~isfield(opts, "impr")); opts.impr = %t; end
if (~isfield(opts, "ninterp")); opts.ninterp = 8; end

if (direct <> "incr" & direct <> "decr" & direct <> "any")
  CL__error("Invalid value for argument direct"); 
end

typeofy = typeof(y); 
if (typeofy <> "constant" & typeofy <> "list")
  CL__error("Invalid type for argument y"); 
end

if (typeofy == "list")
  if (lstsize(y) < 1 | lstsize(y) > 3)
    CL__error("Invalid size for argument y (list of 1 to 3 elements expected)"); 
  end
  // "fptr": Scilab functions, "function": user-defined functions
  if (typeof(y(1)) <> "function" & typeof(y(1)) <> "fptr")
    CL__error("Invalid values in argument y (function expected)"); 
  end
end

// -------------------------------------------------
// Data preparation 
// -------------------------------------------------
// Initializes function for zero detection + associated arguments
Args = struct(); 
Fct = []; // will be a function ! 

if (typeofy == "constant")
  // user provides data to be interpolated
  Fct = detectZero_fct1; 
  Args.x = x; 
  Args.y = y; 
  Args.ninterp = opts.ninterp; 

else // if (typeofy == "list")
  // user provides function
  nb = lstsize(y); 
  if (nb <= 2)
    Fct = detectZero_fct2; 
    Args.f = y(1); 
    Args.nb = nb;  
    if (nb == 2)
      Args.a = y(nb); 
    end
  else // if (nb == 3)
    // optimization: user function called directly
    Fct = y(1);
    Args = y(3); 
  end
  
  // type of y changed! (y will be computed) 
  y = []; 
end

// -------------------------------------------------
// Computation
// -------------------------------------------------
// convert direction
if (direct == "incr")
  idirect = 1; 
elseif (direct == "decr")
  idirect = -1; 
else
  idirect = 0; 
end

[xres, yres] = CL__detectZero(x, y, Fct, Args, idirect, ytol, xbounds, opts.nsub, opts.meth, opts.impr);

endfunction