//  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 [gmst82, gmst82dot] = CL__iers_gmst1982(jd_ut1, cder)
// Greenwich mean sidereal time (IAU 1982)
//
// Calling Sequence
// [gmst82, gmst82dot] = CL__iers_gmst1982(jd_ut1 [,cder])
//
// Description
// <itemizedlist><listitem>
// <p>Compute the Greenwich mean sidereal time (IAU 1982 model)</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// jd_ut1: 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)
// gmst82: Greenwich mean sidereal time (IAU 1982 model) [rad] (1xN)
// gmst82dot: (optional) Time derivative of Greenwich mean sidereal time (IAU 1982 model) [rad/s] (1xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) Technical Note 21, IERS, 1992
// 2) The new definition of universal time, Aoki et al, 1980


// Declarations:


// Code:

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

// conversion from seconds (angle) to radians
// (86400 sec <=> 2*%pi rad) 
SEC2RAD = 2*%pi / 86400;

// Equation for GMST at 0h UT1 (ref2, eq 13, p360): 
// GMST at 0h UT1= 6h41'50.54841" + 8640184.812866" * T + 0.093104" * T^2 - 6.2" x 10^-6 * T^3
// with T = t/36525, t being the number of days elapsed since 2000 January 1, 12h UT1. 
// t takes values:  ... -0.5, 0.5, 1.5, ...
//
c0 = 24110.54841 * SEC2RAD;   // 6h41'50.54841" = 24110.54841" 
c1 = 8640184.812866 * SEC2RAD;
c2 = 0.093104 * SEC2RAD;
c3 = -6.2e-6 * SEC2RAD;

// Number of Julian days since J2000 (J2000: at 12h) 
// NB: jd(1,:) minus 2451545 should be exact for maximum precision
t = ((jd_ut1(1,:) - 2451545) + jd_ut1(2,:));

// delta = difference between t and t0h in radians (1 day <=> 2*%pi rad). 
// t0h = Nearest day/time at 0h UT1, i.e. where the formula applies. 
// NB formula has to be exact mod 2%pi as there is a "modulo" in the formula of gmst82. 
// Formula used equivalent (mod 2*%pi) to: 
// t0h = round(t - 0.5) + 0.5
// delta = 2*%pi * (t - t0h); 
delta = 2*%pi * ((jd_ut1(1,:) - round(jd_ut1(1,:))) + (jd_ut1(2,:) - round(jd_ut1(2,:))) + 0.5); 
 
// Number of Julian centuries from J2000
tc = t / 36525; 

// GMST at any day/time : 
// The formula is initially given for each day at 0h UT1
// if GMST(t1) is the value at t1 (0h), the "real" (continuous) value at t2 (0h) is GMST(t2) + 2*%pi : 
// There is one revolution between t1 and t2. 
// By extension: the "real" (continuous) value at t (t between t1 (0h) and t2 (0h)) is: 
// GMST(t) + 2*%pi * (t-t1)
gmst82 = CL_rMod(c0 + (c1 + (c2 + c3 .* tc) .* tc) .* tc + delta, 2*%pi);

// Time derivative
// gmst82dot = d()/dt = d()/d(tc) * d(tc)/dt, 
// and: d(delta)/dt = 2*%pi / 86400 rad/s
gmst82dot = [];
if (cder)
  gmst82dot = ((c1 + (2*c2 + 3*c3 .* tc) .* tc) / 36525 + 2*%pi) / 86400;
end

endfunction