//  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 [utc_tai] = CL__dat_tai2utc(cjd_tai)
// TAI to UTC conversion
//
// Calling Sequence
// utc_tai = CL__dat_tai2utc(cjd_tai)
//
// Description
// <itemizedlist><listitem>
// <p>Converts the time (from 1950.0) of an event given in the TAI time scale (<b>cjd_tai</b>) into 
// the time (from 1950.0) of the same event in the UTC time scale (cjd_utc = cjd_tai + utc_tai/86400).</p> 
// <p></p></listitem>
// <listitem>
// <p><b>Notes:</b></p>
// <p> - The conversion from TAI to UTC uses internal predefined data that define the leap seconds. </p>
// <p> - If a leap second is added (normally the case), the UTC time between 23h59mn59.999... UTC 
// (included) and 0h0mn0 UTC (excluded) cannot be represented; it 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_tai: Modified (1950.0) julian day in TAI time scale (1xN)
// utc_tai: UTC minus TAI = number of seconds to be added to cjd_tai to get cjd_utc (1xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_dat_utc2tai
//
// Examples
// cjd0 = CL_dat_utc2tai(21550);  // ref time in TAI
// cjd_tai = cjd0 + (-15 : 15) / 864000;  
// utc_tai = CL__dat_tai2utc(cjd_tai);
// [(cjd_tai-cjd0)*86400; (cjd_tai-cjd0)*86400 + utc_tai]'
//

// 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_tai2utc_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


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

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

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

  K = dat_tai2utc_search(cjd_tai, tjump_tai(1:$-1)); 
  utc_tai = -nsec(K); 

  // utc undefined if corresponding tai is in the leap second(s)
  dt_tai = tjump_tai(K+1) - cjd_tai; // > 0
  dn = dnsec(K+1); // > 0 if leap second added

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

endfunction


// Code:

// check size
Nr = size(cjd_tai,1); 
if (Nr >= 2); CL__error("Invalid size (number of rows)"); end
  
// default initialization
utc_tai = %nan * ones(cjd_tai); 

I = find(~isnan(cjd_tai)); 
if (I <> [])
  utc_tai(I) = dat_tai2utc_cv(cjd_tai(I)); 
end


endfunction