extern crate rkepler;
use rkepler::BodyId;
use rkepler::datetime::mjd;
use rkepler::datetime::mjd::Adate;
use rkepler::datetime::mjd::SECDAY;
use rkepler::datetime::time;
use rkepler::frame::constant::CLIGHT;
use rkepler::frame::iau06::iau06_r;
use rkepler::math::angle;
use rkepler::math::angle::{RADASEC, RADDEG};
use rkepler::math::matrix::mat_r;
use rkepler::math::sphe::sph_p;
use rkepler::math::vector::vec_p;
use rkepler::phys_eph::body_prm;
use rkepler::phys_eph::constell_name;
use rkepler::phys_eph::phys_eph;
use rkepler::phys_eph::plngrph;
use rkepler::pos_eph::satell::elp85;
use rkepler::pos_eph::vsop2013::vsop2013_lo;
struct Prm {
mjd_zt: f64,
earth_moon: vec_p::VecP,
sun_moon: vec_p::VecP,
icrf_true: mat_r::MatR,
rot: plngrph::PlnGrph,
}
struct Eph {
mjd_zt: f64,
constell: String,
ra: f64,
dec: f64,
prlx: f64,
phase: f64,
diam: f64,
lon_se: f64,
lat_se: f64,
pa_se: f64,
}
static EMRAT: f64 = 81.30057;
fn main() {
let start_date = Adate {
year: 2025,
month: 7,
day: 1,
hour: 0,
min: 0,
sec: 0,
msec: 0,
};
let end_date = Adate {
year: 2025,
month: 8,
day: 1,
hour: 0,
min: 0,
sec: 0,
msec: 0,
};
let step = 1.0;
let dzt = 3.0 / 24.0;
let dt = 37.0 / SECDAY;
let is_cirs = false;
let is_prnt_pr = false;
let is_dms = true;
let start_mjd = start_date.mjd();
let end_mjd = end_date.mjd();
let mut step_mjd = start_mjd;
let mut prm: Vec<Prm> = Vec::new();
let mut eph: Vec<Eph> = Vec::new();
let phys_moon = body_prm::moon();
let phys_earth = body_prm::earth();
while step_mjd < end_mjd {
let step_zt = step_mjd;
let step_tt = step_zt - dzt + dt + time::tt_m_tai();
let step_tdb = step_tt + time::tdb_m_tt_lo(step_tt);
let mut earth_moon = elp85::get_icrs_p(step_tdb, &BodyId::MOON);
let ligh_time = earth_moon.modul() / CLIGHT;
earth_moon = elp85::get_icrs_p(step_tdb - ligh_time, &BodyId::MOON);
let emb = vsop2013_lo::get_icrs_p(step_tdb - ligh_time, &BodyId::EMB);
let sun_moon = vec_p::add_scaled(&emb, 1.0 - 1.0 / (1.0 + EMRAT), &earth_moon);
let nut = iau06_r::nut00b(step_tt);
let icrf_true = tr_equ(step_tt, nut, is_cirs);
let rot_prm = plngrph::moon(step_tt - ligh_time);
prm.push(Prm {
mjd_zt: step_zt,
earth_moon,
sun_moon,
icrf_true,
rot: rot_prm,
});
step_mjd += step;
}
for i in prm.iter() {
let moon_true = mat_r::mul_p(&i.icrf_true, &i.earth_moon);
let moon_icrs_sph = sph_p::from_cart(&i.earth_moon);
let moon_true_sph = sph_p::from_cart(&moon_true);
let constell = constell_name::constell_name(&i.earth_moon);
let eqv_prl = (phys_earth.rad_equ / moon_icrs_sph.rad).asin();
let ph_ang = phys_eph::phase_angl(&i.sun_moon, &i.earth_moon);
let phase = phys_eph::phase(ph_ang);
let diam = 2.0 * phys_moon.vis_equ_rad(moon_icrs_sph.rad);
let (lon_sep, lat_sep, pa) = i.rot.sub_obs(&moon_true, &i.icrf_true);
eph.push(Eph {
mjd_zt: i.mjd_zt,
constell,
ra: moon_true_sph.lon,
dec: moon_true_sph.lat,
prlx: eqv_prl,
phase,
diam,
lon_se: angle::norm_pi(lon_sep),
lat_se: lat_sep,
pa_se: pa,
});
}
println!();
println!("Body: {}", "Moon");
for i in eph.iter() {
let dat = mjd::mjd_adate(i.mjd_zt);
print!(
" {:#04}-{:#02}-{:#02}T{:#02}:{:#02}:{:#02}.{:#03}",
dat.year, dat.month, dat.day, dat.hour, dat.min, dat.sec, dat.msec
);
if is_prnt_pr {
if is_dms {
let ra_hms = angle::rad_hms(i.ra, 3);
let dec_dms = angle::rad_dms(i.dec, 2);
let prlx_ms = angle::rad_dms(i.prlx, 2);
let diam_ms = angle::rad_dms(i.diam, 2);
println!(
" {} {:#2}h{:#02}m{:#02}.{:#03}s {}{:#2}°{:#02}ʹ{:#02}.{:#02}ʺ {:#02}ʹ{:#02}.{:#02}ʺ {:5.3} {:#02}ʹ{:#02}.{:#02}ʺ {:+5.2} {:+5.2} {:6.2}",
i.constell,
ra_hms.1,
ra_hms.2,
ra_hms.3,
ra_hms.4,
dec_dms.0,
dec_dms.1,
dec_dms.2,
dec_dms.3,
dec_dms.4,
prlx_ms.2,
prlx_ms.3,
prlx_ms.4,
i.phase,
diam_ms.2,
diam_ms.3,
diam_ms.4,
i.lon_se * RADDEG,
i.lat_se * RADDEG,
angle::norm_dblpi(i.pa_se) * RADDEG,
);
} else {
println!(
" {} {:10.7} {:+10.6} {:6.2} {:6.4} {:6.2} {:6.2} {:+6.2} {:6.2}",
i.constell,
i.ra * RADDEG / 15.0,
i.dec * RADDEG,
i.prlx * RADASEC,
i.phase,
i.diam * RADASEC,
i.lon_se * RADDEG,
i.lat_se * RADDEG,
angle::norm_dblpi(i.pa_se) * RADDEG,
);
}} else {
if is_dms {
let ra_hms = angle::rad_hms(i.ra, 0);
let dec_dms = angle::rad_dms(i.dec, -1);
let prlx_ms = angle::rad_dms(i.prlx, 0);
let diam_ms = angle::rad_dms(i.diam, 0);
println!(
" {} {:#2}h{:#02}m{:#02}s {}{:#2}°{:#02}.{}ʹ {:#02}ʹ{:#02}ʺ {:4.2} {:#02}ʹ{:#02}ʺ {:+4.1} {:+4.1} {:5.1}",
i.constell,
ra_hms.1,
ra_hms.2,
ra_hms.3,
dec_dms.0,
dec_dms.1,
dec_dms.2,
dec_dms.3 / 6,
prlx_ms.2,
prlx_ms.3,
i.phase,
diam_ms.2,
diam_ms.3,
i.lon_se * RADDEG,
i.lat_se * RADDEG,
angle::norm_dblpi(i.pa_se) * RADDEG,
);
} else {
println!(
" {} {:7.4} {:+7.3} {:4.0} {:4.2} {:4.0} {:+4.1} {:+4.1} {:5.1}",
i.constell,
i.ra * RADDEG / 15.0,
i.dec * RADDEG,
i.prlx * RADASEC,
i.phase,
i.diam * RADASEC,
i.lon_se * RADDEG,
i.lat_se * RADDEG,
angle::norm_dblpi(i.pa_se) * RADDEG,
);
}}
}
}
fn tr_equ(mjd_tt: f64, nut: (f64, f64), is_cirs: bool) -> mat_r::MatR {
let bpn = iau06_r::gcrs_true(mjd_tt, nut.0, nut.1);
if is_cirs {
let xy = iau06_r::bpn_xy(&bpn);
let s = iau06_r::s_lo(mjd_tt, xy.0, xy.1);
iau06_r::gcrs_cirs(xy.0, xy.1, s)
} else {
bpn
}
}