//  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 [tai_utc] = CL__dat_utc2tai(cjd_utc)
// UTC to TAI conversion
//
// Calling Sequence
// tai_utc = CL__dat_utc2tai(cjd_utc)
//
// Description
// <itemizedlist><listitem>
// <p>Converts the time (from 1950.0) of an event given in the UTC time scale (<b>cjd_utc</b>) into 
// the time (from 1950.0) of the same event in the TAI time scale (cjd_tai = cjd_utc + tai_utc/86400).</p> 
// <p></p></listitem>
// <listitem>
// <p><b>Notes:</b></p>
// <p> - The conversion from UTC to TAI uses internal predefined data that define the leap seconds. </p>
// <p> - If a leap second is removed (not normally the case), the UTC time between 
// 23h59mn58.999... UTC (included) and 0h0mn0 UTC (excluded) does not exist; 
// TAI time is then %nan. </p>
// <p></p></listitem>
// <listitem>
// <p>See <link linkend="Dates and time scales">Dates and time scales</link> for more details.</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// cjd_utc: Modified (1950.0) julian day in UTC time scale (1xN)
// tai_utc: TAI minus UTC = number of seconds to be added to cjd_utc to get cjd_tai (1xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_dat_tai2utc
//
// Examples
// cjd0 = 21550; // ref time in UTC 
// cjd_utc = cjd0 + (-15 : 15) / 864000;  
// tai_utc = CL__dat_utc2tai(cjd_utc);
// [(cjd_utc-cjd0)*86400; (cjd_utc-cjd0)*86400 + tai_utc]'


// Declarations:
TAI_UTC_DATA = CL__dataGetEnv("TAI_UTC_DATA", internal=%t);


// --------------------------------------
// returns k: tjmp(k) <= t < tjmp(k+1) 
// error if t < tjmp(1) or t >= tjmp($)
// k: same size as t
// --------------------------------------
function [k] = dat_utc2tai_search(t, tjmp)
  // NB : dsearch returns interval k if
  //      cjd_utc is in ] tjump_utc(k), tjump_utc(k+1) ]
  k = dsearch(t, tjmp, d='ch'); 

  if (find(k == 0 | k == length(tjmp)) <> [])
    CL__error("Anomaly in date conversion"); 
  end

  // looks for interval "[...[" instead of "]...]"
  I = find(t == tjmp(k+1)); 
  if (I <> []); 
    k(I) = k(I) + 1; 
  end
endfunction


// --------------------------------------
// UTC to TAI conversion
// NB: 
// - cjd_utc should not be %nan 
// - Uses global variable "TAI_UTC_DATA" 
// --------------------------------------
function [tai_utc] = dat_utc2tai_cv(cjd_utc)
  // Algorithm : 
  // We look for k such that : 
  //    t_utc in [ tjump_utc(k), tjump_utc(k+1) [ 
  // Then:
  //    t_tai = t_utc + nsec(k)  

  // Note: one element added (%nan) so that: 
  // tjump_utc(K+1) is always defined (see below)
  tjump_utc = [TAI_UTC_DATA(:,1); TAI_UTC_DATA($,2); %nan]'; 
  nsec = [TAI_UTC_DATA(:,3); %nan; %nan]';  

  // number of leap seconds added
  dnsec = [nsec(1), nsec(2:$) - nsec(1:$-1)];  

  K = dat_utc2tai_search(cjd_utc, tjump_utc(1:$-1)); 
  tai_utc = nsec(K); 

  // tai undefined if corresponding utc does not exist
  dt_utc = tjump_utc(K+1) - cjd_utc; // > 0
  dn = -dnsec(K+1); // > 0 if leap second removed

  // NB: I <> [] only possible if dn > 0
  I = find(dt_utc*86400 <= dn); 
  tai_utc(I) = %nan; 

endfunction


// Code:

// check size
Nr = size(cjd_utc,1); 
if (Nr >= 2); CL__error("Invalid size (number of rows)"); end

// default initialization
tai_utc = %nan * ones(cjd_utc); 

I = find(~isnan(cjd_utc)); 
if (I <> [])
  tai_utc(I) = dat_utc2tai_cv(cjd_utc(I)); 
end


endfunction