//  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, fct, args, sgn, ytol, xbounds, nsub, meth, impr) 
// Intervals of positive or negative values (low level)
//
// Calling Sequence
// [xres, yres] = CL__detectSign(x, y, fct, args, sgn, ytol, meth, impr)
//
// Description
// <itemizedlist><listitem>
// <p>Computes the intervals where fct(x) is positive or negative. </p> 
// <p>x and y are 2 row vectors. If y is empty => computed. y can be computed by: </p>
// <p>y = fct(x, ind, args) [in unused). </p>
// <p></p>
// <p>The other arguments are: </p>
// <p>- <b>ytol</b>: Desired accuracy on f(x). </p>
// <p>- <b>meth</b>: Method used to finely compute the zeros of y(x): "s": secant, "d": dichotomy, 
// "ds": mix of secant, "d": dichotomy (see CL_fsolveb).</p>
// <p>- <b>impr</b>: True (%t) if the improved detection mechanism is active (see below).</p>
// <p></p>
// <p>The method consists in selecting the abscissas for which the sign of y changes, and then computing a more 
// exact value using CL_fsolveb. If the improved detection mechanism is selected, the function first detects when 
// y reaches a extremum, which avoids choosing a very small step. The improved detection works however if the 
// function is twice differentiable. </p> 
// </listitem></itemizedlist>
// <p></p>
// </listitem>
// <listitem>
// <p><b>Note:</b></p>
// <p>- The abscissae x must be strictly increasing.</p>
// <p>- The results may be degraded is the data are not continuous. </p>
// <p></p></listitem>
//
// Parameters
// x: Abscissae. (1xN)
// y: Ordinates, or []. (1xN)
// fct: y = fct(x, ind, args).  
// args: (any type) argument for fct. 
// sgn:  Wanted sign: 1: positive, -1: negative.  
// ytol: Accuracy on ordinates or empty (empty <=> default value).  
// xbounds: Abscissae intervals. (2xQ) 
// nsub: (integer) Number of sub-division points (>= 0).  
// meth: (string) Method for zero computation: "s", "d" or "ds". 
// impr: (boolean) True if improved detection mechanism is selected.  
// xres: intervals of abscissae where value is positive (2xP)
// yres: corresponding ordinates (should be close to zero) (2xP)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_detectZero
// CL_fsolveb
//
// Examples
// x = linspace(0, 10, 30); 
// 
// function [y] = myfunction(x, ind, args)
//   y = sin(x) - args; 
// endfunction
// 
// args = 0.99; 
// sgn = 1; 
// ytol = 1.e-6; 
// xbounds = [-%inf, %inf]; 
// nsub = 0; 
// meth = "ds"; 
// impr = %t; 
// xres = CL__detectSign(x, [], myfunction, args, sgn, ytol, xbounds, nsub, meth, impr) 

// =================================================
// Declarations
// =================================================
// factor for computation of ytol
ATOL_FACTOR = 1.e-8; 
RTOL_FACTOR = 1.e-4; 

// =================================================
// MAIN
// =================================================

// -------------------------------------------------
// Arguments Checking
// -------------------------------------------------
if (argn(2) <> 10)
  CL__error("Invalid number of input arguments (10 expected)"); 
end

if (sgn <> 1 & sgn <> -1)
  CL__error("Invalid value for argument sgn (1 or -1 expected)"); 
end

if (typeof(ytol) <> "constant" | size(ytol,"*") > 1)
  CL__error("Invalid size for argument ytol"); 
end

// xbounds: 0 or 2 rows
if (size(xbounds,1) <> 2 & size(xbounds,1) <> 0)
  CL__error("Invalid size for argument xbounds"); 
end

if (nsub <> round(nsub) | nsub < 0)
  CL__error("Invalid value for argument nsub"); 
end

if (meth <> "s" & meth <> "d" & meth <> "ds")
  CL__error("Invalid value for argument meth"); 
end

if (typeof(impr) <> "boolean")
  CL__error("Invalid type for argument impr"); 
end

if (size(x,1) > 1 | size(y,1) > 1 | (size(x,2) <> size(y,2) & y <> []))
  CL__error("Invalid size for argument x or y"); 
end

// Check x is strictly increasing
if (size(x, 2) >= 2)
  if (find(x(2:$) - x(1:$-1) <= 0) <> [])
    CL__error("Invalid abscissas: not strictly increasing"); 
  end
end

// check there is no Nan value 
// NB: works even if y is []
if (find(isnan(x+y)) <> []) 
  CL__error("Invalid data (data contain Nan values)"); 
end

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

// special cases: x contains less than 2 values
// If x == [] (=> y == []): finished
if (size(x, 2) <= 1 | xbounds == []) 
  return; // <= RETURN
end

// find all zeros
// NB: zeros belong to xbounds 
idirect = 0;
xzero = CL__detectZero(x, y, fct, args, idirect, ytol, xbounds, nsub, meth, impr);

// make xbounds in [x(1), x($)) + eliminates intersections 
// NB: includes validation of xbounds (beginning <= end)
xbounds = CL_intervInters([x(1); x($)], CL_intervUnion(xbounds)); 

xall = gsort(unique([xbounds(1,:), xbounds(2,:), xzero]), "c", "i"); 

// NB: xres: not empty 
xres = CL_intervInters(xbounds, [xall(1:$-1); xall(2:$)]); 

// select intervals with correct sign
ymid = fct((xres(1,:) + xres(2,:))/2, [], args); 
I = find(ymid * sgn > 0); 
xres = CL_intervUnion(xres(:,I));  

// generate y values if necessary
nbargs_out = argn(1); 
if (xres <> [] & nbargs_out >= 2)
  yres = matrix(fct([xres(1,:), xres(2,:)], [], args), 2, -1); 
end

endfunction