rfa 0.5.9

use crate::{rfam::*, utils::*};
use crate::fundamental_args::{fal03::*, faf03::*, falp03::*, fave03::*, fapa03::*, fad03::*, faom03::*, fae03::*, };
///  The CIO locator s, positioning the Celestial Intermediate Origin on
///  the equator of the Celestial Intermediate Pole, given the CIP's X,Y
///  coordinates.  Compatible with IAU 2000A precession-nutation.
///
///  Given:
///     date1,date2   TT as a 2-part Julian Date (Note 1)
///     x,y           CIP coordinates (Note 3)
///
///  Returned (function value):
///                   double    the CIO locator s in radians (Note 2)
///
/// # Notes:
///
///  1) The TT date date1+date2 is a Julian Date, apportioned in any
///     convenient way between the two arguments.  For example,
///     JD(TT)=2450123.7 could be expressed in any of these ways,
///     among others:
///
///     |    date1    |      date2   |                      |
///     |-------------|--------------|----------------------|
///     |2450123.7    |       0.0    |  (JD method)         |
///     |2451545.0    |   -1421.3    |  (J2000 method)      |
///     |2400000.5    |   50123.2    |  (MJD method)        |
///     |2450123.5    |       0.2    | (date & time method) |
///
///     The JD method is the most natural and convenient to use in
///     cases where the loss of several decimal digits of resolution
///     is acceptable.  The J2000 method is best matched to the way
///     the argument is handled internally and will deliver the
///     optimum resolution.  The MJD method and the date & time methods
///     are both good compromises between resolution and convenience.
///
///  2) The CIO locator s is the difference between the right ascensions
///     of the same point in two systems:  the two systems are the GCRS
///     and the CIP,CIO, and the point is the ascending node of the
///     CIP equator.  The quantity s remains below 0.1 arcsecond
///     throughout 1900-2100.
///
///  3) The series used to compute s is in fact for s+XY/2, where X and Y
///     are the x and y components of the CIP unit vector;  this series
///     is more compact than a direct series for s would be.  This
///     function requires X,Y to be supplied by the caller, who is
///     responsible for providing values that are consistent with the
///     supplied date.
///
///  4) The model is consistent with the IAU 2000A precession-nutation.
///
/// # Called:
///    * fal03     mean anomaly of the Moon
///    * falp03    mean anomaly of the Sun
///    * faf03     mean argument of the latitude of the Moon
///    * fad03     mean elongation of the Moon from the Sun
///    * faom03    mean longitude of the Moon's ascending node
///    * fave03    mean longitude of Venus
///    * fae03     mean longitude of Earth
///    * fapa03    general accumulated precession in longitude
///
/// # References:
///
///     Capitaine, N., Chapront, J., Lambert, S. and Wallace, P.,
///     "Expressions for the Celestial Intermediate Pole and Celestial
///     Ephemeris Origin consistent with the IAU 2000A precession-
///     nutation model", Astron.Astrophys. 400, 1145-1154 (2003)
///
///     n.b. The celestial ephemeris origin (CEO) was renamed "celestial
///          intermediate origin" (CIO) by IAU 2006 Resolution 2.
///
///     McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
///     IERS Technical Note No. 32, BKG (2004)
///
///  This revision:  2021 May 11
pub fn s00(date1: f64, date2: f64, x: f64, y: f64)->f64
{

/* Fundamental arguments */
   let mut fa = [0.0; 8];

/* --------------------- */
/* The series for s+XY/2 */
/* --------------------- */

   struct TERM {
      nfa: [i32; 8],      /* coefficients of l,l',F,D,Om,LVe,LE,pA */
      s: f64, c: f64,      /* sine and cosine coefficients */
   }
   impl TERM {
    pub const fn new (nfa: [i32; 8], s: f64, c: f64,) ->Self{
        TERM { nfa: nfa, s: s, c: c }
    }
       
   } 

/* Polynomial coefficients */
   const SP: [f64; 6] = [

   /* 1-6 */
          94.00e-6,
        3808.35e-6,
        -119.94e-6,
      -72574.09e-6,
          27.70e-6,
          15.61e-6
   ];

/* Terms of order t^0 */
   const S0: [TERM; 33] = [

   /* 1-10 */
TERM::new([ 0,  0,  0,  0,  1,  0,  0,  0], -2640.73e-6,   0.39e-6 ), 
TERM::new([ 0,  0,  0,  0,  2,  0,  0,  0],   -63.53e-6,   0.02e-6 ), 
TERM::new([ 0,  0,  2, -2,  3,  0,  0,  0],   -11.75e-6,  -0.01e-6 ), 
TERM::new([ 0,  0,  2, -2,  1,  0,  0,  0],   -11.21e-6,  -0.01e-6 ), 
TERM::new([ 0,  0,  2, -2,  2,  0,  0,  0],     4.57e-6,   0.00e-6 ), 
TERM::new([ 0,  0,  2,  0,  3,  0,  0,  0],    -2.02e-6,   0.00e-6 ), 
TERM::new([ 0,  0,  2,  0,  1,  0,  0,  0],    -1.98e-6,   0.00e-6 ), 
TERM::new([ 0,  0,  0,  0,  3,  0,  0,  0],     1.72e-6,   0.00e-6 ), 
TERM::new([ 0,  1,  0,  0,  1,  0,  0,  0],     1.41e-6,   0.01e-6 ), 
TERM::new([ 0,  1,  0,  0, -1,  0,  0,  0],     1.26e-6,   0.01e-6 ), 

   /* 11-20 */
TERM::new([ 1,  0,  0,  0, -1,  0,  0,  0],     0.63e-6,   0.00e-6 ), 
TERM::new([ 1,  0,  0,  0,  1,  0,  0,  0],     0.63e-6,   0.00e-6 ), 
TERM::new([ 0,  1,  2, -2,  3,  0,  0,  0],    -0.46e-6,   0.00e-6 ), 
TERM::new([ 0,  1,  2, -2,  1,  0,  0,  0],    -0.45e-6,   0.00e-6 ), 
TERM::new([ 0,  0,  4, -4,  4,  0,  0,  0],    -0.36e-6,   0.00e-6 ), 
TERM::new([ 0,  0,  1, -1,  1, -8, 12,  0],     0.24e-6,   0.12e-6 ), 
TERM::new([ 0,  0,  2,  0,  0,  0,  0,  0],    -0.32e-6,   0.00e-6 ), 
TERM::new([ 0,  0,  2,  0,  2,  0,  0,  0],    -0.28e-6,   0.00e-6 ), 
TERM::new([ 1,  0,  2,  0,  3,  0,  0,  0],    -0.27e-6,   0.00e-6 ), 
TERM::new([ 1,  0,  2,  0,  1,  0,  0,  0],    -0.26e-6,   0.00e-6 ), 

   /* 21-30 */
TERM::new([ 0,  0,  2, -2,  0,  0,  0,  0],     0.21e-6,   0.00e-6 ), 
TERM::new([ 0,  1, -2,  2, -3,  0,  0,  0],    -0.19e-6,   0.00e-6 ), 
TERM::new([ 0,  1, -2,  2, -1,  0,  0,  0],    -0.18e-6,   0.00e-6 ), 
TERM::new([ 0,  0,  0,  0,  0,  8,-13, -1],     0.10e-6,  -0.05e-6 ), 
TERM::new([ 0,  0,  0,  2,  0,  0,  0,  0],    -0.15e-6,   0.00e-6 ), 
TERM::new([ 2,  0, -2,  0, -1,  0,  0,  0],     0.14e-6,   0.00e-6 ), 
TERM::new([ 0,  1,  2, -2,  2,  0,  0,  0],     0.14e-6,   0.00e-6 ), 
TERM::new([ 1,  0,  0, -2,  1,  0,  0,  0],    -0.14e-6,   0.00e-6 ), 
TERM::new([ 1,  0,  0, -2, -1,  0,  0,  0],    -0.14e-6,   0.00e-6 ), 
TERM::new([ 0,  0,  4, -2,  4,  0,  0,  0],    -0.13e-6,   0.00e-6 ), 

   /* 31-33 */
TERM::new([ 0,  0,  2, -2,  4,  0,  0,  0],     0.11e-6,   0.00e-6 ), 
TERM::new([ 1,  0, -2,  0, -3,  0,  0,  0],    -0.11e-6,   0.00e-6 ), 
TERM::new([ 1,  0, -2,  0, -1,  0,  0,  0],    -0.11e-6,   0.00e-6 )
   ];

/* Terms of order t^1 */
   const S1: [TERM; 3] =[

   /* 1-3 */
TERM::new([ 0,  0,  0,  0,  2,  0,  0,  0],    -0.07e-6,   3.57e-6 ), 
TERM::new([ 0,  0,  0,  0,  1,  0,  0,  0],     1.71e-6,  -0.03e-6 ), 
TERM::new([ 0,  0,  2, -2,  3,  0,  0,  0],     0.00e-6,   0.48e-6 )
   ];

/* Terms of order t^2 */
   const S2: [TERM; 25] =[

   /* 1-10 */
TERM::new([ 0,  0,  0,  0,  1,  0,  0,  0],   743.53e-6,  -0.17e-6 ), 
TERM::new([ 0,  0,  2, -2,  2,  0,  0,  0],    56.91e-6,   0.06e-6 ), 
TERM::new([ 0,  0,  2,  0,  2,  0,  0,  0],     9.84e-6,  -0.01e-6 ), 
TERM::new([ 0,  0,  0,  0,  2,  0,  0,  0],    -8.85e-6,   0.01e-6 ), 
TERM::new([ 0,  1,  0,  0,  0,  0,  0,  0],    -6.38e-6,  -0.05e-6 ), 
TERM::new([ 1,  0,  0,  0,  0,  0,  0,  0],    -3.07e-6,   0.00e-6 ), 
TERM::new([ 0,  1,  2, -2,  2,  0,  0,  0],     2.23e-6,   0.00e-6 ), 
TERM::new([ 0,  0,  2,  0,  1,  0,  0,  0],     1.67e-6,   0.00e-6 ), 
TERM::new([ 1,  0,  2,  0,  2,  0,  0,  0],     1.30e-6,   0.00e-6 ), 
TERM::new([ 0,  1, -2,  2, -2,  0,  0,  0],     0.93e-6,   0.00e-6 ), 

   /* 11-20 */
TERM::new([ 1,  0,  0, -2,  0,  0,  0,  0],     0.68e-6,   0.00e-6 ), 
TERM::new([ 0,  0,  2, -2,  1,  0,  0,  0],    -0.55e-6,   0.00e-6 ), 
TERM::new([ 1,  0, -2,  0, -2,  0,  0,  0],     0.53e-6,   0.00e-6 ), 
TERM::new([ 0,  0,  0,  2,  0,  0,  0,  0],    -0.27e-6,   0.00e-6 ), 
TERM::new([ 1,  0,  0,  0,  1,  0,  0,  0],    -0.27e-6,   0.00e-6 ), 
TERM::new([ 1,  0, -2, -2, -2,  0,  0,  0],    -0.26e-6,   0.00e-6 ), 
TERM::new([ 1,  0,  0,  0, -1,  0,  0,  0],    -0.25e-6,   0.00e-6 ), 
TERM::new([ 1,  0,  2,  0,  1,  0,  0,  0],     0.22e-6,   0.00e-6 ), 
TERM::new([ 2,  0,  0, -2,  0,  0,  0,  0],    -0.21e-6,   0.00e-6 ), 
TERM::new([ 2,  0, -2,  0, -1,  0,  0,  0],     0.20e-6,   0.00e-6 ), 

   /* 21-25 */
TERM::new([ 0,  0,  2,  2,  2,  0,  0,  0],     0.17e-6,   0.00e-6 ), 
TERM::new([ 2,  0,  2,  0,  2,  0,  0,  0],     0.13e-6,   0.00e-6 ), 
TERM::new([ 2,  0,  0,  0,  0,  0,  0,  0],    -0.13e-6,   0.00e-6 ), 
TERM::new([ 1,  0,  2, -2,  2,  0,  0,  0],    -0.12e-6,   0.00e-6 ), 
TERM::new([ 0,  0,  2,  0,  0,  0,  0,  0],    -0.11e-6,   0.00e-6 )
   ];

/* Terms of order t^3 */
   const S3: [TERM; 4] =[

   /* 1-4 */
TERM::new([ 0,  0,  0,  0,  1,  0,  0,  0],     0.30e-6, -23.51e-6 ), 
TERM::new([ 0,  0,  2, -2,  2,  0,  0,  0],    -0.03e-6,  -1.39e-6 ), 
TERM::new([ 0,  0,  2,  0,  2,  0,  0,  0],    -0.01e-6,  -0.24e-6 ), 
TERM::new([ 0,  0,  0,  0,  2,  0,  0,  0],     0.00e-6,   0.22e-6 )
   ];

/* Terms of order t^4 */
   const S4: [TERM; 1] =[

   /* 1-1 */
TERM::new([ 0,  0,  0,  0,  1,  0,  0,  0],    -0.26e-6,  -0.01e-6 )
   ];

/* ------------------------------------------------------------------ */

/* Interval between fundamental epoch J2000.0 and current date (JC). */
   let t = ((date1 - URSA_DJ00) + date2) / URSA_DJC;

/* Fundamental Arguments (from IERS Conventions 2003) */

/* Mean anomaly of the Moon. */
   fa[0] = fal03(t);

/* Mean anomaly of the Sun. */
   fa[1] = falp03(t);

/* Mean longitude of the Moon minus that of the ascending node. */
   fa[2] = faf03(t);

/* Mean elongation of the Moon from the Sun. */
   fa[3] = fad03(t);

/* Mean longitude of the ascending node of the Moon. */
   fa[4] = faom03(t);

/* Mean longitude of Venus. */
   fa[5] = fave03(t);

/* Mean longitude of Earth. */
   fa[6] = fae03(t);

/* General precession in longitude. */
   fa[7] = fapa03(t);

/* Evaluate s. */
   let mut w0 = SP[0];
   let mut w1 = SP[1];
   let mut w2 = SP[2];
   let mut w3 = SP[3];
   let mut w4 = SP[4];
   let  w5 = SP[5];

   for s0_i in S0.iter() {
      let mut a = 0.0;
      for j in 0..8 {
          a += s0_i.nfa[j] as f64 * fa[j];
      }
      w0 += s0_i.s * sin(a) + s0_i.c * cos(a);
   }

   for s1_i in S1.iter() {
      let mut a = 0.0;
      for j in 0..8 {
          a += s1_i.nfa[j] as f64 * fa[j];
      }
      w1 += s1_i.s * sin(a) + s1_i.c * cos(a);
   }

   for s2_i in S2.iter() {
      let mut a = 0.0;
      for j in 0..8 {
          a += s2_i.nfa[j] as f64 * fa[j];
      }
      w2 += s2_i.s * sin(a) + s2_i.c * cos(a);
   }

   for s3_i in S3.iter() {
      let mut a = 0.0;
      for j in 0..8 {
          a += s3_i.nfa[j] as f64 * fa[j];
      }
      w3 += s3_i.s * sin(a) + s3_i.c * cos(a);
   }

   for s4_i in S4.iter() {
    let mut a = 0.0;
      for j in 0..8 {
          a += s4_i.nfa[j] as f64 * fa[j];
      }
      w4 += s4_i.s * sin(a) + s4_i.c * cos(a);
   }

       (w0 +
       (w1 +
       (w2 +
       (w3 +
       (w4 +
        w5 * t) * t) * t) * t) * t) * URSA_DAS2R - x*y/2.0

/* Finished. */

}