#![allow(non_snake_case)]
#![allow(clippy::excessive_precision)]
use erfa_sys::{
eraAnp, eraC2s, eraEpj, eraEpj2jd, eraEpv00, eraGmst06, eraIr, eraP06e, eraPdp, eraPmat06,
eraPn, eraPnm06a, eraRx, eraRxp, eraRxpv, eraRxr, eraRz, eraS2c, ERFA_CMPS, ERFA_DAU,
ERFA_DAYSEC, ERFA_DJM0,
};
const ERFA_AULT: f64 = ERFA_DAU / ERFA_CMPS;
pub fn palGmst(ut1: f64) -> f64 {
unsafe { eraGmst06(ERFA_DJM0, ut1, ERFA_DJM0, ut1) }
}
pub unsafe fn palDcs2c(a: f64, b: f64, v: *mut f64) {
eraS2c(a, b, v)
}
pub unsafe fn palDcc2s(v: *mut f64, a: &mut f64, b: &mut f64) {
eraC2s(v, a, b)
}
pub fn palDranrm(angle: f64) -> f64 {
unsafe { eraAnp(angle) }
}
pub unsafe fn palDvn(v: *mut f64, uv: *mut f64, vm: *mut f64) {
eraPn(v, vm, uv)
}
pub unsafe fn palDvdv(va: *mut f64, vb: *mut f64) -> f64 {
eraPdp(va, vb)
}
pub unsafe fn palDmxv(dm: *mut [f64; 3], va: *mut f64, vb: *mut f64) {
eraRxp(dm, va, vb)
}
pub unsafe fn palEvp(
date: f64,
deqx: f64,
dvb: *mut f64,
dpb: *mut f64,
dvh: *mut f64,
dph: *mut f64,
) {
let mut pvh = [[0.0; 3]; 2];
let mut pvb = [[0.0; 3]; 2];
let mut d1 = 0.0;
let mut d2 = 0.0;
let mut r = [[0.0; 3]; 3];
eraEpv00(2400000.5, date, pvh.as_mut_ptr(), pvb.as_mut_ptr());
if deqx > 0.0 {
eraEpj2jd(deqx, &mut d1, &mut d2);
eraPmat06(d1, d2, r.as_mut_ptr());
eraRxpv(r.as_mut_ptr(), pvh.as_mut_ptr(), pvh.as_mut_ptr());
eraRxpv(r.as_mut_ptr(), pvb.as_mut_ptr(), pvb.as_mut_ptr());
}
let dvb = std::slice::from_raw_parts_mut(dvb, 3);
let dpb = std::slice::from_raw_parts_mut(dpb, 3);
let dvh = std::slice::from_raw_parts_mut(dvh, 3);
let dph = std::slice::from_raw_parts_mut(dph, 3);
for i in 0..3 {
dvh[i] = pvh[1][i] / ERFA_DAYSEC;
dvb[i] = pvb[1][i] / ERFA_DAYSEC;
dph[i] = pvh[0][i];
dpb[i] = pvb[0][i];
}
}
pub unsafe fn palPrenut(epoch: f64, date: f64, rmatpn: *mut [f64; 3]) {
let mut bpa = 0.0;
let mut bpia = 0.0;
let mut bqa = 0.0;
let mut chia = 0.0;
let mut d1 = 0.0;
let mut d2 = 0.0;
let mut eps0 = 0.0;
let mut epsa = 0.0;
let mut gam = 0.0;
let mut oma = 0.0;
let mut pa = 0.0;
let mut phi = 0.0;
let mut pia = 0.0;
let mut psi = 0.0;
let mut psia = 0.0;
let mut r1 = [[0.0; 3]; 3];
let mut r2 = [[0.0; 3]; 3];
let mut thetaa = 0.0;
let mut za = 0.0;
let mut zetaa = 0.0;
eraEpj2jd(epoch, &mut d1, &mut d2);
eraP06e(
d1,
d2,
&mut eps0,
&mut psia,
&mut oma,
&mut bpa,
&mut bqa,
&mut pia,
&mut bpia,
&mut epsa,
&mut chia,
&mut za,
&mut zetaa,
&mut thetaa,
&mut pa,
&mut gam,
&mut phi,
&mut psi,
);
eraIr(r1.as_mut_ptr());
eraRz(-chia, r1.as_mut_ptr());
eraRx(oma, r1.as_mut_ptr());
eraRz(psia, r1.as_mut_ptr());
eraRx(-eps0, r1.as_mut_ptr());
eraPnm06a(ERFA_DJM0, date, r2.as_mut_ptr());
eraRxr(r2.as_mut_ptr(), r1.as_mut_ptr(), rmatpn);
}
pub unsafe fn palMappa(eq: f64, date: f64, amprms: *mut f64) {
let gr2: f64 = 2.0 * 9.87063e-9;
let mut ebd = [0.0; 3];
let mut ehd = [0.0; 3];
let mut eh = [0.0; 3];
let mut e = 0.0;
let mut vn = [0.0; 3];
let mut vm = 0.0;
let amprms = std::slice::from_raw_parts_mut(amprms, 21);
amprms.fill(0.0);
amprms[0] = eraEpj(ERFA_DJM0, date) - eq;
palEvp(
date,
eq,
ebd.as_mut_ptr(),
&mut amprms[1],
ehd.as_mut_ptr(),
eh.as_mut_ptr(),
);
eraPn(eh.as_mut_ptr(), &mut e, &mut amprms[4]);
amprms[7] = gr2 / e;
for i in 0..3 {
amprms[i + 8] = ebd[i] * ERFA_AULT;
}
eraPn(&mut amprms[8], &mut vm, vn.as_mut_ptr());
amprms[11] = (1.0 - vm * vm).sqrt();
palPrenut(eq, date, amprms.as_mut_ptr().add(12) as _);
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_abs_diff_eq;
use ndarray::prelude::*;
#[test]
fn gmst() {
assert_abs_diff_eq!(palGmst(53736.), 1.754174971870091203, epsilon = 1e-12);
}
#[test]
fn vecmat() {
let mut dv1 = [0_f64; 3];
let mut dv2 = [0_f64; 3];
let mut dv3 = [0_f64; 3];
let mut dv5 = [0.006889040510209034, -1.577473205461961, 0.5201843672856759];
let mut dv6 = [0_f64; 3];
let mut dvm = 0_f64;
let mut drm1: [[f64; 3]; 3] = [
[0.9930075842721269, 0.05902743090199868, -0.1022335560329612],
[
-0.07113807138648245,
0.9903204657727545,
-0.1191836812279541,
],
[0.09420887631983825, 0.1256229973879967, 0.9875948309655174],
];
let mut drm2: [[f64; 3]; 3] = [
[-0.1681574770810878, 0.1981362273264315, 0.9656423242187410],
[-0.2285369373983370, 0.9450659587140423, -0.2337117924378156],
[
-0.9589024617479674,
-0.2599853247796050,
-0.1136384607117296,
],
];
unsafe { palDcs2c(3.0123, -0.999, dv1.as_mut_ptr()) };
assert_abs_diff_eq!(dv1[0], -0.5366267667260525, epsilon = 1e-12);
assert_abs_diff_eq!(dv1[1], 0.06977111097651444, epsilon = 1e-12);
assert_abs_diff_eq!(dv1[2], -0.8409302618566215, epsilon = 1e-12);
unsafe { palDmxv(drm1.as_mut_ptr(), dv1.as_mut_ptr(), dv2.as_mut_ptr()) };
unsafe { palDmxv(drm2.as_mut_ptr(), dv2.as_mut_ptr(), dv3.as_mut_ptr()) };
assert_abs_diff_eq!(dv3[0], -0.7267487768696160, epsilon = 1e-10);
assert_abs_diff_eq!(dv3[1], 0.5011537352639822, epsilon = 1e-12);
assert_abs_diff_eq!(dv3[2], 0.4697671220397141, epsilon = 1e-12);
dv5.iter_mut().for_each(|x| *x *= 1000.0);
unsafe { palDvn(dv5.as_mut_ptr(), dv6.as_mut_ptr(), &mut dvm) };
assert_abs_diff_eq!(dv6[0], 0.004147420704640065, epsilon = 1e-12);
assert_abs_diff_eq!(dv6[1], -0.9496888606842218, epsilon = 1e-12);
assert_abs_diff_eq!(dv6[2], 0.3131674740355448, epsilon = 1e-12);
assert_abs_diff_eq!(dvm, 1661.042127339937, epsilon = 1e-9);
}
#[test]
fn dcc2s() {
let mut dv = [100_f64, -50_f64, 25_f64];
let mut da = 0_f64;
let mut db = 0_f64;
unsafe { palDcc2s(dv.as_mut_ptr(), &mut da, &mut db) };
assert_abs_diff_eq!(da, -0.4636476090008061, epsilon = 1e-12);
assert_abs_diff_eq!(db, 0.2199879773954594, epsilon = 1e-12);
}
#[test]
fn dranrm() {
assert_abs_diff_eq!(palDranrm(-0.1), 6.183185307179587, epsilon = 1e-12);
}
#[test]
fn evp() {
let mut dvb = Array1::zeros(3);
let mut dpb = Array1::zeros(3);
let mut dvh = Array1::zeros(3);
let mut dph = Array1::zeros(3);
let vbex = array![
1.6957348127008098514e-07,
-9.1093446116039685966e-08,
-3.9528532243991863036e-08,
];
let pbex = array![
-0.49771075259730546136,
-0.80273812396332311359,
-0.34851593942866060383,
];
let vhex = array![
1.6964379181455713805e-07,
-9.1147224045727438391e-08,
-3.9553158272334222497e-08,
];
let phex = array![
-0.50169124421419830639,
-0.80650980174901798492,
-0.34997162028527262212,
];
unsafe {
palEvp(
2010.0,
2012.0,
dvb.as_mut_ptr(),
dpb.as_mut_ptr(),
dvh.as_mut_ptr(),
dph.as_mut_ptr(),
)
};
assert_abs_diff_eq!(dvb, vbex, epsilon = 1e-12);
assert_abs_diff_eq!(dpb, pbex, epsilon = 1e-12);
assert_abs_diff_eq!(dvh, vhex, epsilon = 1e-12);
assert_abs_diff_eq!(dph, phex, epsilon = 1e-12);
}
#[test]
fn prenut() {
}
#[test]
fn mappa() {
let mut amprms = Array1::zeros(21);
let expected = array![
1.9986310746064646082_f64,
-0.1728200754134739392,
0.88745394651412767839,
0.38472374350184274094,
-0.17245634725219796679,
0.90374808622520386159,
0.3917884696321610738,
2.0075929387510784968e-08,
-9.9464149073251757597e-05,
-1.6125306981057062306e-05,
-6.9897255793245634435e-06,
0.99999999489900059935,
0.99999983777998024959,
-0.00052248206600935195865,
-0.00022683144398381763045,
0.00052248547063364874764,
0.99999986339269864022,
1.4950491424992534218e-05,
0.00022682360163333854623,
-1.5069005133483779417e-05,
0.99999997416198904698,
];
unsafe { palMappa(2010.0, 55927.0, amprms.as_mut_ptr()) };
assert_abs_diff_eq!(amprms, expected, epsilon = 1e-12);
}
}