extern crate rkepler;

use rkepler::datetime::mjd;
use rkepler::frame::constant::*;
use rkepler::frame::iau06::iau06_r;
use rkepler::math::angle::*;
use rkepler::math::sphe::sph;
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;

#[test]
fn t_constel_name() {
    //10Uma Lyn
    let mut star_j2000 = sph::Sph {
        lon: hms_rad(9.0, 0.6, 0.0),
        lat: dms_rad('+', 41.0, 47.0, 0.0),
    };
    let mut star_j2000_p = star_j2000.to_cart();
    let mut constell = constell_name::constell_name(&star_j2000_p);
    assert_eq!(constell, "Lyn".to_string());
    //alph CMa
    star_j2000 = sph::Sph {
        lon: hms_rad(6.0, 45.1, 0.0),
        lat: dms_rad('-', 16.0, 43.0, 0.0),
    };
    star_j2000_p = star_j2000.to_cart();
    constell = constell_name::constell_name(&star_j2000_p);
    assert_eq!(constell, "CMa".to_string());
    //51 And
    star_j2000 = sph::Sph {
        lon: hms_rad(1.0, 38.0, 0.0),
        lat: dms_rad('+', 48.0, 38.0, 0.0),
    };
    star_j2000_p = star_j2000.to_cart();
    constell = constell_name::constell_name(&star_j2000_p);
    assert_eq!(constell, "And".to_string());
}

#[test]
fn t_sun() {
    let mjd_tt = mjd::jd_mjd(2460000.5);
    let sph_es = sph_p::SphP {
        lon: hms_rad(22.0, 31.0, 31.2914),
        lat: dms_rad('-', 9.0, 16.0, 17.259),
        rad: 0.989713972019,
    };

    let earth_sun = sph_es.to_cart();
    let ligtime_es = sph_es.rad / CLIGHT;
    let ligtime_es_m = ligtime_es * 1440.0;
    assert!(
        (ligtime_es_m - 8.2312).abs() < 1e-6,
        "Light time: {}",
        ligtime_es_m
    );

    let rot_prm = plngrph::sun(mjd_tt - ligtime_es);
    let ra0p = rot_prm.rap * RADDEG;
    let de0p = rot_prm.dep * RADDEG;
    let w0 = norm_dblpi(rot_prm.w0) * RADDEG;
    assert!((ra0p - 286.13).abs() < 1e-7, "Ra0: {}", ra0p);
    assert!((de0p - 63.87).abs() < 1e-7, "De0: {}", de0p);
    assert!((w0 - 140.28912).abs() < 1e-5, "W0: {}", w0);

    let (dpsi, deps) = iau06_r::nut00b(mjd_tt);
    let bpn = iau06_r::gcrs_true(mjd_tt, dpsi, deps);
    let (lon_sep, lat_sep, pa) = rot_prm.sub_obs(&earth_sun, &bpn);
    let lon_sep_d = lon_sep * RADDEG;
    let lat_sep_d = lat_sep * RADDEG;
    let pa_d = norm_dblpi(pa) * RADDEG;
    assert!((lon_sep_d - 1.03).abs() < 1e-2, "Lon sep: {}", lon_sep_d);
    assert!((lat_sep_d - -7.14).abs() < 1e-2, "Lat sep: {}", lat_sep_d);
    assert!((pa_d - 339.64).abs() < 1e-2, "PA: {}", pa_d);

    let phys_prm = body_prm::sun();
    let f = phys_prm.flatt_pol;
    assert!((f - 0.0).abs() < 1e-6, "f: {}", f);
    let reqv = phys_prm.vis_equ_rad(sph_es.rad);
    let reqv_s = reqv * RADASEC;
    assert!((reqv_s - 969.19).abs() < 1e-3, "Reqv: {}", reqv_s);
    let dpole = plngrph::d_pol(reqv_s, lat_sep);
    assert!((dpole - -961.67).abs() < 1e-3, "dPole: {}", dpole);

    let vmag = phys_prm.vis_mag(sph_es.rad, 1.0);
    assert!((vmag - -26.80).abs() < 1e-2, "Magnetude: {}", vmag); //-26.76
    let sbrt = phys_eph::surf_bright(vmag, reqv_s, reqv_s, 1.0);
    assert!((sbrt - -10.63).abs() < 1e-2, "S-brt: {}", sbrt); //-10.589
}

#[test]
fn t_venus() {
    let mjd_tt = mjd::jd_mjd(2460000.5);
    let sph_es = sph_p::SphP {
        lon: hms_rad(22.0, 31.0, 31.2914),
        lat: dms_rad('-', 9.0, 16.0, 17.259),
        rad: 0.989713972019,
    };
    let sph_sp = sph_p::SphP {
        lon: hms_rad(3.0, 5.0, 33.3726),
        lat: dms_rad('+', 15.0, 44.0, 26.125),
        rad: 0.722727967638,
    };
    let sph_ep = sph_p::SphP {
        lon: hms_rad(0.0, 22.0, 24.0626),
        lat: dms_rad('+', 1.0, 30.0, 19.983),
        rad: 1.391939037454,
    };

    let earth_sun = sph_es.to_cart();
    let sun_pln = sph_sp.to_cart();
    let earth_pln = sph_ep.to_cart();
    let ligtime_sp = sph_sp.rad / CLIGHT;
    let ligtime_sp_m = ligtime_sp * 1440.0;
    let ligtime_ep = sph_ep.rad / CLIGHT;
    let ligtime_ep_m = ligtime_ep * 1440.0;
    assert!(
        (ligtime_sp_m - 6.010745).abs() < 1e-6,
        "Light time: {}",
        ligtime_sp_m
    );
    assert!(
        (ligtime_ep_m - 11.576404).abs() < 1e-6,
        "Light time: {}",
        ligtime_ep_m
    );

    let rot_prm = plngrph::venus(mjd_tt - ligtime_ep);
    let ra0p = rot_prm.rap * RADDEG;
    let de0p = rot_prm.dep * RADDEG;
    let w0 = rot_prm.w0 * RADDEG;
    assert!((ra0p - 272.76).abs() < 1e-7, "Ra0: {}", ra0p);
    assert!((de0p - 67.16).abs() < 1e-7, "De0: {}", de0p);
    assert!((w0 - -125.5019794).abs() < 1e-7, "W0: {}", w0);

    let (dpsi, deps) = iau06_r::nut00b(mjd_tt);
    let bpn = iau06_r::gcrs_true(mjd_tt, dpsi, deps);
    let (lon_sep, lat_sep, pa) = rot_prm.sub_obs(&earth_pln, &bpn);
    let lon_sep_d = lon_sep * RADDEG;
    let lat_sep_d = lat_sep * RADDEG;
    let pa_d = norm_dblpi(pa) * RADDEG;
    assert!((lon_sep_d - 308.71).abs() < 1e-2, "Lon sep: {}", lon_sep_d);
    assert!((lat_sep_d - -0.28).abs() < 1e-2, "Lat sep: {}", lat_sep_d);
    assert!((pa_d - 337.19).abs() < 1e-2, "PA: {}", pa_d);

    let (lon_ssp, lat_ssp, pas) = rot_prm.sub_obs(&sun_pln, &bpn);
    let lon_sep_d = lon_ssp * RADDEG;
    let lat_sep_d = lat_ssp * RADDEG;
    let pa_d = norm_dblpi(pas) * RADDEG;
    assert!((lon_sep_d - 351.34).abs() < 1e-2, "Lon sep: {}", lon_sep_d);
    assert!((lat_sep_d - 0.45).abs() < 1e-2, "Lat sep: {}", lat_sep_d);
    assert!((pa_d - 343.69).abs() < 1e-2, "PA: {}", pa_d);

    let phys_prm = body_prm::venus();
    let f = phys_prm.flatt_pol;
    assert!((f - 0.0).abs() < 1e-6, "f: {}", f);
    let reqv = phys_prm.vis_equ_rad(sph_ep.rad);
    let reqv_s = reqv * RADASEC;
    assert!((reqv_s - 5.995).abs() < 1e-3, "Reqv: {}", reqv_s);
    let dpole = plngrph::d_pol(reqv_s, lat_sep);
    assert!((dpole - -5.994).abs() < 1e-3, "dPole: {}", dpole);

    let elong = phys_eph::elongation(&earth_pln, &earth_sun) * RADDEG;
    assert!((elong - 29.6419).abs() < 1e-3, "Elongation: {}", elong);
    let pha_ang = phys_eph::phase_angl(&sun_pln, &earth_pln);
    let pha_ang_d = pha_ang * RADDEG;
    assert!(
        (pha_ang_d - 42.6415).abs() < 1e-3, //42.6353
        "Phase angle: {}",
        pha_ang_d
    );
    let phase = phys_eph::phase(pha_ang);
    assert!((phase - 0.868).abs() < 1e-3, "Phase: {}", phase);

    let vmag = body_prm::vis_mag_venus(sph_sp.rad, sph_ep.rad, pha_ang);
    assert!((vmag - -3.964).abs() < 1e-2, "Magnetude: {}", vmag); //-3.934
    let sbrt = phys_eph::surf_bright(vmag, reqv_s, reqv_s, phase);
    assert!((sbrt - 1.013).abs() < 1e-2, "S-brt: {}", sbrt); //1.044
}

#[test]
fn t_saturn() {
    let mjd_tt = mjd::jd_mjd(2460000.5);
    let earth_sun = vec_p::VecP {
        x: 0.904892857278,
        y: -0.367799882801,
        z: -0.159455039018,
    };
    let sun_pln = vec_p::VecP {
        x: 8.327302479796,
        y: -4.677826508638,
        z: -2.288720156389,
    };
    let earth_pln = vec_p::VecP {
        x: 9.231660146091,
        y: -5.046416175108,
        z: -2.448514001383,
    };

    let r_es = earth_sun.modul();
    let r_sp = sun_pln.modul();
    let r_ep = earth_pln.modul();
    let ligtime_sp = r_sp / CLIGHT;
    let ligtime_sp_m = ligtime_sp * 1440.0;
    let ligtime_ep = r_ep / CLIGHT;
    let ligtime_ep_m = ligtime_ep * 1440.0;
    assert!((r_es - 0.989713972019).abs() < 1e-7, "Res: {}", r_es);
    assert!((r_sp - 9.821622441796).abs() < 1e-7, "Rsp: {}", r_sp);
    assert!((r_ep - 10.802087116862).abs() < 1e-7, "Rep: {}", r_ep);
    assert!(
        (ligtime_sp_m - 81.683943).abs() < 1e-6,
        "Light time: {}",
        ligtime_sp_m
    );
    assert!(
        (ligtime_ep_m - 89.83822).abs() < 1e-6,
        "Light time: {}",
        ligtime_ep_m
    );

    let rot_prm = plngrph::saturn(mjd_tt - ligtime_ep);
    let ra0p = rot_prm.rap * RADDEG;
    let de0p = rot_prm.dep * RADDEG;
    let w0 = rot_prm.w0 * RADDEG;
    assert!((ra0p - 40.5806661).abs() < 1e-7, "Ra0: {}", ra0p);
    assert!((de0p - 83.5360740).abs() < 1e-7, "De0: {}", de0p);
    assert!((w0 - 176.1582152).abs() < 1e-7, "W0: {}", w0);

    let (dpsi, deps) = iau06_r::nut00b(mjd_tt);
    let bpn = iau06_r::gcrs_true(mjd_tt, dpsi, deps);
    let (lon_sep, lat_sep, pa) = rot_prm.sub_obs(&earth_pln, &bpn);
    let lon_sep_d = lon_sep * RADDEG;
    let lat_sep_d = lat_sep * RADDEG;
    let pa_d = norm_dblpi(pa) * RADDEG;
    assert!((lon_sep_d - 204.84).abs() < 1e-2, "Lon sep: {}", lon_sep_d);
    assert!((lat_sep_d - 10.88).abs() < 1e-2, "Lat sep: {}", lat_sep_d);
    assert!((pa_d - 6.10).abs() < 1e-2, "PA: {}", pa_d);

    let phys_prm = body_prm::saturn();
    let lat_plngr = plngrph::plngr_lat(phys_prm.flatt_pol, lat_sep);
    let lat_plngr_d = lat_plngr * RADDEG;
    assert!(
        (lat_plngr_d - 13.29).abs() < 1e-2,
        "Lat plngr: {}",
        lat_plngr_d
    );

    let (lon_ssp, lat_ssp, pas) = rot_prm.sub_obs(&sun_pln, &bpn);
    let lon_sep_d = lon_ssp * RADDEG;
    let lat_sep_d = lat_ssp * RADDEG;
    let pa_d = norm_dblpi(pas) * RADDEG;
    assert!((lon_sep_d - 204.23).abs() < 1e-2, "Lon sep: {}", lon_sep_d);
    assert!((lat_sep_d - 11.32).abs() < 1e-2, "Lat sep: {}", lat_sep_d);
    assert!((pa_d - 6.135).abs() < 1e-2, "PA: {}", pa_d);

    let f = phys_prm.flatt_pol;
    assert!((f - 0.0979624).abs() < 1e-6, "f: {}", f);
    let reqv = phys_prm.vis_equ_rad(r_ep);
    let reqv_s = reqv * RADASEC;
    assert!((reqv_s - 7.6924).abs() < 1e-3, "Reqv: {}", reqv_s);
    let rpol_s = phys_prm.vis_pol_rad(reqv_s, lat_sep);
    //let rpol_s = rpol * RADASEC;
    assert!((rpol_s - 6.966).abs() < 1e-3, "Rpol: {}", rpol_s);
    let dpole = plngrph::d_pol(rpol_s, lat_sep);
    assert!((dpole - 6.84).abs() < 1e-3, "dPole: {}", dpole); //6.81

    let elong = phys_eph::elongation(&earth_pln, &earth_sun) * RADDEG;
    assert!((elong - 7.4733).abs() < 1e-3, "Elongation: {}", elong);
    let pha_ang = phys_eph::phase_angl(&sun_pln, &earth_pln);
    let pha_ang_d = pha_ang * RADDEG;
    assert!(
        (pha_ang_d - 0.7453).abs() < 1e-3,
        "Phase angle: {}",
        pha_ang_d
    );
    let phase = phys_eph::phase(pha_ang);
    assert!((phase - 1.000).abs() < 1e-3, "Phase: {}", phase);

    let vmag = body_prm::vis_mag_saturn(r_sp, r_ep, pha_ang, lat_sep);
    assert!((vmag - 0.83).abs() < 1e-2, "Magnetude: {}", vmag); //1.21,  0.869
    let sbrt = phys_eph::surf_bright(vmag, reqv_s, reqv_s, phase);
    assert!((sbrt - 6.508).abs() < 1e-2, "S-brt: {}", sbrt); //6.537
}

#[test]
fn t_moon() {
    let date = mjd::Adate {
        year: 2023,
        month: 2,
        day: 25,
        hour: 0,
        min: 0,
        sec: 0,
        msec: 0,
    };
    let mjd_tt = date.mjd();
    let earth_sun = vec_p::VecP {
        x: 0.904892857278,
        y: -0.367799882801,
        z: -0.159455039018,
    };
    let sun_pln = vec_p::VecP {
        x: -0.902895674191,
        y: 0.369255074186,
        z: 0.16009761677,
    };
    let earth_pln = vec_p::VecP {
        x: 0.0019943952,
        y: 0.001458294977,
        z: 0.000644251044,
    };

    let r_es = earth_sun.modul();
    let r_sp = sun_pln.modul();
    let r_ep = earth_pln.modul();
    let ligtime_ep = r_ep / CLIGHT;
    let ligtime_ep_m = ligtime_ep * 1440.0;
    assert!((r_es - 0.989713972019).abs() < 1e-7, "Res: {}", r_es);
    assert!((r_sp - 0.98853485279).abs() < 1e-7, "Rsp: {}", r_sp);
    assert!((r_ep - 0.002553291182).abs() < 1e-7, "Rep: {}", r_ep);
    assert!(
        (ligtime_ep_m - 0.021235).abs() < 1e-6,
        "Light time: {}",
        ligtime_ep_m
    );

    let rot_prm = plngrph::moon(mjd_tt - ligtime_ep);
    let ra0p = rot_prm.rap * RADDEG;
    let de0p = rot_prm.dep * RADDEG;
    let w0 = rot_prm.w0 * RADDEG;
    assert!((ra0p - 267.6206851).abs() < 1e-7, "Ra0: {}", ra0p);
    assert!((de0p - 67.7971485).abs() < 1e-7, "De0: {}", de0p);
    assert!((w0 - 213.2266121).abs() < 1e-7, "W0: {}", w0);

    let (dpsi, deps) = iau06_r::nut00b(mjd_tt);
    let bpn = iau06_r::gcrs_true(mjd_tt, dpsi, deps);
    let (lon_sep, lat_sep, pa) = rot_prm.sub_obs(&earth_pln, &bpn);
    let lon_sep_d = lon_sep * RADDEG;
    let lat_sep_d = lat_sep * RADDEG;
    let pa_d = norm_dblpi(pa) * RADDEG;
    assert!((lon_sep_d - 7.61).abs() < 1e-2, "Lon sep: {}", lon_sep_d);
    assert!((lat_sep_d - -0.32).abs() < 1e-2, "Lat sep: {}", lat_sep_d);
    assert!((pa_d - 342.81).abs() < 1e-2, "PA: {}", pa_d);

    let (lon_ssp, lat_ssp, pas) = rot_prm.sub_obs(&sun_pln, &bpn);
    let lon_sep_d = lon_ssp * RADDEG;
    let lat_sep_d = lat_ssp * RADDEG;
    let pa_d = norm_dblpi(pas) * RADDEG;
    assert!((lon_sep_d - 124.97).abs() < 1e-2, "Lon sep: {}", lon_sep_d);
    assert!((lat_sep_d - -1.34).abs() < 1e-2, "Lat sep: {}", lat_sep_d);
    assert!((pa_d - 20.83).abs() < 1e-2, "PA: {}", pa_d);

    let phys_prm = body_prm::moon();
    let reqv = phys_prm.vis_equ_rad(r_ep);
    let reqv_s = reqv * RADASEC;
    assert!((reqv_s - 938.202).abs() < 1e-3, "Reqv: {}", reqv_s);
    let dpole = plngrph::d_pol(reqv_s, lat_sep);
    assert!((dpole - -938.187).abs() < 1e-3, "dPole: {}", dpole);

    let elong = phys_eph::elongation(&earth_pln, &earth_sun) * RADDEG;
    assert!((elong - 62.5314).abs() < 1e-3, "Elongation: {}", elong);
    let pha_ang = phys_eph::phase_angl(&sun_pln, &earth_pln);
    let pha_ang_d = pha_ang * RADDEG;
    assert!(
        (pha_ang_d - 117.3374).abs() < 1e-3, //117.3431
        "Phase angle: {}",
        pha_ang_d
    );
    let phase = phys_eph::phase(pha_ang);
    assert!((phase - 0.270).abs() < 1e-3, "Phase: {}", phase);

    let vmag = body_prm::vis_mag_moon(r_sp, r_ep, pha_ang);
    assert!((vmag - -8.69).abs() < 1e-2, "Magnetude: {}", vmag); //-8.89
    let sbrt = phys_eph::surf_bright(vmag, reqv_s, reqv_s, phase);
    assert!((sbrt - 5.997).abs() < 1e-2, "S-brt: {}", sbrt); //5.734
}
