use neptunes::biosphere::aerochemistry::{
    cloudlayerchemosynthesisefficiency, hypotheticalchemotroph,
};
use neptunes::biosphere::ecosystems::hypotheticalatmosphericbiome;
use neptunes::biosphere::fauna::{
    Aeroplankton, CloudLayerInteraction, hypotheticalchemoautotroph, hypotheticalh2floater,
};
use neptunes::geodata::bathymetry::InteriorLayerData;
use neptunes::geodata::coordinates::{
    LatLon, NEPTUNEFLATTENING, NEPTUNESEMIMAJORM, NEPTUNESEMIMINORM,
};
use neptunes::geodata::elevation::ElevationProvider;
use neptunes::geodata::regions::{RegionDatabase, RegionType};
use neptunes::lighting::day_night::{DayNightCycle, DaylightState};
use neptunes::lighting::seasons::{
    AXIALTILTDEG as SEASONSTILT, NEPTUNEORBITALPERIODDAYS, Season, VERNALEQUINOXJD, seasonat,
    subsolarpoint,
};
use neptunes::lighting::solar_position::{NEPTUNEAXIALTILTDEG, SolarPosition};
use neptunes::physics::collisions::{cometimpact, kuiperbeltimpactor};
use neptunes::physics::orbit::{
    ARGUMENTPERIHELIONDEG, ECCENTRICITY, INCLINATIONDEG, LONGITUDEASCENDINGNODEDEG, NeptuneOrbit,
    SEMIMAJORAXIS,
};
use neptunes::physics::rotation::{
    AXIALTILTDEG, AXIALTILTRAD, NeptuneRotation, PRECESSIONPERIODYEARS, SIDEREALDAYS,
};
use neptunes::physics::tides::{
    TidalForce, neaptideamplitude, springtideamplitude, tritontosolartideratio,
};
use neptunes::rendering::atmosphere_scattering::AtmosphereParams;
use neptunes::rendering::clouds::{CloudLayer, CloudSystem};
use neptunes::rendering::materials::PbrMaterial;
use neptunes::rendering::ocean_rendering::OceanParams;
use neptunes::rendering::shaders::ShaderData;
use neptunes::satellites::artificial::{ArtificialSatellite, Constellation, OrbitType};
use neptunes::satellites::triton::{
    NEPTUNETRITONDISTANCE, TRITONMASS, TRITONORBITALPERIOD, TRITONRADIUS, TritonSource, TritonState,
};
use neptunes::temporal::calendar::{
    DateTime, NEPTUNEDAYSECONDS, NEPTUNEYEARSECONDS, SECONDSPERDAY, UNIXEPOCHJD,
};
use neptunes::temporal::epoch::{Epoch, J1950EPOCH, J2000EPOCH, MJDOFFSET};
use neptunes::temporal::time_scale::TimeScale;
use neptunes::terrain::heightmap::Heightmap;
use neptunes::terrain::lod::{LodConfig, LodTerrain};
use neptunes::terrain::mesh::TerrainMesh;
use neptunes::terrain::texturing::CloudBiomeClassifier;

fn main() {
    let orbit = NeptuneOrbit::new();
    let perioddays = orbit.orbitalperioddays();
    let vperihelion = orbit.velocityatdistance(orbit.perihelionm());
    let vaphelion = orbit.velocityatdistance(orbit.aphelionm());
    let energy = orbit.specific_orbital_energy();
    let angmom = orbit.specific_angular_momentum();
    let vesc = NeptuneOrbit::escape_velocity_at_surface();
    let fsun = orbit.gravitational_force_sun();
    let r = orbit.current_radius();
    let vmean = orbit.mean_orbital_velocity();
    println!(
        "Orbit: P={:.2}d  vperi={:.0}m/s  vaph={:.0}m/s  vmean={:.0}m/s",
        perioddays, vperihelion, vaphelion, vmean
    );
    println!(
        "  E={:.2e}J/kg  L={:.2e}m²/s  vesc={:.0}m/s  Fsun={:.2e}N  r={:.3e}m",
        energy, angmom, vesc, fsun, r
    );
    println!(
        "  a={:.3e}m  e={}  i={}°  Ω={}°  ω={}°",
        SEMIMAJORAXIS,
        ECCENTRICITY,
        INCLINATIONDEG,
        LONGITUDEASCENDINGNODEDEG,
        ARGUMENTPERIHELIONDEG
    );

    let rot = NeptuneRotation::new();
    let vequator = rot.surfacevelocityatlatitude(0.0);
    let v45 = rot.surfacevelocityatlatitude(45.0);
    let accel = rot.centripetalaccelerationatlatitude(0.0);
    let coriolis = rot.coriolisparameter(45.0);
    let moi = rot.momentofinertia();
    let ke = rot.rotationalkineticenergy();
    let lrot = rot.angular_momentum();
    let prec = rot.precession_rate_rad_per_year();
    let daylensummer = rot.day_length_variation_due_to_tilt(172, 60.0);
    let daylenwinter = rot.day_length_variation_due_to_tilt(355, 60.0);
    println!(
        "Rotation: veq={:.0}m/s  v45={:.0}m/s  ac={:.4}m/s²  f45={:.5e}",
        vequator, v45, accel, coriolis
    );
    println!(
        "  I={:.3e}kg·m²  KE={:.3e}J  L={:.3e}kg·m²/s  prec={:.2e}rad/yr",
        moi, ke, lrot, prec
    );
    println!(
        "  day@60°N: summer={:.1}h  winter={:.1}h",
        daylensummer, daylenwinter
    );
    println!(
        "  sidereal={:.4}s  tilt={}°={:.4}rad  precperiod={}yr",
        SIDEREALDAYS, AXIALTILTDEG, AXIALTILTRAD, PRECESSIONPERIODYEARS
    );

    let tritontide = TidalForce::fromtriton();
    let suntide = TidalForce::fromsun();
    let atriton = tritontide.tidalacceleration();
    let asun = suntide.tidalacceleration();
    let pot = tritontide.tidalpotential(0.0);
    let bulgetriton = tritontide.tidalbulgeheight();
    let bulgesun = suntide.tidalbulgeheight();
    let grav = tritontide.gravitationalattraction();
    let spring = springtideamplitude();
    let neap = neaptideamplitude();
    let ratio = tritontosolartideratio();
    println!(
        "Tides: aTriton={:.2e}  aSun={:.2e}  ratio={:.2}",
        atriton, asun, ratio
    );
    println!(
        "  bulgeTriton={:.3}m  bulgeSun={:.3}m  spring={:.3}m  neap={:.3}m",
        bulgetriton, bulgesun, spring, neap
    );
    println!("  potential(0)={:.2e}  Fgrav={:.2e}N", pot, grav);

    let comet = cometimpact();
    let kuiper = kuiperbeltimpactor();
    let kecomet = comet.kineticenergymt();
    let crater = comet.craterdiameterm(1600.0);
    let fireball = comet.fireballradiusm();
    let ejecta = comet.ejectavolumem3(1600.0);
    let vimp = comet.impactvelocity();
    let kekuiper = kuiper.kineticenergymt();
    println!(
        "Comet: {:.2e}Mt  crater={:.0}m  fireball={:.0}m  ejecta={:.2e}m³  vimp={:.0}m/s",
        kecomet, crater, fireball, ejecta, vimp
    );
    println!("Kuiper: {:.2e}Mt", kekuiper);

    let mut triton = TritonState::new();
    let source = match triton.source.get() {
        TritonSource::Binary => "binary",
        TritonSource::Simulation => "simulated",
    };
    let pos0 = triton.position();
    triton.step(3600.0);
    let pos1 = triton.position();
    let gtriton = triton.gravityat(NEPTUNETRITONDISTANCE);
    let dtriton = triton.distance();
    println!(
        "Triton({}): pos0=({:.0},{:.0},{:.0})  pos1h=({:.0},{:.0},{:.0})",
        source, pos0.0, pos0.1, pos0.2, pos1.0, pos1.1, pos1.2
    );
    println!("  g@dist={:.4e}m/s²  dist={:.3e}m", gtriton, dtriton);
    println!(
        "  M={:.3e}kg  R={:.4e}m  d={:.3e}m  P={:.0}s",
        TRITONMASS, TRITONRADIUS, NEPTUNETRITONDISTANCE, TRITONORBITALPERIOD
    );

    let low = ArtificialSatellite::loworbit("NepLow", 500.0, 500000.0);
    let polar = ArtificialSatellite::polarorbit("NepPolar", 300.0, 800000.0);
    let custom = ArtificialSatellite::new("NepHEO", 1500.0, 2000000.0, 0.5, 1.0);
    let plow = low.orbitalperiods();
    let vlow = low.orbitalvelocityms();
    let ppolar = polar.orbitalperiods();
    let gsurf = low.gravityatsurface();
    let poslow = low.position();
    let pospolar = polar.position();

    let orbittype = |sat: &ArtificialSatellite| -> &str {
        match sat.orbittype {
            OrbitType::Low => "Low",
            OrbitType::Medium => "Medium",
            OrbitType::High => "High",
            OrbitType::Custom => "Custom",
        }
    };
    println!(
        "Low({}): P={:.0}s  v={:.0}m/s  gsurf={:.2}  pos=({:.0},{:.0},{:.0})",
        orbittype(&low),
        plow,
        vlow,
        gsurf,
        poslow.0,
        poslow.1,
        poslow.2
    );
    println!(
        "Polar({}): P={:.0}s  pos=({:.0},{:.0},{:.0})",
        orbittype(&polar),
        ppolar,
        pospolar.0,
        pospolar.1,
        pospolar.2
    );
    println!(
        "Custom({}): e={:.1}  i={:.1}rad",
        orbittype(&custom),
        custom.eccentricity,
        custom.inclinationrad
    );

    let mut constellation = Constellation::new("NepNet");
    constellation.add(ArtificialSatellite::loworbit("N1", 260.0, 600000.0));
    constellation.add(ArtificialSatellite::loworbit("N2", 260.0, 600000.0));
    let mut stepped = low;
    stepped.step(3600.0);
    let posstepped = stepped.position();
    constellation.stepall(1800.0);
    let cpositions = constellation.positions();
    let ccount = cpositions.len();
    println!(
        "Stepped: pos=({:.0},{:.0},{:.0})  constellation: {} sats",
        posstepped.0, posstepped.1, posstepped.2, ccount
    );

    let dnc = DayNightCycle::new(2451545.0);
    let state = dnc.stateat(0.0, 0.0);
    let ambient = dnc.ambientlight(45.0, 90.0);
    let terminator = dnc.terminatorpoints(36);
    let tcount = terminator.len();
    let statelabel = match state {
        DaylightState::Day => "Day",
        DaylightState::CivilTwilight => "Civil",
        DaylightState::NauticalTwilight => "Nautical",
        DaylightState::AstronomicalTwilight => "Astro",
        DaylightState::Night => "Night",
    };
    println!(
        "DayNight: state={}  ambient(45,90)={:.3}  terminator={}pts",
        statelabel, ambient, tcount
    );

    let sp = SolarPosition::compute(2451545.0, 30.0, 0.0);
    let sdist = sp.distancem();
    let above = sp.isabovehorizon();
    let flux = sp.solarfluxwm2();
    println!(
        "Solar: dist={:.3e}m  above={}  flux={:.4}W/m²  tilt={}°",
        sdist, above, flux, NEPTUNEAXIALTILTDEG
    );

    let ss = seasonat(2451545.0, 45.0);
    let (sublat, sublon) = subsolarpoint(2451545.0);
    let northlabel = match ss.seasonnorth {
        Season::Spring => "Spring",
        Season::Summer => "Summer",
        Season::Autumn => "Autumn",
        Season::Winter => "Winter",
    };
    println!(
        "Season@45°N: {}  decl={:.2}°  daylen={:.2}h  subsolar=({:.2},{:.2})",
        northlabel, ss.solardeclinationdeg, ss.daylengthhours, sublat, sublon
    );
    println!(
        "  tilt={}°  Porbit={}d  vernaleq={}JD",
        SEASONSTILT, NEPTUNEORBITALPERIODDAYS, VERNALEQUINOXJD
    );

    let dt = DateTime::new(2025, 4, 3, 12, 0, 0.0);
    let jd = dt.tojuliandate();
    let nd = dt.toneptuneday();
    let ny = dt.toneptuneyear();
    let unix = dt.tounixtimestamp();
    let dt2 = DateTime::fromjuliandate(jd);
    let dt3 = DateTime::fromunixtimestamp(unix);
    println!(
        "DateTime: {}-{}-{} JD={:.2}  NepDay={:.4}  NepYear={:.6}",
        dt.year, dt.month, dt.day, jd, nd, ny
    );
    println!(
        "  roundtrip JD: {}-{}-{}  from unix: {}-{}-{}",
        dt2.year, dt2.month, dt2.day, dt3.year, dt3.month, dt3.day
    );
    println!(
        "  constants: UNIX_JD={} SPD={} NepDay={}s NepYear={:.2e}s",
        UNIXEPOCHJD, SECONDSPERDAY, NEPTUNEDAYSECONDS, NEPTUNEYEARSECONDS
    );

    let ep = Epoch::j2000();
    let ep2 = Epoch::fromjd(2460000.0);
    let mjd = ep2.tomjd();
    let cent = ep2.centuriessincej2000();
    let gmst = ep2.neptunegmstdegrees();
    println!(
        "Epoch: J2000={:.1}JD  ep2={:.1}JD  MJD={:.1}  centuries={:.4}  GMST={:.2}°",
        ep.juliandate, ep2.juliandate, mjd, cent, gmst
    );
    println!(
        "  constants: J2000={} J1950={} MJDoff={}",
        J2000EPOCH, J1950EPOCH, MJDOFFSET
    );

    let mut ts = TimeScale::realtime();
    ts.step(86400.0);
    let simhrs = ts.simhours();
    let simdays = ts.simdays();
    let simnepdays = ts.simneptunedays();
    let simnepyrs = ts.simneptuneyears();

    let mut ff = TimeScale::fastforward(10.0);
    ff.step(1000.0);
    let ffhrs = ff.simhours();
    let simdt = ff.simulationdt(1.0);

    let mut slow = TimeScale::slowmotion(0.1);
    slow.step(10000.0);
    slow.pause();
    slow.resume();
    slow.togglepause();
    slow.setspeed(2.0);
    let slowhrs = slow.simhours();
    println!(
        "TimeScale: real={}h/{}d/{:.4}Nd/{:.8}Ny  ff={}h  simdt={:.1}  slow={}h",
        simhrs, simdays, simnepdays, simnepyrs, ffhrs, simdt, slowhrs
    );

    let atm = AtmosphereParams::default();
    let rayleigh = atm.rayleighdensity(50000.0);
    let mie = atm.miedensity(10000.0);
    let sray = atm.scatterrayleigh(0.5);
    let smie = atm.scattermie(0.5);
    let sky = atm.skycolor([0.0, 1.0, 0.0], [0.0, 0.5, 0.866], 0.0);
    println!(
        "Atmo: rayleigh@50km={:.4}  mie@10km={:.4}  scatR={:.4}  scatM={:.4}",
        rayleigh, mie, sray, smie
    );
    println!("  sky=({:.3},{:.3},{:.3})", sky[0], sky[1], sky[2]);

    let ch4 = CloudLayer::methaneice();
    let h2s = CloudLayer::hydrogensulfide();
    let nh4sh = CloudLayer::ammoniahydrosulfide();
    let deep = CloudLayer::deepwater();
    let haze = CloudLayer::highaltitudehaze();
    println!(
        "Clouds: CH4={:?}@{:.0}m  H2S={:?}@{:.0}m  NH4SH={:?}@{:.0}m",
        ch4.cloudtype,
        ch4.basealtitudem,
        h2s.cloudtype,
        h2s.basealtitudem,
        nh4sh.cloudtype,
        nh4sh.basealtitudem
    );
    println!(
        "  Deep={:?}@{:.0}m  Haze={:?}@{:.0}m",
        deep.cloudtype, deep.basealtitudem, haze.cloudtype, haze.basealtitudem
    );

    let mut cs = CloudSystem::neptunedefault();
    cs.step(3600.0);
    let d1 = cs.sampledensity(50000.0, 30.0, 90.0);
    let d2 = cs.sampledensity(100000.0, -45.0, 180.0);
    println!(
        "CloudSystem: density@50km={:.4}  density@100km={:.4}",
        d1, d2
    );

    let mch4 = PbrMaterial::methanecloud();
    let mh2s = PbrMaterial::h2scloud();
    let mdeep = PbrMaterial::deepatmosphere();
    let mhaze = PbrMaterial::hazelayer();
    let mrock = PbrMaterial::corerock();
    let msup = PbrMaterial::superionicice();
    let mstorm = PbrMaterial::stormregion();
    println!(
        "PBR: CH4 rough={:.2} metal={:.2}  H2S rough={:.2}  deep rough={:.2}",
        mch4.roughness, mch4.metallic, mh2s.roughness, mdeep.roughness
    );
    println!(
        "  haze rough={:.2}  rock metal={:.2}  supIce rough={:.2}  storm rough={:.2}",
        mhaze.roughness, mrock.metallic, msup.roughness, mstorm.roughness
    );

    let ocean = OceanParams::default();
    let spec = ocean.generatespectrum();
    let hcount = spec.heights.len();
    let ncount = spec.normals.len();
    println!(
        "Ocean: spectrum grid={}  heights={}  normals={}",
        spec.gridsize, hcount, ncount
    );

    let terrain_shader = ShaderData::terrain();
    let atm_shader = ShaderData::atmosphere();
    let ocean_shader = ShaderData::interiorocean();
    let tu_count = terrain_shader.uniforms.len();
    let au_count = atm_shader.uniforms.len();
    let ou_count = ocean_shader.uniforms.len();
    println!(
        "Shaders: terrain({} uniforms)  atmo({} uniforms)  ocean({} uniforms)",
        tu_count, au_count, ou_count
    );

    let mut hm = Heightmap::new(64);
    hm.generateneptunecloudtopography();
    let h1 = hm.sample(0.0, 0.0);
    let h2 = hm.sample(45.0, 90.0);
    let r1 = hm.radiusat(0.0, 0.0);
    println!(
        "Heightmap(64): h(0,0)={:.0}m  h(45,90)={:.0}m  r(0,0)={:.0}m",
        h1, h2, r1
    );

    let mesh = TerrainMesh::fromregion(0.0, 0.0, 10.0, 10.0, 16, &|lat, lon| hm.sample(lat, lon));
    println!(
        "Mesh: {} vertices  {} triangles",
        mesh.vertexcount(),
        mesh.trianglecount()
    );

    let classifier = CloudBiomeClassifier::default();
    let biome1 = classifier.classify(50000.0, 30.0, 100.0);
    let biome2 = classifier.classify(200000.0, 80.0, 10.0);
    let splat = classifier.splat(80000.0, 45.0, 60.0);
    let wlen = splat.weights.len();
    println!(
        "Biome: 50km/30°={:?}  200km/80°={:?}  splat weights={}",
        biome1, biome2, wlen
    );

    let config = LodConfig::default();
    let mut lod = LodTerrain::new(config);
    lod.update([0.0, 0.0, 1e8]);
    println!("LOD: updated at [0,0,1e8]");

    let ll = LatLon::new(30.0, 60.0, 0.0);
    let ecef = ll.toecef();
    let cart = ll.tocartesiansimple();
    let back = ecef.tolatlon();
    let dist = ll.distanceto(&LatLon::new(0.0, 0.0, 0.0));
    println!(
        "LatLon(30,60): ecef=({:.0},{:.0},{:.0})  cart=({:.0},{:.0},{:.0})",
        ecef.x, ecef.y, ecef.z, cart[0], cart[1], cart[2]
    );
    println!(
        "  roundtrip: ({:.2},{:.2})  dist={:.0}m",
        back.latdeg, back.londeg, dist
    );
    println!(
        "  constants: f={} a={}m b={}m",
        NEPTUNEFLATTENING, NEPTUNESEMIMAJORM, NEPTUNESEMIMINORM
    );

    let bathy = InteriorLayerData::global(3600.0);
    let bval = bathy.sample(0.0, 0.0);
    let isice = bathy.isicemantle(0.0, 0.0);
    let stats = InteriorLayerData::interiorstats();
    println!(
        "Bathymetry: sample(0,0)={:.0}m  isIce={}  atmo={:.0}km  ice={:.0}km  core={:.0}km",
        bval, isice, stats.atmospheredepthkm, stats.icemantledepthkm, stats.coreradiuskm
    );

    let elev = ElevationProvider::global(3600.0);
    let ev = elev.sample(45.0, 90.0);
    let features = ElevationProvider::notablefeatures();
    let fcount = features.len();
    println!(
        "Elevation: sample(45,90)={:.0}m  {} notable features",
        ev, fcount
    );

    let rdb = RegionDatabase::atmosphericbands();
    let regions = rdb.regions;
    let rcount = regions.len();
    let typecounts: Vec<usize> = [
        RegionType::AtmosphericBand,
        RegionType::VortexRegion,
        RegionType::PolarCap,
        RegionType::StormSystem,
        RegionType::CloudFeature,
    ]
    .iter()
    .map(|rt| {
        regions
            .iter()
            .filter(|r| std::mem::discriminant(&r.regiontype) == std::mem::discriminant(rt))
            .count()
    })
    .collect();
    println!(
        "Regions: {} total  bands={}  vortex={}  polar={}  storm={}  cloud={}",
        rcount, typecounts[0], typecounts[1], typecounts[2], typecounts[3], typecounts[4]
    );

    let chemotroph = hypotheticalchemotroph();
    let rate = chemotroph.chemosynthesisrate(100.0, 50.0, 72.0);
    let colprod = chemotroph.columnproduction(rate);
    let vflux = chemotroph.volatilefluxmmday(0.01, 1e-5);
    let npp = chemotroph.nppkgcm2yr(rate);
    let turnover = chemotroph.reservoirturnoveryr(npp);
    let cloudeff = cloudlayerchemosynthesisefficiency(2.0, 0.5);
    println!(
        "Chemotroph: rate={:.6}  colprod={:.6}  vflux={:.6}  npp={:.6}  turnover={:.2}yr",
        rate, colprod, vflux, npp, turnover
    );
    println!("  cloudEfficiency={:.4}", cloudeff);

    let eco = hypotheticalatmosphericbiome();
    let species = eco.totalspecies();
    let individuals = eco.totalindividuals();
    let shannon = eco.shannonindex();
    let simpson = eco.simpsondiversity();
    let expspec = eco.expectedspeciesfromarea(0.25, 10.0);
    let ecnpp = eco.totalnppgcyr();
    let bioturnover = eco.biomassturnovertimeyr(50.0);
    println!(
        "Ecosystem: species={}  individuals={}  H'={:.3}  D={:.3}",
        species, individuals, shannon, simpson
    );
    println!(
        "  expSpecies={:.1}  NPP={:.2}gC/yr  bioTurnover={:.2}yr",
        expspec, ecnpp, bioturnover
    );

    let floater = hypotheticalh2floater();
    let autotr = hypotheticalchemoautotroph();
    let fgrowth = floater.growthrate();
    let fproj = floater.projectforward(1.0);
    let fmetab = floater.metabolicratew();
    let fdouble = floater.doublingtimeyr();
    let adrag = Aeroplankton::dragcoefficient(100.0);
    let aterm = autotr.terminalvelocityms(1e-4, 0.45);
    println!(
        "H2Floater: growth={:.4}  proj1yr={:.1}  metab={:.4e}W  double={:.2}yr",
        fgrowth, fproj, fmetab, fdouble
    );
    println!(
        "Autotroph: drag(Re100)={:.4}  terminal={:.6}m/s",
        adrag, aterm
    );

    let mut prey_pred = CloudLayerInteraction {
        producers: hypotheticalh2floater(),
        consumers: hypotheticalchemoautotroph(),
        interactionrate: 0.001,
        conversionefficiency: 0.1,
        consumerlossrate: 0.05,
    };
    let pg = prey_pred.producergrowthrate();
    let cg = prey_pred.consumergrowthrate();
    prey_pred.step(0.1);
    println!(
        "PredPrey: prodGrowth={:.4}  consGrowth={:.4}  after step: prod={:.1}  cons={:.2}",
        pg, cg, prey_pred.producers.count as f64, prey_pred.consumers.count as f64
    );
}
