//  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 [x, y, s, xdot, ydot, sdot] = CL__iers_xys_interp(jd, cder)
// X,Y,s with interpolation in pre-computed values
//
// Calling Sequence
// [x, y, s, xdot, ydot, sdot] = CL__iers_xys_interp(jd [,cder])
//
// Description
// <itemizedlist><listitem>
// <p>This function returns x,y,s values interpolated using pre-computed x,y,s values</p>
// <p>The pre-computed x,y,s values were obtained using IAU2006 precession and IAU2000AR06 nutation, using classical angles.</p>
// <p>Inteprolation parameters: time step: 12h, degree of polynomial: 5, time span includes [2000/1/1 0h, 2100/1/1 0h]</p>
// <p>The interpolation ensures a precision of 1 microarcsecond</p>
// <p>If a date is outside the time span, the IAU2006 precession and IAU2000AR06 nutation models (using classical angles) are used.</p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// jd: Two-part julian day (Time scale: TT) (2xN)
// cder: (boolean, optional) Option to compute derivatives. If cder is %f, the derivatives will be set to []. Default is %t. (1x1)
// x: X coordinate of CIP [rad] (1xN)
// y: Y coordinate of CIP [rad] (1xN)
// s: CIO locator [rad] (1xN)
// xdot: (optional) Time derivative of X coordinate of CIP [rad/s] (1xN)
// ydot: (optional) Time derivative of Y coordinate of CIP [rad/s] (1xN)
// sdot: (optional) Time derivative of CIO locator [rad/s] (1xN)
//
// Authors
// CNES - DCT/SB

// Declarations:

// Internal function
// xys and xysdot are structures containing :
// - tmin : the minimum date on which interpolation is performed = 1 jan 2000 0h TT
// - tmax : the maximum date on which interpolation is performed = 1 jan 2100 0h TT
// - tref : the array (1xNref) of reference dates
// - xref (or xdotref) : the array (1xNref) of reference X values (or the time derivative of X)
// - yref (or ydotref) : the array (1xNref) of reference Y values (or the time derivative of Y)
// - sref (or sdotref) : the array (1xNref) of reference s values (or the time derivative of Z)
// Note: 
// The dates are "modified julian dates" from 2000/1/1 0h, in TT time scale.

function [xys, xysdot] = CL__iers_xys_load(cder)
  // NB : argument cder is mandatory in this subfunction !
  
  fname = fullfile(CL_home(), "data", "frame", "iers_xys_2006_2000AR06_classic.dat");
  
  if (isfile(fname))
    load(fname,"xys");
  else
    CL__error("File: " + fname + ": not found"); 
  end
  
  xysdot = [];
  if (cder)
    fname = fullfile(CL_home(), "data", "frame", "iers_xysdot_2006_2000AR06_classic.dat");
    
    if (isfile(fname))
      load(fname,"xysdot");
    else
      CL__error("File: " + fname + ": not found"); 
    end
  end

endfunction

// Internal function
// Interpolates x,y,s and optionnaly xdot,ydot,sdot at dates t.
// t: Modified julian date from 1 jan 2000 0h (TT time scale)
// xys, xysdot: structures containing pre-computed results (see internal function above)
function [x, y, s, xdot, ydot, sdot] = CL__iers_xys_interp1(t, xys, xysdot, cder)
  // Number of points for the interpolation (The degree of the polynomial is Ni-1)
  Ni = 6;
  
  [xys_interp] = CL_interpLagrange(xys.tref, [xys.xref;xys.yref;xys.sref], t, Ni, opt=1);
  x = xys_interp(1,:);
  y = xys_interp(2,:);
  s = xys_interp(3,:);
  
  xdot = [];
  ydot = [];
  sdot = [];
  if (cder)
    [xys_interp] = CL_interpLagrange(xysdot.tref, [xysdot.xdotref;xysdot.ydotref;xysdot.sdotref], t, Ni, opt=1);
    xdot = xys_interp(1,:);
    ydot = xys_interp(2,:);
    sdot = xys_interp(3,:);
  end

endfunction


// Code:

if (~exists("cder","local")); cder = %t; end;
if (argn(1) <= 3); cder = %f; end;

N = size(jd,2);

// t = modified julian date from 1 jan 2000 0h TT
t = (jd(1,:) - 2451544) + jd(2,:) - 0.5;

x = zeros(1,N);
y = zeros(1,N);
s = zeros(1,N);
xdot = [];
ydot = [];
sdot = [];

// Load pre-computed X,Y,s values
[xys, xysdot] = CL__iers_xys_load(cder);

t_interp = (t >= xys.tmin & t <= xys.tmax);

// Compute only X,Y,s
if (~cder)
  
  // Interpolation for dates inside pre-computed values
  I = find(t_interp);
  [x(I), y(I), s(I)] = CL__iers_xys_interp1(t(I), xys, cder);

  // IAU 2006 precession + IAU2000AR06 nutation, using classical angles for dates outside 
  I = find(~t_interp);
  [x(I), y(I), s(I)] = CL__iers_xys_block(jd(:,I), "classic", cder);
  
  
  
// Compute X,Y,s and their time derivatives
else
  xdot = zeros(1,N);
  ydot = zeros(1,N);
  sdot = zeros(1,N);

  // Interpolation for dates inside pre-computed values
  I = find(t_interp);
  [x(I), y(I), s(I), xdot(I), ydot(I), sdot(I)] = CL__iers_xys_interp1(t(I), xys, xysdot, cder);

  // IAU 2006 precession + IAU2000AR06 nutation, using classical angles for dates outside 
  I = find(~t_interp);
  [x(I), y(I), s(I), xdot(I), ydot(I), sdot(I)] = CL__iers_xys_block(jd(:,I), "classic", cder);
end

endfunction