//  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 [y] = CL_unMod(y, a, x, slope)
// Makes data continuous (removes "modulo")
//
// Calling Sequence
// yc = CL_unMod(y, a, [, x, slope])
//
// Description
// <itemizedlist><listitem>
// <p>Makes data continuous by adding an integer number of "a" to data when necessary. </p>
// <p>Abscissae and slope(s) may be provided in case the "distance" between consecutive
// data is too large (more than a/2). </p> 
// <p>Notes: </p>
// <p>- y may contain %nan: %nan values are ignored. </p>
// <p>- The slope does not have to be exact. It only has to be precise enough so that the  
// prediction error between consecutive data is small enough. </p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// y: Vector of real values (ordinates) (1xN)
// a: Positive value (= "modulo")  (1x1)
// x: (optional) Vector of real values (abscissae). Default is []. (1xN)
// slope: (optional) Approximate value of the slope = dy/dx. Default is []. (1x1 or 1xN)
// yc : y made continuous (1xN)
//
// Authors
// CNES - DCT/SB
//
// Examples
// // x: finely sampled
// x = 0:0.2:30; 
// y = CL_rMod(x, 0, 2*%pi); 
// CL_unMod(y, 2*%pi) - x  // => 0
// 
// // x: loosely sampled
// x = x(1:20:$);
// y = y(1:20:$);
// CL_unMod(y, 2*%pi, x, 1) - x // => 0

// Declarations:

// Internal function
// Hypothesis: y does not contain %nan
// a = scalar, other vectors: same size. 
function [y] = unMod(y, a, x, slope)
  if (size(y,2) < 2)
    return; // nothing to do
  end
  
  if (slope == [])
    // Detects jumps (jump are between indices I and I+1)
    I = find(abs(y(:,2:$) - y(:,1:$-1)) > a/2); 
    delta = zeros(y); 
    delta(I+1) = -round((y(I+1)-y(I))/a) * a;
    y = y + cumsum(delta);  
  else
    // adjust each "point" (except first one)
    I = 1:length(x)-1; 
    d = -y(I+1) + y(I) + (x(I+1)-x(I)) .* slope(I); 
    delta = zeros(y); 
    delta(I+1) = round(d/a) * a;
    y = y + cumsum(delta);  
  end
endfunction


// Code:
if (~exists("x", "local")); x = []; end
if (~exists("slope", "local")); slope = []; end

if (a <= 0)
  CL__error("Invalid value for ''a'' (should be >0)"); 
end

[y, x, slope] = CL__checkInputs(y, 1, x, 1, slope, 1); 

// x and slope must be both empty or not empty
if ((x == [] & slope <> []) | (x <> [] & slope == []))
  CL__error("Inconsistent values for ''x'' and ''slope''"); 
end

// check that either x is [] or x is strictly increasing
if (x <> [])
  if (find(x(:,2:$) - x(:,1:$-1) <= 0) <> [])
    CL__error("Invalid value for ''x'': not strictly increasing");
  end
end

// makes continuous without %nan
// and restore initial vector size. 
// NB: code optimized for the case there are no %nan
// (also works if all elements are %nan)

// nanex: %t if there is at least one %nan
nanex = (find(isnan(y)) <> []); 

if (nanex)
  Innan = find(~isnan(y)); // not %nan
  N = size(y,2); 

  y = y(Innan); 
  if (x <> [])
    x = x(Innan);
    slope = slope(Innan); 
  end
end

[y] = unMod(y, a, x, slope); 

if (nanex)
  ynnan = y; 
  y = %nan * ones(1,N); 
  y(Innan) = ynnan; 
end


endfunction