sofars 0.6.0

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])
}