//  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'.


// Cuts a 2D curves into continuous segments. 
// (add %nan + intermediate values in x and y)
// x: array of x coordinates (typically longitude)
// y: array of y coordinates (typically latitude)
// xmin/xmax : target interval for x. 
// NB: x is supposed defined modulo xmod = xmax-xmin. 
// Algorithm: 
// x is first translated modulo xmod into [xmin,xmax]
// if |x(k+1) - x(k)| > xmod/2, intermediate points are added. 
// output: x1, y1 = continuous segments separated by %nan
// (always at least 1 %nan => x1 and y1 never empty)
// NB: also works if x and y already contain %nan

function [x1, y1] = CL__plot_genxy(x, y, xmin, xmax)

  // --------------------------------------------------
  // array insertion utility for inserting nb elmements 
  // after each element of index I.
  // I = set of indices in 1:n 
  // (increasing order, unique values, length(I) <= n)
  // nb = number of elements to insert after each index
  // Kn: indices where to insert "old" elements in new array
  // KI: indices of elements "I" in new array
  // (insertion positions in new array are: 
  // KI+1, KI+2 .. KI+nb
  // NB: new array has size = n + nb * length(I)
  // --------------------------------------------------
  function [Kn, KI] = insert_after_ind(n, I, nb); 
    // note1: n == 0 => Kn = [] 
    // note2: (0:min(0,n-1)) == 0 if n >= 1, else []
    d = zeros(1, n);  
    d(I) = nb; 
    Kn = (1:n) + cumsum([(0:min(0,n-1)), d(1:n-1)]);
    KI = (0 : length(I)-1) * nb + I; 
  endfunction


  if (xmax <= xmin)
     CL__error('xmax should be strictly greater than xmin');
  end

  if (size(x,1) > 1 | size(y,1) > 1)
     CL__error('wrong dimension for x or y: should be 1xN');
  end

  if (size(x,2) <> size(y,2))
     CL__error('wrong dimension for x or y: should have the same length');
  end

  // case : no data
  if (x == [])
    x1 = %nan; 
    y1 = %nan;
    return; // *** RETURN *** 
  end

  A = (xmax-xmin)/2; 
  x = CL_rMod(x, xmin, xmax); 

  delta_x = x(2:$) - x(1:$-1);  // OK if length(x) < 2 (-> [])  
  I = find(abs(delta_x) > A); // discontinuities
  nI = length(I); 
  
  // case : no discontinuity
  // special case as formulas that follows would not work 
  if (nI == 0)
    x1 = x; 
    y1 = y; 
    return; // *** RETURN *** 
  end


  // next point after discontinuity -> close to previous point
  // NOTE : X1 <> X(I) (in principle) 
  X1 = CL_rMod(x(I+1), x(I)-A, x(I)+A); 
 
  // 1st intermediate (added) point
  XA = xmin + (2*A) * round((x(I)-xmin)/(2*A));  // xmin or xmax 
  YA = y(I) + ((y(I+1)-y(I))./(X1-x(I))) .* (XA-x(I)); 

  // 2nd intermediate (added) point
  XB = xmin + xmax - XA; 
  YB = YA; 


  // built arrays with inserted values
  // NB: 3 values inserted at each discontinuity 
  // (including %nan)
  n = length(x); // same as length(y)
  x1 = %nan * ones(1, n + 3*nI);
  y1 = x1;

  // insert values

  [Kn, KI] = insert_after_ind(n, I, 3); 
 
  x1(Kn) = x; 
  x1(KI+1) = XA; 
  x1(KI+3) = XB; // NB: KI+2 => %nan

  y1(Kn) = y; 
  y1(KI+1) = YA; 
  y1(KI+3) = YB; // NB: KI+2 => %nan


endfunction