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::datetime::time::tt_m_tai;
use rkepler::frame::constant::CLIGHT;
use rkepler::frame::iau06::iau06_r;
use rkepler::math::angle;
use rkepler::math::angle::RADASEC;
use rkepler::math::angle::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::plngrph;
use rkepler::phys_eph::plngrph::PlnGrph;
use rkepler::pos_eph::satell::elp85;
use rkepler::pos_eph::vsop2013::vsop2013_hi;
use rkepler::pos_eph::vsop2013::vsop2013_lo;
struct Prm {
mjd_zt: f64,
sun_icrs: vec_p::VecP,
icrs_true: mat_r::MatR,
icrs_ecl: mat_r::MatR,
rot: PlnGrph,
}
struct Eph {
mjd_zt: f64,
lon: f64,
lat: f64,
rec: f64,
dec: f64,
dist: f64,
diam: f64,
constell: String,
lon_se: f64,
lat_se: f64,
pa_se: f64,
}
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_alg_hi = true;
let is_cirs = false;
let is_prnt_pr = false;
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_prm = body_prm::sun();
while step_mjd < end_mjd {
let step_zt = step_mjd;
let step_tt = step_zt - dzt + dt + tt_m_tai();
let step_tdb = step_tt + time::tdb_m_tt_lo(step_tt);
let mut earth_sun = sun_p(step_tdb, is_alg_hi);
let ligh_time = earth_sun.modul() / CLIGHT;
earth_sun = sun_p(step_tdb - ligh_time, is_alg_hi);
let nut = iau06_r::nut00b(step_tt);
let tr_equ = tr_equ(step_tt, nut, is_cirs);
let tr_ecl = iau06_r::gcrs_true_ecl(step_tt, nut.0);
let rot_prm = plngrph::sun(step_tt - ligh_time);
prm.push(Prm {
mjd_zt: step_zt,
sun_icrs: earth_sun,
icrs_true: tr_equ,
icrs_ecl: tr_ecl,
rot: rot_prm,
});
step_mjd += step;
}
for i in prm.iter() {
let sun_true = mat_r::mul_p(&i.icrs_true, &i.sun_icrs);
let sun_ecl = mat_r::mul_p(&i.icrs_ecl, &i.sun_icrs);
let sun_icrs_sph = sph_p::from_cart(&i.sun_icrs);
let sun_true_sph = sph_p::from_cart(&sun_true);
let sun_ecl_sph = sph_p::from_cart(&sun_ecl);
let cnstll = constell_name::constell_name(&i.sun_icrs);
let sun_diam = 2.0 * phys_prm.vis_equ_rad(sun_icrs_sph.rad);
let (lon_sep, lat_sep, pa) = i.rot.sub_obs(&sun_true, &i.icrs_true);
eph.push(Eph {
mjd_zt: i.mjd_zt,
lon: sun_ecl_sph.lon,
lat: sun_ecl_sph.lat,
rec: sun_true_sph.lon,
dec: sun_true_sph.lat,
dist: sun_icrs_sph.rad,
diam: sun_diam,
constell: cnstll,
lon_se: lon_sep,
lat_se: lat_sep,
pa_se: pa,
});
}
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 {
println!(
" {:10.6} {:+10.6} {} {:10.7} {:+10.6} {:10.8} {:6.2} {:6.2} {:+6.2} {:6.2}",
i.lon * RADDEG,
i.lat * RADDEG,
i.constell,
i.rec * RADDEG / 15.0,
i.dec * RADDEG,
i.dist,
i.diam * RADASEC,
i.lon_se * RADDEG,
i.lat_se * RADDEG,
angle::norm_dblpi(i.pa_se) * RADDEG,
);
} else {
println!(
" {:7.3} {:+7.3} {} {:7.4} {:+7.3} {:7.5} {:4.0} {:5.1} {:+5.1} {:5.1}",
i.lon * RADDEG,
i.lat * RADDEG,
i.constell,
i.rec * RADDEG / 15.0,
i.dec * RADDEG,
i.dist,
i.diam * RADASEC,
i.lon_se * RADDEG,
i.lat_se * RADDEG,
angle::norm_dblpi(i.pa_se) * RADDEG,
);
}
}
}
fn sun_p(mjd_tdb: f64, alg: bool) -> vec_p::VecP {
static EMRAT: f64 = 81.30057;
if alg {
let emb = vsop2013_hi::get_icrs_p(mjd_tdb, &BodyId::EMB);
let moon = elp85::get_icrs_p(mjd_tdb, &BodyId::MOON);
let earth = vec_p::add_scaled(&emb, -1.0 / (1.0 + EMRAT), &moon);
vec_p::scale(-1.0, &earth)
} else {
let emb = vsop2013_lo::get_icrs_p(mjd_tdb, &BodyId::EMB);
let moon = elp85::get_icrs_p(mjd_tdb, &BodyId::MOON);
let earth = vec_p::add_scaled(&emb, -1.0 / (1.0 + EMRAT), &moon);
vec_p::scale(-1.0, &earth)
}
}
fn tr_equ(mjd_tt: f64, nut: (f64, f64), is_cirs: bool) -> mat_r::MatR {
if is_cirs {
iau06_r::gcrs_cirs(mjd_tt, nut.0, nut.1)
} else {
iau06_r::gcrs_true(mjd_tt, nut.0, nut.1)
}
}