//  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_interpLin(xref,yref,x)
// Linear interpolation
//
// Calling Sequence
// [y] = CL_interpLin(xref,yref,x)
//
// Description
// <itemizedlist><listitem>
// <p>Linear interpolation.</p>
// <p>Given reference abscissae <b>(xref)</b> sorted in strictly increasing order and corresponding ordinates <b>(yref)</b>, 
// the function computes interpolated ordinates for the given abscissae <b>(x)</b>.</p>
// <p></p></listitem>
// <listitem>
// <p>Notes:</p>  
// <p> - y = %nan for values of x that are outside [xref(1),xref($)]</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// xref: Reference abscissae, sorted in strictly increasing order (1xN)
// yref: Corresponding ordinates (PxN)
// x: Abscissae, where to interpolate (1xM)
// y: Interpolated values (PxM)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_interpLagrange
//
// Examples
// // Example 1: Exact interpolation 
// xref = 1:6; 
// yref = [3*xref; -1+8*xref]; 
// x = 0:0.5:4; 
//
// [y] = CL_interpLin(xref, yref, x); 
// y - [3*x; -1+8*x]
//
// // Example 2: Interpolation of sin and cos
// xref = linspace(0,1,101); 
// yref = [sin(2*%pi*xref); cos(2*%pi*xref)];  
// x = linspace(0,1,1001); 
// 
// [y] = CL_interpLin(xref, yref, x); 
// scf();
// plot(x, y - [sin(2*%pi*x); cos(2*%pi*x)]);


N = size(xref, 2);
P = size(yref, 1);
M = size(x, 2);
y = %nan * ones(P,M);
  
if (size(xref,1) <> 1 | size(yref,2) <> N)
  CL__error("Wrong input argument size (xref or yref)"); 
end

if (M == 0) 
  y = []; 
  return; // <== RETURN
end

if (N < 2)
  CL__error("Insufficient number of values in xref (should be >= 2)"); 
end

if (find(xref(2:$)-xref(1:$-1) <= 0) <> [])
  CL__error("Abscissae not strictly increasing"); 
end



// y will be set to Nan for absissae outside [xref(1),xref($)]
y = %nan * ones(P,M);
K = find(x >= xref(1) & x <= xref($));
x = x(K);
  
// I = index such that xref(I) < x <= xref(I+1)
// Special case: if x = xref(1), then I = 1
I = dsearch(x, xref);
    
// Linear interpolation: y = y1 + (y2-y1) / (x2-x1) * (x-x1) 
// Do not create intermediate variable x1,x2,y1,y2 in order to save memory
y(:,K) = yref(:,I) + (yref(:,I+1)-yref(:,I)) ./       ...  // y1 + (y2-y1) /
                (ones(P,1) * (xref(I+1)-xref(I))) .*  ...  // (x2-x1) *
                (ones(P,1) * (x - xref(I)));               // (x-x1)

endfunction