//  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, fct, args, direct, ytol, xbounds, nsub, meth, impr)
// Zero crossings computation (low level)
//
// Calling Sequence
// [xres, yres] = CL__detectZero(x, y, fct, args, direct, ytol, xbounds, nsub, meth, impr)
//
// Description
// <itemizedlist><listitem>
// <p>Computes zero crossings, that is abscissas x such that y(x) = 0. </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) [ind unused). </p>
// <p></p>
// <p>The other arguments are: </p>
// <p>- <b>direct</b>: Specifies the direction of y: 1=increasing, -1=decreasing, 0=both (default).</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 and dichotomy (see CL_fsolveb).</p>
// <p>- <b>impr</b>: True (%t) if the improved detection mechanism is active.</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 active, the function first detects when 
// y reaches a extremum, which improves the zero detection capability. 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. 
// direct: 1=increasing, -1=decreasing, 0=both.  
// 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: Abscissae of solutions (1xP)
// yres: Ordinates of solutions (should be close to zero) (1xP)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_detectSign
// CL_fsolveb
//
// Examples
// x = linspace(0, 10, 30); 
// 
// function [y] = myfunction(x,ind,args)
//   y = sin(x) - a; 
// endfunction
// 
// args = 0.99; 
// direct = "any"; 
// ytol = 1.e-6; 
// xbounds = [-%inf, %inf]; 
// nsub = 0; 
// meth = "ds"; 
// impr = %t; 
// xres = CL_detectZero(x, [], myfunction, args, direct, ytol, xbounds, nsub, meth, impr); 


// =================================================
// Declarations
// =================================================
// factor for computation of ytol
ATOL_FACTOR = 1.e-8; 
RTOL_FACTOR = 1.e-4; // not too small to avoid non-convergence


// generate simulation abscissae => in xbounds (including margin) 
// at least 2 points added at beginning/end of each interval (for CL__detectLocExtr) 
// %nan inserted between sequences => avoid detection of "false" extrema (in CL__detectLocExtr) or zeros
// interv: no intersection between one another and assumed included in [xref(1), xref($)]
// xref: reference abscissae (1xN) 
// interv: abscissa intervals (2xP)
// nsub: number of added intermediate abscissae (1x1)
function [x] = detectZero_genx(xref, interv, nsub)
	
	N = size(xref,2); 

  // i1: min index such that:  interv(2,:) <= xref(i1) 
	i1 = dsearch(interv(1,:), xref) + 1; // => in 2 .. N 
  I = find(interv(1,:) == xref(i1-1)); 
  if (I <> []); i1(I) = i1(I) - 1; end

  // i2: max index such that:  xref(i2) <= interv(1,:) 
	i2 = dsearch(interv(2,:), xref); // => in 1 .. N-1
  I = find(interv(2,:) == xref(i2+1)); 
  if (I <> []); i2(I) = i2(I) + 1; end
  
  nmargin = 2; 
  if (nsub > 0); nmargin = 1; end
  
	ind_interv = CL_intervUnion([max(1, i1 - nmargin); min(N, i2 + nmargin)]);
  
  ind = zeros(1, N+1); 
  ind(ind_interv(1,:)) = 1; 
  ind(ind_interv(2,:) + 1) = -1; 
  indcum = cumsum(ind); 
  
  x = %nan * ones(ind); 
  I = find(indcum > 0); 
  x(I) = xref(I); 
  
  I = find(indcum > 0 | ind == -1); 
  x = x(I);
  	
  if (size(x, 2) >= 2 & nsub > 0)
	  x_mat = CL_intervLinspace([x(1:$-1); x(2:$)], 2+nsub);  
	  x = [matrix(x_mat(1:$-1,:),1,-1), x_mat($,$)]; 
    I = find(isnan(x(1:$-1)) & isnan(x(2:$))); 
    x(I) = []; 
  end
  
endfunction


// select values in interv
// interv: no intersection between one another 
function [xres, yres] = detectZero_selectInterv(xres, yres, interv)
  val = matrix(interv, 1, -1); 
  ind = dsearch(xres, val); 
  I = find(ind <> 0 & modulo(ind, 2) == 0 & xres == val(ind+1)); 
  if (I <> []); ind(I) = ind(I) + 1; end
  I = find(ind <> 0 & modulo(ind, 2) == 1); 
  xres = xres(I); 
  yres = yres(I);
endfunction


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

if (direct <> round(direct) | abs(direct) > 1)
  CL__error("Invalid value for argument direct"); 
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 
if (find(isnan(x)) <> [] | find(isnan(y)) <> []) 
  CL__error("Invalid data (data contain Nan values)"); 
end

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

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

// special case: x contains 0 or 1 value or xbounds is []
if (size(x, 2) <= 1 | xbounds ==  []) 
  return; // <= RETURN
end

// computes abscissae if necessary (in xbounds, including margin) 
// => x may include Nan values
if (xbounds(1,1) > x(1) | xbounds(2,$) < x($) | nsub <> 0)
  x = detectZero_genx(x, xbounds, nsub); 
  y = []; // forces recomputation of y
end

// computes y if empty
// if x is %nan => y = %nan; 
if (y == [])
  y = %nan * ones(x); 
  I = find(~isnan(x)); 
  if (I <> []); y(I) = fct(x(I), [], args); end
end

// computes ytol if not initialized
if (ytol == [])
  ytol = min(max(abs(y)) * ATOL_FACTOR, (max(y)-min(y)) * RTOL_FACTOR); 
end


// -------------------------------------------------
// Zero detection
// -------------------------------------------------
// if improved detection: update data (local extrema)
if (impr & size(x, 2) >= 3)
  n = 5; 
  dytol = ytol / 5; 
  [xextr, yextr, I] = CL__detectLocExtr(x, y, fct, args, dytol, n); 
  if (I <> [])
    x(I) = xextr; 
    y(I) = yextr; 
  end
end

if (direct == 1) // increasing
  I = find(y(1:$-1) .* y(2:$) <= 0 & y(1:$-1) < y(2:$));
  
elseif (direct == -1) // decreasing
  I = find(y(1:$-1) .* y(2:$) <= 0 & y(1:$-1) > y(2:$));
  
else // if (direct == 0) 
  I = find(y(1:$-1) .* y(2:$) <= 0 & y(1:$-1) <> y(2:$));
end
  
if (I <> [])
  [xres, yres] = CL_fsolveb(fct, x(I), x(I+1), args, ytol=ytol, meth=meth, y1=y(I), y2=y(I+1));
end

// keep solution in xbounds only
[xres, yres] = detectZero_selectInterv(xres, yres, xbounds); 


endfunction