//  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] = CL_intervSelectIn(interv, x, opt)
// Values inside intervals
//
// Calling Sequence
// xres = CL_intervSelectIn(interv, x [, opt])
//
// Description
// <itemizedlist><listitem>
// <p>Returns values that fall in at least one interval.</p>
// <p>The returned values are sorted in increasing order. </p>
// <p>The argument <b>opt</b> can be:  </p>
// <p>- "d": (d as "default") only values in x that are inside intervals are returned. </p>
// <p>- "b": same as "d" + bounds of intervals (containing at least one value) are added. </p>
// <p>- "bs": same as "b" + "%nan" separators are added between intervals. </p>
// <p></p>
// </listitem>
// <listitem>
// <p><b>Notes:</b></p>
// <p> - Only the union of intervals is considered. The returned bounds are the bounds of these intervals
// and not of the input ones. </p>
// <p> - The output values are all different. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// interv: Set of intervals [start; end] (2xN)
// x: Set of values (1xP)
// opt: (string, optional) output option: "d", "b", "bs". Default is "d"
// xres: selected values (+ optional bounds and separators) (1xQ)
//
// See also
// CL_intervUnion
//
// Authors
// CNES - DCT/SB
//
// Examples
// interv = [[0.9; 2.9], [5.1; 6]];
// 
// // Example 1: 
// // Note: First bound is not included 
// x = 2 : 0.5 : 8; 
// CL_intervSelectIn(interv, x)
// CL_intervSelectIn(interv, x, opt="b")
// CL_intervSelectIn(interv, x, opt="bs")
// 
// // Example 2: 
// // Note: First bound is included 
// x = 0 : 0.5 : 8; 
// CL_intervSelectIn(interv, x)
// CL_intervSelectIn(interv, x, opt="b")
// CL_intervSelectIn(interv, x, opt="bs")

// Declarations:


// Code:
if (~exists("opt", "local")); opt = "d"; end

// NB: size 0: is accepted ([])
if (size(interv,1) <> 0 & size(interv,1) <> 2)
  CL__error("Invalid intervals"); 
end

// NB: size 0: is accepted ([])
if (size(x,1) > 1)
  CL__error("Invalid values for argument x"); 
end

if (opt <> "d" & opt <> "b" & opt <> "bs")
  CL__error("Invalid value for argument opt"); 
end

// If empty inputs: returns []
if (interv == [] | x == [])
  xres = []; 
  return; // <= RETURN
end

// union of intervals => no intersection
// (=> sorted so that lower bound is increasing)  
intervu = CL_intervUnion(interv); 

// bounds of intervals in increasing order
vals = gsort(matrix(intervu, 1, -1), "c", "i"); 

// search values: odd intervals (# 1, 3, ...) are OK
// NB: adjust dsearch results in include bounds  
// condition "k < N" added by security (should be always true)
k = dsearch(x, vals); // => in 1 .. N-1; 0 => not found
I = find(modulo(k,2) == 0 & x == vals(k+1) & k <> 0 & k < size(vals,2)); 
if (I <> []); k(I) = k(I) + 1; end

// ind: indices of intervals in intervu
// I: indices of values in intervals
I = find(modulo(k,2) == 1 & k <> 0); 
ind = unique((k(I)+1)/2); 

// solutions that belong to interv
xres = x(I); 

if (opt == "d")
  // sorted (and unique) values 
  xres = unique(xres);   
end

// add bounds
if (opt == "b" | opt == "bs")
  if (length(ind) > 1)
    xres = [xres, intervu(2,ind(1:$-1)), intervu(1,ind(2:$))];  
  end
  if (length(ind) >= 1)
    bmin = intervu(1,ind(1)); 
    bmax = intervu(2,ind($)); 
    if (min(x) <= bmin); xres = [xres, bmin]; end
    if (max(x) >= bmax); xres = [xres, bmax]; end
  end
  xres = unique(xres); 
end

// add %nan between intervals
if (opt == "bs")
  if (length(ind) > 1)
    valm = mean([intervu(2,ind(1:$-1)); intervu(1,ind(2:$))], "r"); 
    y = [xres; zeros(xres)]; 
    y = [y, [valm; ones(valm)]]; 
    y = CL_sortMat(y, y(1,:)); 
    I = find(y(2,:) == 1); 
    xres = y(1,:); 
    xres(1,I) = %nan; 
  end
end

endfunction