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)
    }
}
