rkepler 0.4.2

/*
Copyright (c) 2022 Peter Hristozov / Петър Христозов

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/

//! Transformation between astronomical coordinate systems.

use crate::math::angle;
use crate::math::matrix::mat_r;
use crate::math::matrix::mat_r::MatR;
use crate::math::sphe::sph::Sph;
use crate::math::vector::vec_p::VecP;

/// L2,B2 system of galactic coordinates in the form presented in the
/// Hipparcos Catalogue.
///
/// In degrees:
/// * P = 192.85948    right ascension of the Galactic north pole in ICRS
/// - Q =  27.12825    declination of the Galactic north pole in ICRS
/// * R =  32.93192    Galactic longitude of the ascending node of
///                    the Galactic equator on the ICRS equator
///
/// ICRS to galactic rotation matrix, obtained by computing
/// R_3(-R) R_1(pi/2-Q) R_3(pi/2+P) to the full precision.
/// # Reference
/// Perryman M.A.C. & ESA, 1997, ESA SP-1200, The Hipparcos and Tycho
/// catalogues.  Astrometric and photometric star catalogues
/// derived from the ESA Hipparcos Space Astrometry Mission.  ESA
/// Publications Division, Noordwijk, Netherlands.
pub const R_ICRS_GAL: MatR = MatR {
    r00: -0.054875560416215368492398900454,
    r01: -0.873437090234885048760383168409,
    r02: -0.483835015548713226831774175116,
    r10: 0.494109427875583673525222371358,
    r11: -0.444829629960011178146614061616,
    r12: 0.746982244497218890527388004556,
    r20: -0.867666149019004701181616534570,
    r21: -0.198076373431201528180486091412,
    r22: 0.455983776175066922272100478348,
};

/// Equatorial to horizon coordinates:  transform hour angle and
/// declination to azimuth and altitude.
/// # Arguments
/// * `hd` local hour angle and declination in radians
/// * `phi` site latitude in radians
/// # Returns
/// Azimuth is returned in the range 0-2pi;  north is zero, and east
/// is +pi/2.  Altitude is returned in the range +/- pi/2.
pub fn eqv_hor(hd: &Sph, phi: f64) -> Sph {
    let (sh, ch) = hd.lon.sin_cos();
    let (sd, cd) = hd.lat.sin_cos();
    let (sp, cp) = phi.sin_cos();
    let x = -ch * cd * sp + sd * cp;
    let y = -sh * cd;
    let z = ch * cd * cp + sd * sp;
    let r = (x * x + y * y).sqrt();
    Sph {
        lon: (if r != 0.0 {
            angle::norm_dblpi(y.atan2(x))
        } else {
            0.0
        }),
        lat: (z.atan2(r)),
    }
}

/// Horizon to equatorial coordinates:  transform azimuth and altitude
/// to hour angle and declination.
/// # Arguments
/// * `ae` azimuth and latitude in radians
/// * `phi` site latitude in radians
/// # Returns
/// Local hour angle is returned in the range +/-pi.
/// Declination is returned in the range +/-pi/2.
pub fn hor_eqv(ae: &Sph, phi: f64) -> Sph {
    let (sa, ca) = ae.lon.sin_cos();
    let (se, ce) = ae.lat.sin_cos();
    let (sp, cp) = phi.sin_cos();
    let x = -ca * ce * sp + se * cp;
    let y = -sa * ce;
    let z = ca * ce * cp + se * sp;
    let r = (x * x + y * y).sqrt();
    Sph {
        lon: (if r != 0.0 { y.atan2(x) } else { 0.0 }),
        lat: (z.atan2(r)),
    }
}

/// Parallactic angle for a given hour angle and declination.
/// # Arguments
/// * `hd` hour angle and declination in radians
/// * `phi` site latitude in radians
/// # Returns
/// parallactic angle in the range -pi to +pi
/// # Reference
/// Smart, W.M., "Spherical Astronomy", Cambridge University Press,
/// 6th edition (Green, 1977), p49.
pub fn parallactic(hd: &Sph, phi: f64) -> f64 {
    let (sh, ch) = hd.lon.sin_cos();
    let (sd, cd) = hd.lat.sin_cos();
    let (sp, cp) = phi.sin_cos();
    let sqsz = cp * sh;
    let cqsz = sp * cd - cp * sd * ch;
    if sqsz != 0.0 || cqsz != 0.0 {
        sqsz.atan2(cqsz)
    } else {
        0.0
    }
}

/// Transformation from ICRS to Galactic Coordinates.
/// # Arguments
/// * `icrs` position vector in ICRS
/// # Returns
/// position vector in galactic coordinates
/// # References
/// Perryman M.A.C. & ESA, 1997, ESA SP-1200, The Hipparcos and Tycho
/// catalogues.  Astrometric and photometric star catalogues
/// derived from the ESA Hipparcos Space Astrometry Mission.  ESA
/// Publications Division, Noordwijk, Netherlands.
pub fn icrs_gal(icrs: &VecP) -> VecP {
    mat_r::mul_p(&R_ICRS_GAL, &icrs)
}

/// Transformation from Galactic Coordinates to ICRS.
/// # Arguments
/// * `gal` position vector in galactic coordinates
/// # Returns
/// position vector in ICRS
/// # References
/// Perryman M.A.C. & ESA, 1997, ESA SP-1200, The Hipparcos and Tycho
/// catalogues.  Astrometric and photometric star catalogues
/// derived from the ESA Hipparcos Space Astrometry Mission.  ESA
/// Publications Division, Noordwijk, Netherlands.
pub fn gal_icrs(gal: &VecP) -> VecP {
    mat_r::multr_p(&R_ICRS_GAL, &gal)
}

/// Equatorial to ecliptic rotation matrix.
/// # Arguments
/// * `obl` J2000, mean or true obliquity of ecliptic in radians
/// # Returns
/// equatorial to ecliptic rotation matrix
pub fn r_eqv_ecl(obl: f64) -> mat_r::MatR {
    mat_r::rot_x(obl)
}