sofars 0.6.0

Pure Rust implementation of the IAU SOFA library
Documentation
use crate::vm::{anp, c2s, s2c, trxp, zp};

use super::{IauAstrom, IauLdBody, ab, ldn};

///  Quick CIRS −> ICRS, multiple deflections
///
///  Quick CIRS to ICRS astrometric place transformation, given the star-
///  independent astrometry parameters plus a list of light-deflecting
///  bodies.
///
///  Use of this function is appropriate when efficiency is important and
///  where many star positions are all to be transformed for one date.
///  The star-independent astrometry parameters can be obtained by
///  calling one of the functions iauApci[13], iauApcg[13], iauApco[13]
///  or iauApcs[13].
///
///  If the only light-deflecting body to be taken into account is the
///  Sun, the iauAticq function can be used instead.
///
///  This function is part of the International Astronomical Union's
///  SOFA (Standards of Fundamental Astronomy) software collection.
///
///  Status:  support function.
///
///  Given:
///  ```
///     ri,di  double      CIRS RA,Dec (radians)
///     astrom iauASTROM*  star-independent astrometry parameters:
///      pmt    double       PM time interval (SSB, Julian years)
///      eb     double[3]    SSB to observer (vector, au)
///      eh     double[3]    Sun to observer (unit vector)
///      em     double       distance from Sun to observer (au)
///      v      double[3]    barycentric observer velocity (vector, c)
///      bm1    double       sqrt(1-|v|^2): reciprocal of Lorenz factor
///      bpn    double[3][3] bias-precession-nutation matrix
///      along  double       longitude + s' (radians)
///      xpl    double       polar motion xp wrt local meridian (radians)
///      ypl    double       polar motion yp wrt local meridian (radians)
///      sphi   double       sine of geodetic latitude
///      cphi   double       cosine of geodetic latitude
///      diurab double       magnitude of diurnal aberration vector
///      eral   double       "local" Earth rotation angle (radians)
///      refa   double       refraction constant A (radians)
///      refb   double       refraction constant B (radians)
///     n      int          number of bodies (Note 3)
///     b      iauLDBODY[n] data for each of the n bodies (Notes 3,4):
///      bm     double       mass of the body (solar masses, Note 5)
///      dl     double       deflection limiter (Note 6)
///      pv     [2][3]       barycentric PV of the body (au, au/day)
///  ```
///  Returned:
///  ```
///     rc,dc  double     ICRS astrometric RA,Dec (radians)
///  ```
///  Notes:
///
///  1) Iterative techniques are used for the aberration and light
///     deflection corrections so that the functions iauAticqn and
///     iauAtciqn are accurate inverses; even at the edge of the Sun's
///     disk the discrepancy is only about 1 nanoarcsecond.
///
///  2) If the only light-deflecting body to be taken into account is the
///     Sun, the iauAticq function can be used instead.
///
///  3) The struct b contains n entries, one for each body to be
///     considered.  If n = 0, no gravitational light deflection will be
///     applied, not even for the Sun.
///
///  4) The struct b should include an entry for the Sun as well as for
///     any planet or other body to be taken into account.  The entries
///     should be in the order in which the light passes the body.
///
///  5) In the entry in the b struct for body i, the mass parameter
///     b[i].bm can, as required, be adjusted in order to allow for such
///     effects as quadrupole field.
///
///  6) The deflection limiter parameter b[i].dl is phi^2/2, where phi is
///     the angular separation (in radians) between star and body at
///     which limiting is applied.  As phi shrinks below the chosen
///     threshold, the deflection is artificially reduced, reaching zero
///     for phi = 0.   Example values suitable for a terrestrial
///     observer, together with masses, are as follows:
///  ```
///        body i     b[i].bm        b[i].dl
///
///        Sun        1.0            6e-6
///        Jupiter    0.00095435     3e-9
///        Saturn     0.00028574     3e-10
///  ```
///  7) For efficiency, validation of the contents of the b array is
///     omitted.  The supplied masses must be greater than zero, the
///     position and velocity vectors must be right, and the deflection
///     limiter greater than zero.
///
///  Called:
///  ```
///     iauS2c       spherical coordinates to unit vector
///     iauTrxp      product of transpose of r-matrix and p-vector
///     iauZp        zero p-vector
///     iauAb        stellar aberration
///     iauLdn       light deflection by n bodies
///     iauC2s       p-vector to spherical
///     iauAnp       normalize angle into range +/- pi
///  ```
pub fn aticqn(ri: f64, di: f64, astrom: &mut IauAstrom, n: i32, b: &[IauLdBody]) -> (f64, f64) {
    let pi: [f64; 3];
    let mut ppr = [0.0; 3];
    let mut pnat = [0.0; 3];
    let mut pco = [0.0; 3];
    let mut w;
    let mut d = [0.0; 3];
    let mut before = [0.0; 3];
    let mut r2;
    let mut r;
    let mut after: [f64; 3];

    // CIRS RA,Dec to Cartesian.
    pi = s2c(ri, di);

    // Bias-precession-nutation, giving GCRS proper direction.
    trxp(&astrom.bpn, &pi, &mut ppr);

    // Aberration, giving GCRS natural direction.
    for _ in 0..2 {
        r2 = 0.0;
        for i in 0..3 {
            w = ppr[i] - d[i];
            before[i] = w;
            r2 += w * w;
        }
        r = r2.sqrt();
        for i in 0..3 {
            before[i] /= r;
        }
        after = ab(&before, &astrom.v, astrom.em, astrom.bm1);
        r2 = 0.0;
        for i in 0..3 {
            d[i] = after[i] - before[i];
            w = ppr[i] - d[i];
            pnat[i] = w;
            r2 += w * w;
        }
        r = r2.sqrt();
        for i in 0..3 {
            pnat[i] /= r;
        }
    }

    // Light deflection, giving BCRS coordinate direction.
    zp(&mut d);
    for _ in 0..5 {
        r2 = 0.0;
        for i in 0..3 {
            w = pnat[i] - d[i];
            before[i] = w;
            r2 += w * w;
        }
        r = r2.sqrt();
        for i in 0..3 {
            before[i] /= r;
        }
        after = ldn(n, b, &astrom.eb, &before);
        r2 = 0.0;
        for i in 0..3 {
            d[i] = after[i] - before[i];
            w = pnat[i] - d[i];
            pco[i] = w;
            r2 += w * w;
        }
        r = r2.sqrt();
        for i in 0..3 {
            pco[i] /= r;
        }
    }

    // ICRS astrometric RA,Dec.
    let (w, dc) = c2s(&pco);

    (anp(w), dc)
}