use crate::{

    consts::{DAU, DAYSEC, DC, DJY, DR2AS},

    vm::{anp, pdp, pm, pmp, pn, ppp, pv2s, sxp},

};

///  Space motion pv−vector to star catalog data

///

///  Convert star position+velocity vector to catalog coordinates.

///

///  This function is part of the International Astronomical Union's

///  SOFA (Standards of Fundamental Astronomy) software collection.

///

///  Status:  support function.

///

///  Given (Note 1):

///  ```

///     pv     double[2][3]   pv-vector (au, au/day)

///  ```

///  Returned (Note 2):

///  ```

///     ra     double         right ascension (radians)

///     dec    double         declination (radians)

///     pmr    double         RA proper motion (radians/year)

///     pmd    double         Dec proper motion (radians/year)

///     px     double         parallax (arcsec)

///     rv     double         radial velocity (km/s, positive = receding)

///  ```

///  Returned (function value):

///  ```

///            int            status:

///                              0 = OK

///                             -1 = superluminal speed (Note 5)

///                             -2 = null position vector

///  ```

///  Notes:

///

///  1) The specified pv-vector is the coordinate direction (and its rate

///     of change) for the date at which the light leaving the star

///     reached the solar-system barycenter.

///

///  2) The star data returned by this function are "observables" for an

///     imaginary observer at the solar-system barycenter.  Proper motion

///     and radial velocity are, strictly, in terms of barycentric

///     coordinate time, TCB.  For most practical applications, it is

///     permissible to neglect the distinction between TCB and ordinary

///     "proper" time on Earth (TT/TAI).  The result will, as a rule, be

///     limited by the intrinsic accuracy of the proper-motion and

///     radial-velocity data;  moreover, the supplied pv-vector is likely

///     to be merely an intermediate result (for example generated by the

///     function iauStarpv), so that a change of time unit will cancel

///     out overall.

///

///     In accordance with normal star-catalog conventions, the object's

///     right ascension and declination are freed from the effects of

///     secular aberration.  The frame, which is aligned to the catalog

///     equator and equinox, is Lorentzian and centered on the SSB.

///

///     Summarizing, the specified pv-vector is for most stars almost

///     identical to the result of applying the standard geometrical

///     "space motion" transformation to the catalog data.  The

///     differences, which are the subject of the Stumpff paper cited

///     below, are:

///

///     (i) In stars with significant radial velocity and proper motion,

///     the constantly changing light-time distorts the apparent proper

///     motion.  Note that this is a classical, not a relativistic,

///     effect.

///

///     (ii) The transformation complies with special relativity.

///

///  3) Care is needed with units.  The star coordinates are in radians

///     and the proper motions in radians per Julian year, but the

///     parallax is in arcseconds; the radial velocity is in km/s, but

///     the pv-vector result is in au and au/day.

///

///  4) The proper motions are the rate of change of the right ascension

///     and declination at the catalog epoch and are in radians per Julian

///     year.  The RA proper motion is in terms of coordinate angle, not

///     true angle, and will thus be numerically larger at high

///     declinations.

///

///  5) Straight-line motion at constant speed in the inertial frame is

///     assumed.  If the speed is greater than or equal to the speed of

///     light, the function aborts with an error status.

///

///  6) The inverse transformation is performed by the function iauStarpv.

///

///  Called:

///  ```

///     iauPn        decompose p-vector into modulus and direction

///     iauPdp       scalar product of two p-vectors

///     iauSxp       multiply p-vector by scalar

///     iauPmp       p-vector minus p-vector

///     iauPm        modulus of p-vector

///     iauPpp       p-vector plus p-vector

///     iauPv2s      pv-vector to spherical

///     iauAnp       normalize angle into range 0 to 2pi

///  ```

///  Reference:

///

///     Stumpff, P., 1985, Astron.Astrophys. 144, 232-240.

pub fn pvstar(pv: &[[f64; 3]; 2]) -> Result<[f64; 6], i32> {

    let mut pv = *pv;

    let d: f64;

    let w: f64;

    let del: f64;

    // Isolate the radial component of the velocity (au/day, inertial).

    let (_, pu) = pn(&pv[0]);

    let vr = pdp(&pu, &pv[1]);

    let ur = sxp(vr, &pu);

    // Isolate the transverse component of the velocity (au/day, inertial).

    let ut = pmp(&pv[1], &ur);

    let vt = pm(ut);

    // Special-relativity dimensionless parameters.

    let bett = vt / DC;

    let betr = vr / DC;

    // The observed-to-inertial correction terms.

    d = 1.0 + betr;

    w = betr * betr + bett * bett;

    if d == 0.0 || w > 1.0 {

        return Err(-1);

    }

    del = -w / ((1.0 - w).sqrt() + 1.0);

    // Scale inertial tangential velocity vector into observed (au/d).

    let ust = sxp(1.0 / d, &ut);

    // Compute observed radial velocity vector (au/d).

    let usr = sxp(DC * (betr - del) / d, &pu);

    // Combine the two to obtain the observed velocity vector.

    pv[1] = ppp(&usr, &ust);

    // Cartesian to spherical.

    let (mut ra, dec, r, rad, decd, rd) = pv2s(&pv);

    if r == 0.0 {

        return Err(-2);

    }

    // Return RA in range 0 to 2pi.

    ra = anp(ra);

    // Return proper motions in radians per year.

    let pmr = rad * DJY;

    let pmd = decd * DJY;

    // Return parallax in arcsec.

    let px = DR2AS / r;

    // Return radial velocity in km/s.

    let rv = 1e-3 * rd * DAU / DAYSEC;

    // Success.

    Ok([ra, dec, pmr, pmd, px, rv])

}