//  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_detectSign(x, y, sgn, ytol, xbounds, opts)
// Intervals of positive or negative values 
//
// Calling Sequence
// [xres, yres] = CL_detectSign(x, y, sgn [, ytol, xbounds, opts])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the abscissa intervals where ordinate values are positive or negative. </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>sgn</b>: "+" => intervals where y is positive, "-" => intervals where y is negative.</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. 
// sgn: (string) Wanted sign for y: "+" or "-".  
// 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: Intervals of abscissae where y has the requested sign. (2xP)
// yres: Corresponding ordinates. (2xP)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_detectZero
// CL_fsolveb
//
// Examples
// x = 0 : 10; 
// a = 0.99; 
// 
// // Example 1: Data interpolation 
// y = sin(x) - a; 
// xres = CL_detectSign(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_detectSign(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] = detectSign_fct1(x, ind, args)
  y = CL_interpLagrange(args.x, args.y, x, args.ninterp);
endfunction

// -------------------------------------------------
// Function for zero detection => interface to user defined function - 1
// args contains: fct (+ nb + a if necessary; nb = number of arguments)
// NB: ind: not used but necessary for CL_fsolveb
// -------------------------------------------------
function [y] = detectSign_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) < 3)
  CL__error("Invalid number of input arguments (at least 3 expected)"); 
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 (sgn <> "+" & sgn <> "-")
  CL__error("Invalid value for argument sgn"); 
end

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

if (size(x,1) > 1)
  CL__error("Invalid size for argument x"); 
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
  if (typeof(y(1)) <> "function" & typeof(y(1)) <> "fptr")
    CL__error("Invalid values in argument y (function expected)"); 
  end
end

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

if (typeofy == "constant")
  // user provides data to be interpolated
  Fct = detectSign_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 = detectSign_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
  
  // y will be computed (and type change!) 
  y = []; 

end

// -------------------------------------------------
// Computation
// -------------------------------------------------
xres = [];
yres = []; 

xbounds = CL_intervInters([x(1); x($)], CL_intervUnion(xbounds)); 

if (xbounds == [])
  return; // <= RETURN 
end

// convert sign
if (sgn == "+")
  isgn = 1; 
else
  isgn = -1; 
end

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