//  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 [tsid,tsiddot] = CL__mod_sidTimeVeis(jd, cder)
// Sidereal time: Gamma50 (Veis) frame to PEF
//
// Calling Sequence
// [tsid,tsiddot] = CL__mod_sidTimeVeis(jd [,cder])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the sidereal time relative to Gamma50 Veis reference frame (G50 Veis to PEF).</p> 
// <p>It also optionally computes the first derivative of the sidereal time : <b>tsiddot</b>, 
// which is a constant value. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// jd: Two-part julian day (Time scale: UT1) (2xN)
// cder: (boolean, optional) Option to compute derivative. If cder is %f, the derivative will be set to []. Default is %t. (1x1)
// tsid: Gamma50 Veis sidereal time [rad] (1xN)
// tsiddot: Gamma50 Veis sidereal time first derivative [rad/s] (1xN)
//
// Bibliography
// 1) Georges Veis - the system reference, in Sao special report 200, 1st volume, 1966
//
// See also
// CL_mod_sidTime
//
// Authors
// CNES - DCT/SB
//

// --------------------------------------------------------
//
// Reference: Georges Veis - the system reference,
// in Sao special report 200, 1st volume, 1966
// gives the values:  
// a1=100.075542 degrees
// a2=360.985612288 degrees per day
//
// (veis) Sideral time: 
// gmst = ( a1 + a2*t ) CL_rMod 2*pi
// 
// t: number of days since 1950.0 (UT1)
//
// Let's write: t = jj + frac,  with jj: integer
//
// gmst  = a1 + a2 * t                            [mod 2*pi]
//       = a1 + a2 * (jj + frac)                  [mod 2*pi]
//       = a1 + (a2-2*pi)* jj + a2 * frac         [mod 2*pi]
//       = a1 + (a2-2*pi)* (t-frac) + a2 * frac   [mod 2*pi]
//       = a1 + (a2-2*pi)* t + 2*pi * frac        [mod 2*pi]
//       = a1 + eps*t + 2*pi * frac               [mod 2*pi]
//         (with eps = a2 - 2*pi rad/day)
//
// --------------------------------------------------------


// Declarations:


// Code:

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

// a1 = 100.075542 degrees
a1 = 1.7466477086178711333726904629689443812586; // rad 

// a2 = 360.985612288 deg/day
a2 = 6.3003874867533010742986650751563511671780; // rad/day

// eps = a2 - 2*pi = 0.985612288 deg/day
eps = 0.172021795737145973733783085973453987844E-1; // rad 

// a3 = a2 / 86400; 
a3 = 0.729211514670520957673456605920873977E-4;  // rad/s

// Reference epoch (J1950.0), Julian Date = 2433282.5 UT1
DJ1950 = 2433282.0;  // NB: real value = DJ1950 + 0.5

// Days since fundamental epoch
// NB: jd(1,:) - DJ1950: supposed exact! 
t = (jd(1,:) - DJ1950) + jd(2,:) - 0.5; 

// Fractional part of t (days)
// NB: computation is modulo 1 day (see CL_rMod in tsid)
f = (jd(1,:) - floor(jd(1,:))) + (jd(2,:) - floor(jd(2,:))) - 0.5;

// Sidereal time at t
tsid = CL_rMod(a1 + 2*%pi * f + eps * t, 2*%pi);

// Time derivative
tsiddot = [];
if (cder)
  tsiddot = a3 * ones(t);  // rad/s
end

endfunction