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;
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,
earth_sun: vec_p::VecP,
icrs_true: mat_r::MatR,
icrs_ecl: mat_r::MatR,
rot: PlnGrph,
}
struct Eph {
mjd_zt: f64,
lon: f64,
lat: f64,
ra: 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 = true;
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_prm = body_prm::sun();
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_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 icrs_true = tr_equ(step_tt, nut, is_cirs);
let icrs_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,
earth_sun,
icrs_true,
icrs_ecl,
rot: rot_prm,
});
step_mjd += step;
}
for i in prm.iter() {
let sun_true = mat_r::mul_p(&i.icrs_true, &i.earth_sun);
let sun_ecl = mat_r::mul_p(&i.icrs_ecl, &i.earth_sun);
let sun_icrs_sph = sph_p::from_cart(&i.earth_sun);
let sun_true_sph = sph_p::from_cart(&sun_true);
let sun_ecl_sph = sph_p::from_cart(&sun_ecl);
let constell = constell_name::constell_name(&i.earth_sun);
let 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,
ra: sun_true_sph.lon,
dec: sun_true_sph.lat,
dist: sun_icrs_sph.rad,
diam,
constell,
lon_se: lon_sep,
lat_se: lat_sep,
pa_se: pa,
});
}
println!();
println!("Body: {}", "Sun");
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 lon_dms = angle::rad_dms(i.lon, 2);
let lat_dms = angle::rad_dms(i.lat, 2);
let ra_hms = angle::rad_hms(i.ra, 3);
let dec_dms = angle::rad_dms(i.dec, 2);
let diam_ms = angle::rad_dms(i.diam, 2);
println!(
" {:#3}°{:#02}ʹ{:#02}.{:#02}ʺ {}{:#02}°{:#02}ʹ{:#02}.{:#02}ʺ {} {:#2}h{:#02}m{:#02}.{:#03}s {}{:#2}°{:#02}ʹ{:#02}.{:#02}ʺ {:10.8} {:#02}ʹ{:#02}.{:#02}ʺ {:6.2} {:+5.2} {:6.2}",
lon_dms.1,
lon_dms.2,
lon_dms.3,
lon_dms.4,
lat_dms.0,
lat_dms.1,
lat_dms.2,
lat_dms.3,
lat_dms.4,
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,
i.dist,
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.6} {:+10.6} {} {:10.7} {:+10.6} {:10.8} {:6.2} {:6.2} {:+5.2} {:6.2}",
i.lon * RADDEG,
i.lat * RADDEG,
i.constell,
i.ra * 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 {
if is_dms {
let lon_dms = angle::rad_dms(i.lon, -1);
let lat_dms = angle::rad_dms(i.lat, -1);
let ra_hms = angle::rad_hms(i.ra, 0);
let dec_dms = angle::rad_dms(i.dec, -1);
let diam_ms = angle::rad_dms(i.diam, 0);
println!(
" {:#3}°{:#02}.{}ʹ {}{:#2}°{:#02}.{}ʹ {} {:#2}h{:#02}m{:#02}s {}{:#2}°{:#02}.{}ʹ {:7.5} {:#02}ʹ{:#02}ʺ {:5.1} {:+4.1} {:5.1}",
lon_dms.1,
lon_dms.2,
lon_dms.3 / 6,
lat_dms.0,
lat_dms.1,
lat_dms.2,
lat_dms.3 / 6,
i.constell,
ra_hms.1,
ra_hms.2,
ra_hms.3,
dec_dms.0,
dec_dms.1,
dec_dms.2,
dec_dms.3 / 6,
i.dist,
diam_ms.2,
diam_ms.3,
i.lon_se * RADDEG,
i.lat_se * RADDEG,
angle::norm_dblpi(i.pa_se) * RADDEG,
);
} else {
println!(
" {:7.3}° {:+7.3}° {} {:7.4}h {:+7.3}° {:7.5} {:4.0}ʺ {:5.1} {:+4.1} {:5.1}",
i.lon * RADDEG,
i.lat * RADDEG,
i.constell,
i.ra * 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 {
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
}
}