/* US76 atmosphere — version 0.1.0
Changelog:
  2026-09-09 — 0.1.0: Establish source version and changelog; numerical code unchanged.
*/
/// US76 atmosphere calculations for geometric altitudes from 0 to 84,852 m.
///
/// Copied from old/lib/algorithms/us76.dart and aligned with
/// webapp/public/us76.js. Constants and their original rounding are retained.
/// Requires only the Dart SDK; no Flutter dependencies.
library;
import 'dart:math';

/// Layered standard-atmosphere model with hydrostatic pressure, ideal-gas
/// density, and Sutherland viscosity. Represents standard dry-air conditions,
/// not measured weather. Sound speed uses a fixed specific-heat ratio of 1.4.
///
/// The retained rounded layer pressures are prescribed independently; small
/// pressure discontinuities at layer boundaries are therefore possible.
class US76Atmosphere {
  /// Inclusive geometric input limit in metres, matching the web calculator.
  /// This is not the upper geopotential boundary of the final layer.
  static const double maxAltitudeM = 84852;

  /// Specific gas constant for dry air, in J/(kg·K).
  static const double R = 287.05; // Specific gas constant for dry air (J/(kg·K))
  static const double g0 = 9.80665; // Standard gravity (m/s²)
  static const double gamma = 1.4; // Ratio of specific heats for dry air
  static const double sutherland_mu0 = 1.716e-5; // Reference dynamic viscosity (Pa·s)
  static const double sutherland_t0 = 273.0; // Viscosity reference temperature (K)
  static const double sutherland_S = 110.4; // Sutherland's constant (K)

  /// Layer bases: geopotential altitude (m), temperature (K), pressure (Pa),
  /// and lapse rate (K/m). These altitudes differ from the geometric input.
  static const List<double> hBase = [0, 11000, 20000, 32000, 47000, 51000, 71000, 84852];
  /// Seven layer-base temperatures in kelvin, indexed by layer.
  static const List<double> tBase = [288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65];
  /// Seven prescribed layer-base pressures in pascals.
  static const List<double> pBase = [101325, 22632.1, 5474.89, 868.02, 110.91, 66.94, 3.96];
  /// Temperature gradients with geopotential altitude, in K/m; zero is isothermal.
  static const List<double> lapseRate = [-0.0065, 0.0, 0.001, 0.0028, 0.0, -0.0028, -0.002];

  /// Computes properties from geometric altitude [h] in metres.
  /// Throws [RangeError] for nonfinite inputs or values outside 0..84,852 m.
  /// Both endpoints are accepted. No unit conversion or display rounding occurs.
  ///
  /// Returned keys:
  /// - `geopotentialAltitude`: converted altitude H = r h / (r + h), in m.
  /// - `barometricAltitude`: geometric input h in m (legacy key name).
  /// - `ambientPressure`: static pressure in Pa.
  /// - `ambientTemperature`: static temperature in K.
  /// - `airDensity`: density in kg/m³.
  /// - `speedOfSound`: sound speed in m/s, using gamma = 1.4.
  /// - `dynamicViscosity`: dynamic viscosity in Pa·s.
  /// - `kinematicViscosity`: dynamic viscosity / density, in m²/s.
  static Map<String, double> getAtmosphericProperties(double h) {
    if (!h.isFinite || h < 0 || h > maxAltitudeM) {
      throw RangeError.value(
        h, 'h', 'Geometric altitude must be between 0 and 84852 m.',
      );
    }
    double geopotentialAltitude = _geopotentialAltitude(h);
    int layer = _getLayer(geopotentialAltitude);
    double temperature = _temperature(geopotentialAltitude, layer);
    double pressure = _pressure(geopotentialAltitude, layer, temperature);
    double density = _density(pressure, temperature);
    double speedOfSound = _speedOfSound(temperature);
    double dynamicViscosity = _dynamicViscosity(temperature);
    double kinematicViscosity = _kinematicViscosity(dynamicViscosity, density);

    return {
      'geopotentialAltitude': geopotentialAltitude,
      'barometricAltitude': h,
      'ambientPressure': pressure,
      'ambientTemperature': temperature,
      'airDensity': density,
      'speedOfSound': speedOfSound,
      'dynamicViscosity': dynamicViscosity,
      'kinematicViscosity': kinematicViscosity,
    };
  }

  /// Converts geometric [h] in m to geopotential altitude in m using the
  /// retained model radius of 6,371,000 m.
  static double _geopotentialAltitude(double h) {
    const double earthRadius = 6371000.0; // Earth radius (m)
    return earthRadius * h / (earthRadius + h);
  }

  /// Selects a layer from geopotential [h] in m. Intervals include the lower
  /// boundary and exclude the upper; an exact boundary uses the higher layer.
  /// Eight boundaries define seven layers. Valid public inputs map into them.
  static int _getLayer(double h) {
    for (int i = 0; i < hBase.length - 1; i++) {
      if (h >= hBase[i] && h < hBase[i + 1]) {
        return i;
      }
    }
    // Validated geometric inputs always map into one of the seven layers.
    throw StateError('No atmospheric layer for geopotential altitude $h m.');
  }

  /// Returns T = Tbase + L (h - hbase) in K for geopotential [h] in m.
  static double _temperature(double h, int layer) {
    double T = tBase[layer];
    double L = lapseRate[layer];
    double h0 = hBase[layer];

    if (L == 0.0) {
      // Isothermal layer
      return T;
    } else {
      // Layer with a temperature gradient
      return T + L * (h - h0);
    }
  }

  /// Returns hydrostatic pressure in Pa at geopotential [h] (m) and [T] (K).
  /// Uses the exponential form for an isothermal layer, otherwise the power law.
  static double _pressure(double h, int layer, double T) {
    double p = pBase[layer];
    double tInit = tBase[layer];
    double L = lapseRate[layer];
    double h0 = hBase[layer];

    if (L == 0.0) {
      // Isothermal layer
      return p * exp(-g0 * (h - h0) / (R * T));
    } else {
      // Layer with a temperature gradient
      return p * pow((T / tInit), (-g0 / (L * R)));
    }
  }

  /// Returns density in kg/m³ from pressure [p] in Pa and temperature [T] in K.
  static double _density(double p, double T) {
    return p / (R * T);
  }

  /// Returns sqrt(gamma R T) in m/s for temperature [T] in K.
  static double _speedOfSound(double T) {
    return sqrt(gamma * R * T);
  }

  /// Returns Sutherland dynamic viscosity in Pa·s for temperature [T] in K.
  static double _dynamicViscosity(double T) {
    return sutherland_mu0 * pow((T / sutherland_t0), 1.5) * (sutherland_t0 + sutherland_S) / (T + sutherland_S);
  }

  /// Returns kinematic viscosity in m²/s from [mu] in Pa·s and [rho] in kg/m³.
  static double _kinematicViscosity(double mu, double rho) {
    return mu / rho;
  }
}