//  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_block(jd, model, cder)
// Utility function that computes X,Y,s using loops to limit memory consumption.
//
// Calling Sequence
// [x, y, s, xdot, ydot, sdot] = CL__iers_xys_block(jd, model [,cder])
//
// Description
// <itemizedlist><listitem>
// Utility function to compute X,Y,s based on the IAU 2006 precession and IAU 2000AR06 nutation models. 
// - If model = "classic", the model using classical angles is used. (See function CL__iers_xys2006A_cla)
// - If model = "series", the model using series is used. (See function CL__iers_xys2006A_ser)
//
// The function computes the results for small sub-blocks of dates in order to limit memory usage
// NB: uses stacksize. 
// </itemizedlist>
//
// Parameters
// jd: Two-part julian day (Time scale: TT) (2xN)
// model: (string) Model to be used. ("classic" or "series") (1x1)
// 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: Time derivative of X coordinate of CIP [rad/s] (1xN)
// ydot: Time derivative of Y coordinate of CIP [rad/s] (1xN)
// sdot: Time derivative of CIO locator [rad/s] (1xN)
//
// Authors
// CNES - DCT/SB


// Declarations:

// Code:

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

// Function to be called :
if (model == "classic")
  f = CL__iers_xys2006A_cla;
elseif (model == "series")
  f = CL__iers_xys2006A_ser;
else
  CL__error("Unknown model");
end 

if (CL__scilabVersion() < 600)
  // Maximum number of dates computable with current stacksize
  sz = stacksize() ; 
  Nmax = (sz(1) - sz(2)) / 10000;  
else
  // Block size is constant for Scilab version >= 6
  Nmax = 2000; 
end

N = size(jd,2);

nb_blocks = ceil(N/Nmax);

// I: indices of beginning/end of blocks:
// block 1: I(1) -> I(2)
// block 2: I(2) -> I(3) -- note: index I(2) computed twice
// block ("nb_block"): I(nb_block) -> I(nb_block+1)
I = round(linspace(1,N,nb_blocks+1));

x = zeros(1,N);
y = zeros(1,N);
s = zeros(1,N);

xdot = [];
ydot = [];
sdot = [];

// Compute X,Y,S and time derivatives
if (cder)
  xdot = zeros(1,N);
  ydot = zeros(1,N);
  sdot = zeros(1,N);
  for k = 1 : nb_blocks
    Ib = [I(k):I(k+1)];
    [x(Ib), y(Ib), s(Ib), xdot(Ib), ydot(Ib), sdot(Ib)] = f(jd(:,Ib),cder);
  end
 
// Compute X,Y,S (without time derivatives) 
else
  for k = 1 : nb_blocks
    Ib = [I(k):I(k+1)];
    [x(Ib), y(Ib), s(Ib)] = f(jd(:,Ib),cder);
  end
end

endfunction