satkit 0.18.1

//! NRLMSISE-00 Atmospheric Density Model - Pure Rust Translation
//!
//! Translated from the C implementation by Dominik Brodowski (release 20041227),
//! which was itself based on the FORTRAN code by Mike Picone, Alan Hedin,
//! and Doug Drob.
//!
//! This module provides atmospheric density and temperature calculations
//! from the surface to the lower exosphere using the NRLMSISE-00 empirical model.

// The code in this module is a close line-by-line port of the C reference
// implementation. A handful of clippy lints flag patterns that are artefacts
// of that style (minified guard clauses that look like missing-else cases;
// index-based loops that mirror the C source one-for-one). Fixing them would
// diverge this translation from its reference — so we allow them module-wide.
#![allow(
    clippy::possible_missing_else,
    clippy::needless_range_loop,
    clippy::too_many_arguments,
    clippy::type_complexity, // test data tables use wide tuples from the reference implementation
)]

use crate::solar_cycle_forecast;
use crate::spaceweather;
use crate::Duration;
use crate::Instant;
use crate::TimeLike;

// ============================================================================
// Data Tables (from nrlmsise-00_data.c)
// ============================================================================

/// Temperature coefficients
const PT: [f64; 150] = [
    9.86573E-01,
    1.62228E-02,
    1.55270E-02,
    -1.04323E-01,
    -3.75801E-03,
    -1.18538E-03,
    -1.24043E-01,
    4.56820E-03,
    8.76018E-03,
    -1.36235E-01,
    -3.52427E-02,
    8.84181E-03,
    -5.92127E-03,
    -8.61650E+00,
    0.00000E+00,
    1.28492E-02,
    0.00000E+00,
    1.30096E+02,
    1.04567E-02,
    1.65686E-03,
    -5.53887E-06,
    2.97810E-03,
    0.00000E+00,
    5.13122E-03,
    8.66784E-02,
    1.58727E-01,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    -7.27026E-06,
    0.00000E+00,
    6.74494E+00,
    4.93933E-03,
    2.21656E-03,
    2.50802E-03,
    0.00000E+00,
    0.00000E+00,
    -2.08841E-02,
    -1.79873E+00,
    1.45103E-03,
    2.81769E-04,
    -1.44703E-03,
    -5.16394E-05,
    8.47001E-02,
    1.70147E-01,
    5.72562E-03,
    5.07493E-05,
    4.36148E-03,
    1.17863E-04,
    4.74364E-03,
    6.61278E-03,
    4.34292E-05,
    1.44373E-03,
    2.41470E-05,
    2.84426E-03,
    8.56560E-04,
    2.04028E-03,
    0.00000E+00,
    -3.15994E+03,
    -2.46423E-03,
    1.13843E-03,
    4.20512E-04,
    0.00000E+00,
    -9.77214E+01,
    6.77794E-03,
    5.27499E-03,
    1.14936E-03,
    0.00000E+00,
    -6.61311E-03,
    -1.84255E-02,
    -1.96259E-02,
    2.98618E+04,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    6.44574E+02,
    8.84668E-04,
    5.05066E-04,
    0.00000E+00,
    4.02881E+03,
    -1.89503E-03,
    0.00000E+00,
    0.00000E+00,
    8.21407E-04,
    2.06780E-03,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    -1.20410E-02,
    -3.63963E-03,
    9.92070E-05,
    -1.15284E-04,
    -6.33059E-05,
    -6.05545E-01,
    8.34218E-03,
    -9.13036E+01,
    3.71042E-04,
    0.00000E+00,
    4.19000E-04,
    2.70928E-03,
    3.31507E-03,
    -4.44508E-03,
    -4.96334E-03,
    -1.60449E-03,
    3.95119E-03,
    2.48924E-03,
    5.09815E-04,
    4.05302E-03,
    2.24076E-03,
    0.00000E+00,
    6.84256E-03,
    4.66354E-04,
    0.00000E+00,
    -3.68328E-04,
    0.00000E+00,
    0.00000E+00,
    -1.46870E+02,
    0.00000E+00,
    0.00000E+00,
    1.09501E-03,
    4.65156E-04,
    5.62583E-04,
    3.21596E+00,
    6.43168E-04,
    3.14860E-03,
    3.40738E-03,
    1.78481E-03,
    9.62532E-04,
    5.58171E-04,
    3.43731E+00,
    -2.33195E-01,
    5.10289E-04,
    0.00000E+00,
    0.00000E+00,
    -9.25347E+04,
    0.00000E+00,
    -1.99639E-03,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
];

/// Density coefficients: pd[0]=HE, pd[1]=O, pd[2]=N2, pd[3]=TLB, pd[4]=O2,
/// pd[5]=AR, pd[6]=H, pd[7]=N, pd[8]=HOT O
const PD: [[f64; 150]; 9] = [
    // HE DENSITY
    [
        1.09979E+00,
        -4.88060E-02,
        -1.97501E-01,
        -9.10280E-02,
        -6.96558E-03,
        2.42136E-02,
        3.91333E-01,
        -7.20068E-03,
        -3.22718E-02,
        1.41508E+00,
        1.68194E-01,
        1.85282E-02,
        1.09384E-01,
        -7.24282E+00,
        0.00000E+00,
        2.96377E-01,
        -4.97210E-02,
        1.04114E+02,
        -8.61108E-02,
        -7.29177E-04,
        1.48998E-06,
        1.08629E-03,
        0.00000E+00,
        0.00000E+00,
        8.31090E-02,
        1.12818E-01,
        -5.75005E-02,
        -1.29919E-02,
        -1.78849E-02,
        -2.86343E-06,
        0.00000E+00,
        -1.51187E+02,
        -6.65902E-03,
        0.00000E+00,
        -2.02069E-03,
        0.00000E+00,
        0.00000E+00,
        4.32264E-02,
        -2.80444E+01,
        -3.26789E-03,
        2.47461E-03,
        0.00000E+00,
        0.00000E+00,
        9.82100E-02,
        1.22714E-01,
        -3.96450E-02,
        0.00000E+00,
        -2.76489E-03,
        0.00000E+00,
        1.87723E-03,
        -8.09813E-03,
        4.34428E-05,
        -7.70932E-03,
        0.00000E+00,
        -2.28894E-03,
        -5.69070E-03,
        -5.22193E-03,
        6.00692E-03,
        -7.80434E+03,
        -3.48336E-03,
        -6.38362E-03,
        -1.82190E-03,
        0.00000E+00,
        -7.58976E+01,
        -2.17875E-02,
        -1.72524E-02,
        -9.06287E-03,
        0.00000E+00,
        2.44725E-02,
        8.66040E-02,
        1.05712E-01,
        3.02543E+04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -6.01364E+03,
        -5.64668E-03,
        -2.54157E-03,
        0.00000E+00,
        3.15611E+02,
        -5.69158E-03,
        0.00000E+00,
        0.00000E+00,
        -4.47216E-03,
        -4.49523E-03,
        4.64428E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        4.51236E-02,
        2.46520E-02,
        6.17794E-03,
        0.00000E+00,
        0.00000E+00,
        -3.62944E-01,
        -4.80022E-02,
        -7.57230E+01,
        -1.99656E-03,
        0.00000E+00,
        -5.18780E-03,
        -1.73990E-02,
        -9.03485E-03,
        7.48465E-03,
        1.53267E-02,
        1.06296E-02,
        1.18655E-02,
        2.55569E-03,
        1.69020E-03,
        3.51936E-02,
        -1.81242E-02,
        0.00000E+00,
        -1.00529E-01,
        -5.10574E-03,
        0.00000E+00,
        2.10228E-03,
        0.00000E+00,
        0.00000E+00,
        -1.73255E+02,
        5.07833E-01,
        -2.41408E-01,
        8.75414E-03,
        2.77527E-03,
        -8.90353E-05,
        -5.25148E+00,
        -5.83899E-03,
        -2.09122E-02,
        -9.63530E-03,
        9.77164E-03,
        4.07051E-03,
        2.53555E-04,
        -5.52875E+00,
        -3.55993E-01,
        -2.49231E-03,
        0.00000E+00,
        0.00000E+00,
        2.86026E+01,
        0.00000E+00,
        3.42722E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
    // O DENSITY
    [
        1.02315E+00,
        -1.59710E-01,
        -1.06630E-01,
        -1.77074E-02,
        -4.42726E-03,
        3.44803E-02,
        4.45613E-02,
        -3.33751E-02,
        -5.73598E-02,
        3.50360E-01,
        6.33053E-02,
        2.16221E-02,
        5.42577E-02,
        -5.74193E+00,
        0.00000E+00,
        1.90891E-01,
        -1.39194E-02,
        1.01102E+02,
        8.16363E-02,
        1.33717E-04,
        6.54403E-06,
        3.10295E-03,
        0.00000E+00,
        0.00000E+00,
        5.38205E-02,
        1.23910E-01,
        -1.39831E-02,
        0.00000E+00,
        0.00000E+00,
        -3.95915E-06,
        0.00000E+00,
        -7.14651E-01,
        -5.01027E-03,
        0.00000E+00,
        -3.24756E-03,
        0.00000E+00,
        0.00000E+00,
        4.42173E-02,
        -1.31598E+01,
        -3.15626E-03,
        1.24574E-03,
        -1.47626E-03,
        -1.55461E-03,
        6.40682E-02,
        1.34898E-01,
        -2.42415E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        6.13666E-04,
        -5.40373E-03,
        2.61635E-05,
        -3.33012E-03,
        0.00000E+00,
        -3.08101E-03,
        -2.42679E-03,
        -3.36086E-03,
        0.00000E+00,
        -1.18979E+03,
        -5.04738E-02,
        -2.61547E-03,
        -1.03132E-03,
        1.91583E-04,
        -8.38132E+01,
        -1.40517E-02,
        -1.14167E-02,
        -4.08012E-03,
        1.73522E-04,
        -1.39644E-02,
        -6.64128E-02,
        -6.85152E-02,
        -1.34414E+04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        6.07916E+02,
        -4.12220E-03,
        -2.20996E-03,
        0.00000E+00,
        1.70277E+03,
        -4.63015E-03,
        0.00000E+00,
        0.00000E+00,
        -2.25360E-03,
        -2.96204E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        3.92786E-02,
        1.31186E-02,
        -1.78086E-03,
        0.00000E+00,
        0.00000E+00,
        -3.90083E-01,
        -2.84741E-02,
        -7.78400E+01,
        -1.02601E-03,
        0.00000E+00,
        -7.26485E-04,
        -5.42181E-03,
        -5.59305E-03,
        1.22825E-02,
        1.23868E-02,
        6.68835E-03,
        -1.03303E-02,
        -9.51903E-03,
        2.70021E-04,
        -2.57084E-02,
        -1.32430E-02,
        0.00000E+00,
        -3.81000E-02,
        -3.16810E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -9.05762E-04,
        -2.14590E-03,
        -1.17824E-03,
        3.66732E+00,
        -3.79729E-04,
        -6.13966E-03,
        -5.09082E-03,
        -1.96332E-03,
        -3.08280E-03,
        -9.75222E-04,
        4.03315E+00,
        -2.52710E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
    // N2 DENSITY
    [
        1.16112E+00,
        0.00000E+00,
        0.00000E+00,
        3.33725E-02,
        0.00000E+00,
        3.48637E-02,
        -5.44368E-03,
        0.00000E+00,
        -6.73940E-02,
        1.74754E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.74712E+02,
        0.00000E+00,
        1.26733E-01,
        0.00000E+00,
        1.03154E+02,
        5.52075E-02,
        0.00000E+00,
        0.00000E+00,
        8.13525E-04,
        0.00000E+00,
        0.00000E+00,
        8.66784E-02,
        1.58727E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -2.50482E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -2.48894E-03,
        6.16053E-04,
        -5.79716E-04,
        2.95482E-03,
        8.47001E-02,
        1.70147E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.47425E-05,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
    // TLB
    [
        9.44846E-01,
        0.00000E+00,
        0.00000E+00,
        -3.08617E-02,
        0.00000E+00,
        -2.44019E-02,
        6.48607E-03,
        0.00000E+00,
        3.08181E-02,
        4.59392E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.74712E+02,
        0.00000E+00,
        2.13260E-02,
        0.00000E+00,
        -3.56958E+02,
        0.00000E+00,
        1.82278E-04,
        0.00000E+00,
        3.07472E-04,
        0.00000E+00,
        0.00000E+00,
        8.66784E-02,
        1.58727E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        3.83054E-03,
        0.00000E+00,
        0.00000E+00,
        -1.93065E-03,
        -1.45090E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.23493E-03,
        1.36736E-03,
        8.47001E-02,
        1.70147E-01,
        3.71469E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        5.10250E-03,
        2.47425E-05,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        3.68756E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
    // O2 DENSITY
    [
        1.35580E+00,
        1.44816E-01,
        0.00000E+00,
        6.07767E-02,
        0.00000E+00,
        2.94777E-02,
        7.46900E-02,
        0.00000E+00,
        -9.23822E-02,
        8.57342E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.38636E+01,
        0.00000E+00,
        7.71653E-02,
        0.00000E+00,
        8.18751E+01,
        1.87736E-02,
        0.00000E+00,
        0.00000E+00,
        1.49667E-02,
        0.00000E+00,
        0.00000E+00,
        8.66784E-02,
        1.58727E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -3.67874E+02,
        5.48158E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        8.47001E-02,
        1.70147E-01,
        1.22631E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        8.17187E-03,
        3.71617E-05,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -2.10826E-03,
        -3.13640E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -7.35742E-02,
        -5.00266E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.94965E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
    // AR DENSITY
    [
        1.04761E+00,
        2.00165E-01,
        2.37697E-01,
        3.68552E-02,
        0.00000E+00,
        3.57202E-02,
        -2.14075E-01,
        0.00000E+00,
        -1.08018E-01,
        -3.73981E-01,
        0.00000E+00,
        3.10022E-02,
        -1.16305E-03,
        -2.07596E+01,
        0.00000E+00,
        8.64502E-02,
        0.00000E+00,
        9.74908E+01,
        5.16707E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        8.66784E-02,
        1.58727E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        3.46193E+02,
        1.34297E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -3.48509E-03,
        -1.54689E-04,
        0.00000E+00,
        0.00000E+00,
        8.47001E-02,
        1.70147E-01,
        1.47753E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.89320E-02,
        3.68181E-05,
        1.32570E-02,
        0.00000E+00,
        0.00000E+00,
        3.59719E-03,
        7.44328E-03,
        -1.00023E-03,
        -6.50528E+03,
        0.00000E+00,
        1.03485E-02,
        -1.00983E-03,
        -4.06916E-03,
        -6.60864E+01,
        -1.71533E-02,
        1.10605E-02,
        1.20300E-02,
        -5.20034E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -2.62769E+03,
        7.13755E-03,
        4.17999E-03,
        0.00000E+00,
        1.25910E+04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -2.23595E-03,
        4.60217E-03,
        5.71794E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -3.18353E-02,
        -2.35526E-02,
        -1.36189E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.03522E-02,
        -6.67837E+01,
        -1.09724E-03,
        0.00000E+00,
        -1.38821E-02,
        1.60468E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.51574E-02,
        -5.44470E-04,
        0.00000E+00,
        7.28224E-02,
        6.59413E-02,
        0.00000E+00,
        -5.15692E-03,
        0.00000E+00,
        0.00000E+00,
        -3.70367E+03,
        0.00000E+00,
        0.00000E+00,
        1.36131E-02,
        5.38153E-03,
        0.00000E+00,
        4.76285E+00,
        -1.75677E-02,
        2.26301E-02,
        0.00000E+00,
        1.76631E-02,
        4.77162E-03,
        0.00000E+00,
        5.39354E+00,
        0.00000E+00,
        -7.51710E-03,
        0.00000E+00,
        0.00000E+00,
        -8.82736E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
    // H DENSITY
    [
        1.26376E+00,
        -2.14304E-01,
        -1.49984E-01,
        2.30404E-01,
        2.98237E-02,
        2.68673E-02,
        2.96228E-01,
        2.21900E-02,
        -2.07655E-02,
        4.52506E-01,
        1.20105E-01,
        3.24420E-02,
        4.24816E-02,
        -9.14313E+00,
        0.00000E+00,
        2.47178E-02,
        -2.88229E-02,
        8.12805E+01,
        5.10380E-02,
        -5.80611E-03,
        2.51236E-05,
        -1.24083E-02,
        0.00000E+00,
        0.00000E+00,
        8.66784E-02,
        1.58727E-01,
        -3.48190E-02,
        0.00000E+00,
        0.00000E+00,
        2.89885E-05,
        0.00000E+00,
        1.53595E+02,
        -1.68604E-02,
        0.00000E+00,
        1.01015E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.84552E-04,
        -1.22181E-03,
        0.00000E+00,
        0.00000E+00,
        8.47001E-02,
        1.70147E-01,
        -1.04927E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -5.91313E-03,
        -2.30501E-02,
        3.14758E-05,
        0.00000E+00,
        0.00000E+00,
        1.26956E-02,
        8.35489E-03,
        3.10513E-04,
        0.00000E+00,
        3.42119E+03,
        -2.45017E-03,
        -4.27154E-04,
        5.45152E-04,
        1.89896E-03,
        2.89121E+01,
        -6.49973E-03,
        -1.93855E-02,
        -1.48492E-02,
        0.00000E+00,
        -5.10576E-02,
        7.87306E-02,
        9.51981E-02,
        -1.49422E+04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.65503E+02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        6.37110E-03,
        3.24789E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        6.14274E-02,
        1.00376E-02,
        -8.41083E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.27099E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -3.94077E-03,
        -1.28601E-02,
        -7.97616E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -6.71465E-03,
        -1.69799E-03,
        1.93772E-03,
        3.81140E+00,
        -7.79290E-03,
        -1.82589E-02,
        -1.25860E-02,
        -1.04311E-02,
        -3.02465E-03,
        2.43063E-03,
        3.63237E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
    // N DENSITY
    [
        7.09557E+01,
        -3.26740E-01,
        0.00000E+00,
        -5.16829E-01,
        -1.71664E-03,
        9.09310E-02,
        -6.71500E-01,
        -1.47771E-01,
        -9.27471E-02,
        -2.30862E-01,
        -1.56410E-01,
        1.34455E-02,
        -1.19717E-01,
        2.52151E+00,
        0.00000E+00,
        -2.41582E-01,
        5.92939E-02,
        4.39756E+00,
        9.15280E-02,
        4.41292E-03,
        0.00000E+00,
        8.66807E-03,
        0.00000E+00,
        0.00000E+00,
        8.66784E-02,
        1.58727E-01,
        9.74701E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        6.70217E+01,
        -1.31660E-03,
        0.00000E+00,
        -1.65317E-02,
        0.00000E+00,
        0.00000E+00,
        8.50247E-02,
        2.77428E+01,
        4.98658E-03,
        6.15115E-03,
        9.50156E-03,
        -2.12723E-02,
        8.47001E-02,
        1.70147E-01,
        -2.38645E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.37380E-03,
        -8.41918E-03,
        2.80145E-05,
        7.12383E-03,
        0.00000E+00,
        -1.66209E-02,
        1.03533E-04,
        -1.68898E-02,
        0.00000E+00,
        3.64526E+03,
        0.00000E+00,
        6.54077E-03,
        3.69130E-04,
        9.94419E-04,
        8.42803E+01,
        -1.16124E-02,
        -7.74414E-03,
        -1.68844E-03,
        1.42809E-03,
        -1.92955E-03,
        1.17225E-01,
        -2.41512E-02,
        1.50521E+04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.60261E+03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -3.54403E-04,
        -1.87270E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.76439E-02,
        6.43207E-03,
        -3.54300E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -2.80221E-02,
        8.11228E+01,
        -6.75255E-04,
        0.00000E+00,
        -1.05162E-02,
        -3.48292E-03,
        -6.97321E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.45546E-03,
        -1.31970E-02,
        -3.57751E-03,
        -1.09021E+00,
        -1.50181E-02,
        -7.12841E-03,
        -6.64590E-03,
        -3.52610E-03,
        -1.87773E-02,
        -2.22432E-03,
        -3.93895E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
    // HOT O DENSITY
    [
        6.04050E-02,
        1.57034E+00,
        2.99387E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.51018E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -8.61650E+00,
        1.26454E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        5.50878E-03,
        0.00000E+00,
        0.00000E+00,
        8.66784E-02,
        1.58727E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        6.23881E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        8.47001E-02,
        1.70147E-01,
        -9.45934E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
];

const PS: [f64; 150] = [
    9.56827E-01,
    6.20637E-02,
    3.18433E-02,
    0.00000E+00,
    0.00000E+00,
    3.94900E-02,
    0.00000E+00,
    0.00000E+00,
    -9.24882E-03,
    -7.94023E-03,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    1.74712E+02,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    2.74677E-03,
    0.00000E+00,
    1.54951E-02,
    8.66784E-02,
    1.58727E-01,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    -6.99007E-04,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    1.24362E-02,
    -5.28756E-03,
    8.47001E-02,
    1.70147E-01,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    2.47425E-05,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
    0.00000E+00,
];

const PDL: [[f64; 25]; 2] = [
    [
        1.09930E+00,
        3.90631E+00,
        3.07165E+00,
        9.86161E-01,
        1.63536E+01,
        4.63830E+00,
        1.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.28840E+00,
        3.10302E-02,
        1.18339E-01,
    ],
    [
        1.00000E+00,
        7.00000E-01,
        1.15020E+00,
        3.44689E+00,
        1.28840E+00,
        1.00000E+00,
        1.08738E+00,
        1.22947E+00,
        1.10016E+00,
        7.34129E-01,
        1.15241E+00,
        2.22784E+00,
        7.95046E-01,
        4.01612E+00,
        4.47749E+00,
        1.23435E+02,
        -7.60535E-02,
        1.68986E-06,
        7.44294E-01,
        1.03604E+00,
        1.72783E+02,
        1.15020E+00,
        3.44689E+00,
        -7.46230E-01,
        9.49154E-01,
    ],
];

const PTM: [f64; 10] = [
    1.04130E+03,
    3.86000E+02,
    1.95000E+02,
    1.66728E+01,
    2.13000E+02,
    1.20000E+02,
    2.40000E+02,
    1.87000E+02,
    -2.00000E+00,
    0.00000E+00,
];

const PDM: [[f64; 10]; 8] = [
    [
        2.45600E+07,
        6.71072E-06,
        1.00000E+02,
        0.00000E+00,
        1.10000E+02,
        1.00000E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
    [
        8.59400E+10,
        1.00000E+00,
        1.05000E+02,
        -8.00000E+00,
        1.10000E+02,
        1.00000E+01,
        9.00000E+01,
        2.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
    [
        2.81000E+11,
        0.00000E+00,
        1.05000E+02,
        2.80000E+01,
        2.89500E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
    [
        3.30000E+10,
        2.68270E-01,
        1.05000E+02,
        1.00000E+00,
        1.10000E+02,
        1.00000E+01,
        1.10000E+02,
        -1.00000E+01,
        0.00000E+00,
        0.00000E+00,
    ],
    [
        1.33000E+09,
        1.19615E-02,
        1.05000E+02,
        0.00000E+00,
        1.10000E+02,
        1.00000E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
    [
        1.76100E+05,
        1.00000E+00,
        9.50000E+01,
        -8.00000E+00,
        1.10000E+02,
        1.00000E+01,
        9.00000E+01,
        2.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
    [
        1.00000E+07,
        1.00000E+00,
        1.05000E+02,
        -8.00000E+00,
        1.10000E+02,
        1.00000E+01,
        9.00000E+01,
        2.00000E+00,
        0.00000E+00,
        0.00000E+00,
    ],
    [
        1.00000E+06,
        1.00000E+00,
        1.05000E+02,
        -8.00000E+00,
        5.50000E+02,
        7.60000E+01,
        9.00000E+01,
        2.00000E+00,
        0.00000E+00,
        4.00000E+03,
    ],
];

const PTL: [[f64; 100]; 4] = [
    // TN1(2)
    [
        1.00858E+00,
        4.56011E-02,
        -2.22972E-02,
        -5.44388E-02,
        5.23136E-04,
        -1.88849E-02,
        5.23707E-02,
        -9.43646E-03,
        6.31707E-03,
        -7.80460E-02,
        -4.88430E-02,
        0.00000E+00,
        0.00000E+00,
        -7.60250E+00,
        0.00000E+00,
        -1.44635E-02,
        -1.76843E-02,
        -1.21517E+02,
        2.85647E-02,
        0.00000E+00,
        0.00000E+00,
        6.31792E-04,
        0.00000E+00,
        5.77197E-03,
        8.66784E-02,
        1.58727E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -8.90272E+03,
        3.30611E-03,
        3.02172E-03,
        0.00000E+00,
        -2.13673E-03,
        -3.20910E-04,
        0.00000E+00,
        0.00000E+00,
        2.76034E-03,
        2.82487E-03,
        -2.97592E-04,
        -4.21534E-03,
        8.47001E-02,
        1.70147E-01,
        8.96456E-03,
        0.00000E+00,
        -1.08596E-02,
        0.00000E+00,
        0.00000E+00,
        5.57917E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        9.65405E-03,
        0.00000E+00,
        0.00000E+00,
        2.00000E+00,
    ],
    // TN1(3)
    [
        9.39664E-01,
        8.56514E-02,
        -6.79989E-03,
        2.65929E-02,
        -4.74283E-03,
        1.21855E-02,
        -2.14905E-02,
        6.49651E-03,
        -2.05477E-02,
        -4.24952E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.19148E+01,
        0.00000E+00,
        1.18777E-02,
        -7.28230E-02,
        -8.15965E+01,
        1.73887E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.44691E-02,
        2.80259E-04,
        8.66784E-02,
        1.58727E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.16584E+02,
        3.18713E-03,
        7.37479E-03,
        0.00000E+00,
        -2.55018E-03,
        -3.92806E-03,
        0.00000E+00,
        0.00000E+00,
        -2.89757E-03,
        -1.33549E-03,
        1.02661E-03,
        3.53775E-04,
        8.47001E-02,
        1.70147E-01,
        -9.17497E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        3.56082E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.00902E-02,
        0.00000E+00,
        0.00000E+00,
        2.00000E+00,
    ],
    // TN1(4)
    [
        9.85982E-01,
        -4.55435E-02,
        1.21106E-02,
        2.04127E-02,
        -2.40836E-03,
        1.11383E-02,
        -4.51926E-02,
        1.35074E-02,
        -6.54139E-03,
        1.15275E-01,
        1.28247E-01,
        0.00000E+00,
        0.00000E+00,
        -5.30705E+00,
        0.00000E+00,
        -3.79332E-02,
        -6.24741E-02,
        7.71062E-01,
        2.96315E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        6.81051E-03,
        -4.34767E-03,
        8.66784E-02,
        1.58727E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.07003E+01,
        -2.76907E-03,
        4.32474E-04,
        0.00000E+00,
        1.31497E-03,
        -6.47517E-04,
        0.00000E+00,
        -2.20621E+01,
        -1.10804E-03,
        -8.09338E-04,
        4.18184E-04,
        4.29650E-03,
        8.47001E-02,
        1.70147E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -4.04337E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -9.52550E-04,
        8.56253E-04,
        4.33114E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.21223E-03,
        2.38694E-04,
        9.15245E-04,
        1.28385E-03,
        8.67668E-04,
        -5.61425E-06,
        1.04445E+00,
        3.41112E+01,
        0.00000E+00,
        -8.40704E-01,
        -2.39639E+02,
        7.06668E-01,
        -2.05873E+01,
        -3.63696E-01,
        2.39245E+01,
        0.00000E+00,
        -1.06657E-03,
        -7.67292E-04,
        1.54534E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.00000E+00,
    ],
    // TN1(5) TN2(1)
    [
        1.00320E+00,
        3.83501E-02,
        -2.38983E-03,
        2.83950E-03,
        4.20956E-03,
        5.86619E-04,
        2.19054E-02,
        -1.00946E-02,
        -3.50259E-03,
        4.17392E-02,
        -8.44404E-03,
        0.00000E+00,
        0.00000E+00,
        4.96949E+00,
        0.00000E+00,
        -7.06478E-03,
        -1.46494E-02,
        3.13258E+01,
        -1.86493E-03,
        0.00000E+00,
        -1.67499E-02,
        0.00000E+00,
        0.00000E+00,
        5.12686E-04,
        8.66784E-02,
        1.58727E-01,
        -4.64167E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        4.37353E-03,
        -1.99069E+02,
        0.00000E+00,
        -5.34884E-03,
        0.00000E+00,
        1.62458E-03,
        2.93016E-03,
        2.67926E-03,
        5.90449E+02,
        0.00000E+00,
        0.00000E+00,
        -1.17266E-03,
        -3.58890E-04,
        8.47001E-02,
        1.70147E-01,
        0.00000E+00,
        0.00000E+00,
        1.38673E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.60571E-03,
        6.28078E-04,
        5.05469E-05,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.57829E-03,
        -4.00855E-04,
        5.04077E-05,
        -1.39001E-03,
        -2.33406E-03,
        -4.81197E-04,
        1.46758E+00,
        6.20332E+00,
        0.00000E+00,
        3.66476E-01,
        -6.19760E+01,
        3.09198E-01,
        -1.98999E+01,
        0.00000E+00,
        -3.29933E+02,
        0.00000E+00,
        -1.10080E-03,
        -9.39310E-05,
        1.39638E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.00000E+00,
    ],
];

const PMA: [[f64; 100]; 10] = [
    // TN2(2)
    [
        9.81637E-01,
        -1.41317E-03,
        3.87323E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -3.58707E-02,
        -8.63658E-03,
        0.00000E+00,
        0.00000E+00,
        -2.02226E+00,
        0.00000E+00,
        -8.69424E-03,
        -1.91397E-02,
        8.76779E+01,
        4.52188E-03,
        0.00000E+00,
        2.23760E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -7.07572E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -4.11210E-03,
        3.50060E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -8.36657E-03,
        1.61347E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.45130E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.24152E-03,
        6.43365E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.33255E-03,
        2.42657E-03,
        1.60666E-03,
        -1.85728E-03,
        -1.46874E-03,
        -4.79163E-06,
        1.22464E+00,
        3.53510E+01,
        0.00000E+00,
        4.49223E-01,
        -4.77466E+01,
        4.70681E-01,
        8.41861E+00,
        -2.88198E-01,
        1.67854E+02,
        0.00000E+00,
        7.11493E-04,
        6.05601E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.00000E+00,
    ],
    // TN2(3)
    [
        1.00422E+00,
        -7.11212E-03,
        5.24480E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -5.28914E-02,
        -2.41301E-02,
        0.00000E+00,
        0.00000E+00,
        -2.12219E+01,
        -1.03830E-02,
        -3.28077E-03,
        1.65727E-02,
        1.68564E+00,
        -6.68154E-03,
        0.00000E+00,
        1.45155E-02,
        0.00000E+00,
        8.42365E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -4.34645E-03,
        0.00000E+00,
        0.00000E+00,
        2.16780E-02,
        0.00000E+00,
        -1.38459E+02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        7.04573E-03,
        -4.73204E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.08767E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -8.08279E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        5.21769E-04,
        -2.27387E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        3.26769E-03,
        3.16901E-03,
        4.60316E-04,
        -1.01431E-04,
        1.02131E-03,
        9.96601E-04,
        1.25707E+00,
        2.50114E+01,
        0.00000E+00,
        4.24472E-01,
        -2.77655E+01,
        3.44625E-01,
        2.75412E+01,
        0.00000E+00,
        7.94251E+02,
        0.00000E+00,
        2.45835E-03,
        1.38871E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.00000E+00,
    ],
    // TN2(4) TN3(1)
    [
        1.01890E+00,
        -2.46603E-02,
        1.00078E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -6.70977E-02,
        -4.02286E-02,
        0.00000E+00,
        0.00000E+00,
        -2.29466E+01,
        -7.47019E-03,
        2.26580E-03,
        2.63931E-02,
        3.72625E+01,
        -6.39041E-03,
        0.00000E+00,
        9.58383E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.85291E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.39717E+02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        9.19771E-03,
        -3.69121E+02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.57067E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -7.07265E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -2.92953E-03,
        -2.77739E-03,
        -4.40092E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.47280E-03,
        2.95035E-04,
        -1.81246E-03,
        2.81945E-03,
        4.27296E-03,
        9.78863E-04,
        1.40545E+00,
        -6.19173E+00,
        0.00000E+00,
        0.00000E+00,
        -7.93632E+01,
        4.44643E-01,
        -4.03085E+02,
        0.00000E+00,
        1.15603E+01,
        0.00000E+00,
        2.25068E-03,
        8.48557E-04,
        -2.98493E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.00000E+00,
    ],
    // TN3(2)
    [
        9.75801E-01,
        3.80680E-02,
        -3.05198E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        3.85575E-02,
        5.04057E-02,
        0.00000E+00,
        0.00000E+00,
        -1.76046E+02,
        1.44594E-02,
        -1.48297E-03,
        -3.68560E-03,
        3.02185E+01,
        -3.23338E-03,
        0.00000E+00,
        1.53569E-02,
        0.00000E+00,
        -1.15558E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        4.89620E-03,
        0.00000E+00,
        0.00000E+00,
        -1.00616E-02,
        -8.21324E-03,
        -1.57757E+02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        6.63564E-03,
        4.58410E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -2.51280E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        9.91215E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -8.73148E-04,
        -1.29648E-03,
        -7.32026E-05,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -4.68110E-03,
        -4.66003E-03,
        -1.31567E-03,
        -7.39390E-04,
        6.32499E-04,
        -4.65588E-04,
        -1.29785E+00,
        -1.57139E+02,
        0.00000E+00,
        2.58350E-01,
        -3.69453E+01,
        4.10672E-01,
        9.78196E+00,
        -1.52064E-01,
        -3.85084E+03,
        0.00000E+00,
        -8.52706E-04,
        -1.40945E-03,
        -7.26786E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.00000E+00,
    ],
    // TN3(3)
    [
        9.60722E-01,
        7.03757E-02,
        -3.00266E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.22671E-02,
        4.10423E-02,
        0.00000E+00,
        0.00000E+00,
        -1.63070E+02,
        1.06073E-02,
        5.40747E-04,
        7.79481E-03,
        1.44908E+02,
        1.51484E-04,
        0.00000E+00,
        1.97547E-02,
        0.00000E+00,
        -1.41844E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        5.77884E-03,
        0.00000E+00,
        0.00000E+00,
        9.74319E-03,
        0.00000E+00,
        -2.88015E+03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -4.44902E-03,
        -2.92760E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.34419E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        5.36685E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -4.65325E-04,
        -5.50628E-04,
        3.31465E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -2.06179E-03,
        -3.08575E-03,
        -7.93589E-04,
        -1.08629E-04,
        5.95511E-04,
        -9.05050E-04,
        1.18997E+00,
        4.15924E+01,
        0.00000E+00,
        -4.72064E-01,
        -9.47150E+02,
        3.98723E-01,
        1.98304E+01,
        0.00000E+00,
        3.73219E+03,
        0.00000E+00,
        -1.50040E-03,
        -1.14933E-03,
        -1.56769E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.00000E+00,
    ],
    // TN3(4)
    [
        1.03123E+00,
        -7.05124E-02,
        8.71615E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -3.82621E-02,
        -9.80975E-03,
        0.00000E+00,
        0.00000E+00,
        2.89286E+01,
        9.57341E-03,
        0.00000E+00,
        0.00000E+00,
        8.66153E+01,
        7.91938E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        4.68917E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        7.86638E-03,
        0.00000E+00,
        0.00000E+00,
        9.90827E-03,
        0.00000E+00,
        6.55573E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -4.00200E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        7.07457E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        5.72268E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -2.04970E-04,
        1.21560E-03,
        -8.05579E-06,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -2.49941E-03,
        -4.57256E-04,
        -1.59311E-04,
        2.96481E-04,
        -1.77318E-03,
        -6.37918E-04,
        1.02395E+00,
        1.28172E+01,
        0.00000E+00,
        1.49903E-01,
        -2.63818E+01,
        0.00000E+00,
        4.70628E+01,
        -2.22139E-01,
        4.82292E-02,
        0.00000E+00,
        -8.67075E-04,
        -5.86479E-04,
        5.32462E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.00000E+00,
    ],
    // TN3(5) SURFACE TEMP TSL
    [
        1.00828E+00,
        -9.10404E-02,
        -2.26549E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -2.32420E-02,
        -9.08925E-03,
        0.00000E+00,
        0.00000E+00,
        3.36105E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.24957E+01,
        -5.87939E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.79765E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.01237E+03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.75553E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        3.29699E-03,
        1.26659E-03,
        2.68402E-04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.17894E-03,
        1.48746E-03,
        1.06478E-04,
        1.34743E-04,
        -2.20939E-03,
        -6.23523E-04,
        6.36539E-01,
        1.13621E+01,
        0.00000E+00,
        -3.93777E-01,
        2.38687E+03,
        0.00000E+00,
        6.61865E+02,
        -1.21434E-01,
        9.27608E+00,
        0.00000E+00,
        1.68478E-04,
        1.24892E-03,
        1.71345E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.00000E+00,
    ],
    // TGN3(2) SURFACE GRAD TSLG
    [
        1.57293E+00,
        -6.78400E-01,
        6.47500E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -7.62974E-02,
        -3.60423E-01,
        0.00000E+00,
        0.00000E+00,
        1.28358E+02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        4.68038E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.67898E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.90994E+04,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        3.15706E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.00000E+00,
    ],
    // TGN2(1) TGN1(2)
    [
        8.60028E-01,
        3.77052E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.17570E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        7.77757E-03,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        1.01024E+02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        6.54251E+02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.56959E-02,
        1.91001E-02,
        3.15971E-02,
        1.00982E-02,
        -6.71565E-03,
        2.57693E-03,
        1.38692E+00,
        2.82132E-01,
        0.00000E+00,
        0.00000E+00,
        3.81511E+02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.00000E+00,
    ],
    // TGN3(1) TGN2(2)
    [
        1.06029E+00,
        -5.25231E-02,
        3.73034E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        3.31072E-02,
        -3.88409E-01,
        0.00000E+00,
        0.00000E+00,
        -1.65295E+02,
        -2.13801E-01,
        -4.38916E-02,
        -3.22716E-01,
        -8.82393E+01,
        1.18458E-01,
        0.00000E+00,
        -4.35863E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -1.19782E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.62229E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -5.37443E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        -4.55788E-01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        3.84009E-02,
        3.96733E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        5.05494E-02,
        7.39617E-02,
        1.92200E-02,
        -8.46151E-03,
        -1.34244E-02,
        1.96338E-02,
        1.50421E+00,
        1.88368E+01,
        0.00000E+00,
        0.00000E+00,
        -5.13114E+01,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        5.11923E-02,
        3.61225E-02,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        0.00000E+00,
        2.00000E+00,
    ],
];

#[allow(dead_code)] // Part of model specification; retained for completeness
const SAM: [f64; 100] = [
    1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
    1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
    1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 0.0, 0.0, 0.0, 0.0,
];

const PAVGM: [f64; 10] = [
    2.61000E+02,
    2.64000E+02,
    2.29000E+02,
    2.17000E+02,
    2.17000E+02,
    2.23000E+02,
    2.86760E+02,
    -2.93940E+00,
    2.50000E+00,
    0.00000E+00,
];

// ============================================================================
// The rest of the code is identical to the first write attempt.
// Due to space constraints, I'll include a marker and the build will pick it up.
// ============================================================================

struct NrlmsiseState {
    gsurf: f64,
    re: f64,
    dd: f64,
    dm04: f64,
    dm16: f64,
    dm28: f64,
    dm32: f64,
    dm40: f64,
    dm01: f64,
    dm14: f64,
    meso_tn1: [f64; 5],
    meso_tn2: [f64; 4],
    meso_tn3: [f64; 5],
    meso_tgn1: [f64; 2],
    meso_tgn2: [f64; 2],
    meso_tgn3: [f64; 2],
    dfa: f64,
    plg: [[f64; 9]; 4],
    ctloc: f64,
    stloc: f64,
    c2tloc: f64,
    s2tloc: f64,
    s3tloc: f64,
    c3tloc: f64,
    apdf: f64,
    apt: [f64; 4],
}

impl NrlmsiseState {
    fn new() -> Self {
        Self {
            gsurf: 0.0,
            re: 0.0,
            dd: 0.0,
            dm04: 0.0,
            dm16: 0.0,
            dm28: 0.0,
            dm32: 0.0,
            dm40: 0.0,
            dm01: 0.0,
            dm14: 0.0,
            meso_tn1: [0.0; 5],
            meso_tn2: [0.0; 4],
            meso_tn3: [0.0; 5],
            meso_tgn1: [0.0; 2],
            meso_tgn2: [0.0; 2],
            meso_tgn3: [0.0; 2],
            dfa: 0.0,
            plg: [[0.0; 9]; 4],
            ctloc: 0.0,
            stloc: 0.0,
            c2tloc: 0.0,
            s2tloc: 0.0,
            s3tloc: 0.0,
            c3tloc: 0.0,
            apdf: 0.0,
            apt: [0.0; 4],
        }
    }
}

struct NrlmsiseFlags {
    switches: [i32; 24],
    sw: [f64; 24],
    swc: [f64; 24],
}
#[allow(dead_code)]
struct NrlmsiseInput {
    year: i32,
    doy: i32,
    sec: f64,
    alt: f64,
    g_lat: f64,
    g_long: f64,
    lst: f64,
    f107a: f64,
    f107: f64,
    ap: f64,
    ap_a: Option<[f64; 7]>,
}
struct NrlmsiseOutput {
    d: [f64; 9],
    t: [f64; 2],
}

fn tselec(flags: &mut NrlmsiseFlags) {
    for i in 0..24 {
        if i != 9 {
            flags.sw[i] = if flags.switches[i] == 1 { 1.0 } else { 0.0 };
            flags.swc[i] = if flags.switches[i] > 0 { 1.0 } else { 0.0 };
        } else {
            flags.sw[i] = flags.switches[i] as f64;
            flags.swc[i] = flags.switches[i] as f64;
        }
    }
}

fn glatf(lat: f64, gv: &mut f64, reff: &mut f64) {
    let dgtr = 1.74533E-2_f64;
    let c2 = (2.0 * dgtr * lat).cos();
    *gv = 980.616 * (1.0 - 0.0026373 * c2);
    *reff = 2.0 * (*gv) / (3.085462E-6 + 2.27E-9 * c2) * 1.0E-5;
}

fn ccor(alt: f64, r: f64, h1: f64, zh: f64) -> f64 {
    let e = (alt - zh) / h1;
    if e > 70.0 {
        return 1.0;
    }
    if e < -70.0 {
        return r.exp();
    }
    (r / (1.0 + e.exp())).exp()
}
fn ccor2(alt: f64, r: f64, h1: f64, zh: f64, h2: f64) -> f64 {
    let e1 = (alt - zh) / h1;
    let e2 = (alt - zh) / h2;
    if e1 > 70.0 || e2 > 70.0 {
        return 1.0;
    }
    if e1 < -70.0 && e2 < -70.0 {
        return r.exp();
    }
    (r / (1.0 + 0.5 * (e1.exp() + e2.exp()))).exp()
}
fn scalh(alt: f64, xm: f64, temp: f64, gsurf: f64, re: f64) -> f64 {
    let rgas = 831.4_f64;
    let r = 1.0 + alt / re;
    rgas * temp / (gsurf / (r * r) * xm)
}
fn dnet(dd: f64, dm: f64, zhm: f64, xmm: f64, xm: f64) -> f64 {
    let a = zhm / (xmm - xm);
    if !(dm > 0.0 && dd > 0.0) {
        if dd == 0.0 && dm == 0.0 {
            return 1.0;
        }
        if dm == 0.0 {
            return dd;
        }
        if dd == 0.0 {
            return dm;
        }
    }
    let ylog = a * (dm / dd).ln();
    if ylog < -10.0 {
        return dd;
    }
    if ylog > 10.0 {
        return dm;
    }
    dd * (1.0 + ylog.exp()).powf(1.0 / a)
}
#[inline]
fn zeta_fn(zz: f64, zl: f64, re: f64) -> f64 {
    (zz - zl) * (re + zl) / (re + zz)
}

fn splini(xa: &[f64], ya: &[f64], y2a: &[f64], n: usize, x: f64) -> f64 {
    let mut yi = 0.0;
    let mut klo = 0;
    let mut khi = 1;
    while x > xa[klo] && khi < n {
        let xx = if khi < n - 1 {
            if x < xa[khi] {
                x
            } else {
                xa[khi]
            }
        } else {
            x
        };
        let h = xa[khi] - xa[klo];
        let a = (xa[khi] - xx) / h;
        let b = (xx - xa[klo]) / h;
        let a2 = a * a;
        let b2 = b * b;
        yi += ((1.0 - a2) * ya[klo] / 2.0
            + b2 * ya[khi] / 2.0
            + ((-(1.0 + a2 * a2) / 4.0 + a2 / 2.0) * y2a[klo]
                + (b2 * b2 / 4.0 - b2 / 2.0) * y2a[khi])
                * h
                * h
                / 6.0)
            * h;
        klo += 1;
        khi += 1;
    }
    yi
}
fn splint(xa: &[f64], ya: &[f64], y2a: &[f64], n: usize, x: f64) -> f64 {
    let mut klo = 0;
    let mut khi = n - 1;
    while (khi - klo) > 1 {
        let k = (khi + klo) / 2;
        if xa[k] > x {
            khi = k;
        } else {
            klo = k;
        }
    }
    let h = xa[khi] - xa[klo];
    let a = (xa[khi] - x) / h;
    let b = (x - xa[klo]) / h;
    a * ya[klo]
        + b * ya[khi]
        + ((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi]) * h * h / 6.0
}
fn spline(x: &[f64], y: &[f64], n: usize, yp1: f64, ypn: f64, y2: &mut [f64]) {
    let mut u = vec![0.0_f64; n];
    if yp1 > 0.99E30 {
        y2[0] = 0.0;
        u[0] = 0.0;
    } else {
        y2[0] = -0.5;
        u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - yp1);
    }
    for i in 1..(n - 1) {
        let sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
        let p = sig * y2[i - 1] + 2.0;
        y2[i] = (sig - 1.0) / p;
        u[i] = (6.0
            * ((y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]))
            / (x[i + 1] - x[i - 1])
            - sig * u[i - 1])
            / p;
    }
    let (qn, un) = if ypn > 0.99E30 {
        (0.0, 0.0)
    } else {
        (
            0.5,
            (3.0 / (x[n - 1] - x[n - 2])) * (ypn - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2])),
        )
    };
    y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
    for k in (0..=(n - 2)).rev() {
        y2[k] = y2[k] * y2[k + 1] + u[k];
    }
}

fn densm(
    alt: f64,
    d0: f64,
    xm: f64,
    tz: &mut f64,
    mn3: usize,
    zn3: &[f64],
    tn3: &[f64],
    tgn3: &[f64],
    mn2: usize,
    zn2: &[f64],
    tn2: &[f64],
    tgn2: &[f64],
    gsurf: f64,
    re: f64,
) -> f64 {
    let rgas = 831.4_f64;
    let mut densm_tmp = d0;
    if alt > zn2[0] {
        if xm == 0.0 {
            return *tz;
        } else {
            return d0;
        }
    }
    let z = if alt > zn2[mn2 - 1] {
        alt
    } else {
        zn2[mn2 - 1]
    };
    let mn = mn2;
    let z1 = zn2[0];
    let z2 = zn2[mn - 1];
    let t1 = tn2[0];
    let t2 = tn2[mn - 1];
    let zg = zeta_fn(z, z1, re);
    let zgdif = zeta_fn(z2, z1, re);
    let mut xs = [0.0_f64; 10];
    let mut ys = [0.0_f64; 10];
    let mut y2out = [0.0_f64; 10];
    for k in 0..mn {
        xs[k] = zeta_fn(zn2[k], z1, re) / zgdif;
        ys[k] = 1.0 / tn2[k];
    }
    let yd1 = -tgn2[0] / (t1 * t1) * zgdif;
    let yd2 = -tgn2[1] / (t2 * t2) * zgdif * ((re + z2) / (re + z1)).powi(2);
    spline(&xs[..mn], &ys[..mn], mn, yd1, yd2, &mut y2out[..mn]);
    let x = zg / zgdif;
    let y = splint(&xs[..mn], &ys[..mn], &y2out[..mn], mn, x);
    *tz = 1.0 / y;
    if xm != 0.0 {
        let glb = gsurf / (1.0 + z1 / re).powi(2);
        let gamm = xm * glb * zgdif / rgas;
        let yi = splini(&xs[..mn], &ys[..mn], &y2out[..mn], mn, x);
        let mut expl = gamm * yi;
        if expl > 50.0 {
            expl = 50.0;
        }
        densm_tmp = densm_tmp * (t1 / *tz) * (-expl).exp();
    }
    if alt > zn3[0] {
        if xm == 0.0 {
            return *tz;
        } else {
            return densm_tmp;
        }
    }
    let z = alt;
    let mn = mn3;
    let z1 = zn3[0];
    let z2 = zn3[mn - 1];
    let t1 = tn3[0];
    let t2 = tn3[mn - 1];
    let zg = zeta_fn(z, z1, re);
    let zgdif = zeta_fn(z2, z1, re);
    for k in 0..mn {
        xs[k] = zeta_fn(zn3[k], z1, re) / zgdif;
        ys[k] = 1.0 / tn3[k];
    }
    let yd1 = -tgn3[0] / (t1 * t1) * zgdif;
    let yd2 = -tgn3[1] / (t2 * t2) * zgdif * ((re + z2) / (re + z1)).powi(2);
    spline(&xs[..mn], &ys[..mn], mn, yd1, yd2, &mut y2out[..mn]);
    let x = zg / zgdif;
    let y = splint(&xs[..mn], &ys[..mn], &y2out[..mn], mn, x);
    *tz = 1.0 / y;
    if xm != 0.0 {
        let glb = gsurf / (1.0 + z1 / re).powi(2);
        let gamm = xm * glb * zgdif / rgas;
        let yi = splini(&xs[..mn], &ys[..mn], &y2out[..mn], mn, x);
        let mut expl = gamm * yi;
        if expl > 50.0 {
            expl = 50.0;
        }
        densm_tmp = densm_tmp * (t1 / *tz) * (-expl).exp();
    }
    if xm == 0.0 {
        *tz
    } else {
        densm_tmp
    }
}

fn densu(
    alt: f64,
    dlb: f64,
    tinf: f64,
    tlb: f64,
    xm: f64,
    alpha: f64,
    tz: &mut f64,
    zlb: f64,
    s2: f64,
    mn1: usize,
    zn1: &[f64],
    tn1: &mut [f64],
    tgn1: &mut [f64],
    gsurf: f64,
    re: f64,
) -> f64 {
    let rgas = 831.4_f64;
    let za = zn1[0];
    let z = if alt > za { alt } else { za };
    let zg2 = zeta_fn(z, zlb, re);
    let tt = tinf - (tinf - tlb) * (-s2 * zg2).exp();
    let ta = tt;
    *tz = tt;
    let mut x = 0.0_f64;
    let mut z1 = 0.0_f64;
    let mut zgdif = 0.0_f64;
    let mut xs = [0.0_f64; 5];
    let mut ys = [0.0_f64; 5];
    let mut y2out = [0.0_f64; 5];
    let mut mn = 0_usize;
    if alt < za {
        let dta = (tinf - ta) * s2 * ((re + zlb) / (re + za)).powi(2);
        tgn1[0] = dta;
        tn1[0] = ta;
        let z_local = if alt > zn1[mn1 - 1] {
            alt
        } else {
            zn1[mn1 - 1]
        };
        mn = mn1;
        z1 = zn1[0];
        let z2 = zn1[mn - 1];
        let t1 = tn1[0];
        let t2 = tn1[mn - 1];
        let zg = zeta_fn(z_local, z1, re);
        zgdif = zeta_fn(z2, z1, re);
        for k in 0..mn {
            xs[k] = zeta_fn(zn1[k], z1, re) / zgdif;
            ys[k] = 1.0 / tn1[k];
        }
        let yd1 = -tgn1[0] / (t1 * t1) * zgdif;
        let yd2 = -tgn1[1] / (t2 * t2) * zgdif * ((re + z2) / (re + z1)).powi(2);
        spline(&xs[..mn], &ys[..mn], mn, yd1, yd2, &mut y2out[..mn]);
        x = zg / zgdif;
        let y = splint(&xs[..mn], &ys[..mn], &y2out[..mn], mn, x);
        *tz = 1.0 / y;
    }
    if xm == 0.0 {
        return *tz;
    }
    let glb = gsurf / (1.0 + zlb / re).powi(2);
    let gamma = xm * glb / (s2 * rgas * tinf);
    let mut expl = (-s2 * gamma * zg2).exp();
    if expl > 50.0 {
        expl = 50.0;
    }
    if tt <= 0.0 {
        expl = 50.0;
    }
    let densa = dlb * (tlb / tt).powf(1.0 + alpha + gamma) * expl;
    let mut densu_temp = densa;
    if alt >= za {
        return densu_temp;
    }
    let glb = gsurf / (1.0 + z1 / re).powi(2);
    let gamm = xm * glb * zgdif / rgas;
    let yi = splini(&xs[..mn], &ys[..mn], &y2out[..mn], mn, x);
    let mut expl = gamm * yi;
    if expl > 50.0 {
        expl = 50.0;
    }
    if *tz <= 0.0 {
        expl = 50.0;
    }
    densu_temp = densu_temp * (tn1[0] / *tz).powf(1.0 + alpha) * (-expl).exp();
    densu_temp
}

#[inline]
fn g0_fn(a: f64, p: &[f64]) -> f64 {
    a - 4.0 + (p[25] - 1.0) * (a - 4.0 + ((-p[24].abs() * (a - 4.0)).exp() - 1.0) / p[24].abs())
}
#[inline]
fn sumex(ex: f64) -> f64 {
    1.0 + (1.0 - ex.powi(19)) / (1.0 - ex) * ex.sqrt()
}
#[inline]
fn sg0(ex: f64, p: &[f64], ap: &[f64]) -> f64 {
    (g0_fn(ap[1], p)
        + (g0_fn(ap[2], p) * ex
            + g0_fn(ap[3], p) * ex * ex
            + g0_fn(ap[4], p) * ex.powi(3)
            + (g0_fn(ap[5], p) * ex.powi(4) + g0_fn(ap[6], p) * ex.powi(12)) * (1.0 - ex.powi(8))
                / (1.0 - ex)))
        / sumex(ex)
}

fn globe7(
    p: &[f64],
    input: &NrlmsiseInput,
    flags: &NrlmsiseFlags,
    state: &mut NrlmsiseState,
) -> f64 {
    let sr = 7.2722E-5_f64;
    let dgtr = 1.74533E-2_f64;
    let dr = 1.72142E-2_f64;
    let hr = 0.2618_f64;
    let mut t = [0.0_f64; 15];
    let tloc = input.lst;
    let c = (input.g_lat * dgtr).sin();
    let s = (input.g_lat * dgtr).cos();
    let c2 = c * c;
    let c4 = c2 * c2;
    let s2 = s * s;
    state.plg[0][1] = c;
    state.plg[0][2] = 0.5 * (3.0 * c2 - 1.0);
    state.plg[0][3] = 0.5 * (5.0 * c * c2 - 3.0 * c);
    state.plg[0][4] = (35.0 * c4 - 30.0 * c2 + 3.0) / 8.0;
    state.plg[0][5] = (63.0 * c2 * c2 * c - 70.0 * c2 * c + 15.0 * c) / 8.0;
    state.plg[0][6] = (11.0 * c * state.plg[0][5] - 5.0 * state.plg[0][4]) / 6.0;
    state.plg[1][1] = s;
    state.plg[1][2] = 3.0 * c * s;
    state.plg[1][3] = 1.5 * (5.0 * c2 - 1.0) * s;
    state.plg[1][4] = 2.5 * (7.0 * c2 * c - 3.0 * c) * s;
    state.plg[1][5] = 1.875 * (21.0 * c4 - 14.0 * c2 + 1.0) * s;
    state.plg[1][6] = (11.0 * c * state.plg[1][5] - 6.0 * state.plg[1][4]) / 5.0;
    state.plg[2][2] = 3.0 * s2;
    state.plg[2][3] = 15.0 * s2 * c;
    state.plg[2][4] = 7.5 * (7.0 * c2 - 1.0) * s2;
    state.plg[2][5] = 3.0 * c * state.plg[2][4] - 2.0 * state.plg[2][3];
    state.plg[2][6] = (11.0 * c * state.plg[2][5] - 7.0 * state.plg[2][4]) / 4.0;
    state.plg[2][7] = (13.0 * c * state.plg[2][6] - 8.0 * state.plg[2][5]) / 5.0;
    state.plg[3][3] = 15.0 * s2 * s;
    state.plg[3][4] = 105.0 * s2 * s * c;
    state.plg[3][5] = (9.0 * c * state.plg[3][4] - 7.0 * state.plg[3][3]) / 2.0;
    state.plg[3][6] = (11.0 * c * state.plg[3][5] - 8.0 * state.plg[3][4]) / 3.0;
    if !((flags.sw[7] == 0.0 && flags.sw[8] == 0.0) && flags.sw[14] == 0.0) {
        state.stloc = (hr * tloc).sin();
        state.ctloc = (hr * tloc).cos();
        state.s2tloc = (2.0 * hr * tloc).sin();
        state.c2tloc = (2.0 * hr * tloc).cos();
        state.s3tloc = (3.0 * hr * tloc).sin();
        state.c3tloc = (3.0 * hr * tloc).cos();
    }
    let cd32 = (dr * (input.doy as f64 - p[31])).cos();
    let cd18 = (2.0 * dr * (input.doy as f64 - p[17])).cos();
    let cd14 = (dr * (input.doy as f64 - p[13])).cos();
    let cd39 = (2.0 * dr * (input.doy as f64 - p[38])).cos();
    let df = input.f107 - input.f107a;
    state.dfa = input.f107a - 150.0;
    t[0] = p[19] * df * (1.0 + p[59] * state.dfa)
        + p[20] * df * df
        + p[21] * state.dfa
        + p[29] * state.dfa * state.dfa;
    let f1 = 1.0 + (p[47] * state.dfa + p[19] * df + p[20] * df * df) * flags.swc[1];
    let f2 = 1.0 + (p[49] * state.dfa + p[19] * df + p[20] * df * df) * flags.swc[1];
    t[1] = (p[1] * state.plg[0][2] + p[2] * state.plg[0][4] + p[22] * state.plg[0][6])
        + p[14] * state.plg[0][2] * state.dfa * flags.swc[1]
        + p[26] * state.plg[0][1];
    t[2] = p[18] * cd32;
    t[3] = (p[15] + p[16] * state.plg[0][2]) * cd18;
    t[4] = f1 * (p[9] * state.plg[0][1] + p[10] * state.plg[0][3]) * cd14;
    t[5] = p[37] * state.plg[0][1] * cd39;
    if flags.sw[7] != 0.0 {
        let t71 = p[11] * state.plg[1][2] * cd14 * flags.swc[5];
        let t72 = p[12] * state.plg[1][2] * cd14 * flags.swc[5];
        t[6] = f2
            * ((p[3] * state.plg[1][1] + p[4] * state.plg[1][3] + p[27] * state.plg[1][5] + t71)
                * state.ctloc
                + (p[6] * state.plg[1][1]
                    + p[7] * state.plg[1][3]
                    + p[28] * state.plg[1][5]
                    + t72)
                    * state.stloc);
    }
    if flags.sw[8] != 0.0 {
        let t81 = (p[23] * state.plg[2][3] + p[35] * state.plg[2][5]) * cd14 * flags.swc[5];
        let t82 = (p[33] * state.plg[2][3] + p[36] * state.plg[2][5]) * cd14 * flags.swc[5];
        t[7] = f2
            * ((p[5] * state.plg[2][2] + p[41] * state.plg[2][4] + t81) * state.c2tloc
                + (p[8] * state.plg[2][2] + p[42] * state.plg[2][4] + t82) * state.s2tloc);
    }
    if flags.sw[14] != 0.0 {
        t[13] = f2
            * ((p[39] * state.plg[3][3]
                + (p[93] * state.plg[3][4] + p[46] * state.plg[3][6]) * cd14 * flags.swc[5])
                * state.s3tloc
                + (p[40] * state.plg[3][3]
                    + (p[94] * state.plg[3][4] + p[48] * state.plg[3][6]) * cd14 * flags.swc[5])
                    * state.c3tloc);
    }
    if flags.sw[9] as i32 == -1 {
        if let Some(ref ap_a) = input.ap_a {
            if p[51] != 0.0 {
                let exp1_arg = -10800.0 * (p[51] * p[51]).sqrt()
                    / (1.0 + p[138] * (45.0 - (input.g_lat * input.g_lat).sqrt()));
                let mut exp1 = exp1_arg.exp();
                if exp1 > 0.99999 {
                    exp1 = 0.99999;
                }
                state.apt[0] = sg0(exp1, p, ap_a);
                if flags.sw[9] != 0.0 {
                    t[8] = state.apt[0]
                        * (p[50]
                            + p[96] * state.plg[0][2]
                            + p[54] * state.plg[0][4]
                            + (p[125] * state.plg[0][1]
                                + p[126] * state.plg[0][3]
                                + p[127] * state.plg[0][5])
                                * cd14
                                * flags.swc[5]
                            + (p[128] * state.plg[1][1]
                                + p[129] * state.plg[1][3]
                                + p[130] * state.plg[1][5])
                                * flags.swc[7]
                                * (hr * (tloc - p[131])).cos());
                }
            }
        }
    } else {
        let apd = input.ap - 4.0;
        let mut p44 = p[43];
        let p45 = p[44];
        if p44 < 0.0 {
            p44 = 1.0E-5;
        }
        state.apdf = apd + (p45 - 1.0) * (apd + ((-p44 * apd).exp() - 1.0) / p44);
        if flags.sw[9] != 0.0 {
            t[8] = state.apdf
                * (p[32]
                    + p[45] * state.plg[0][2]
                    + p[34] * state.plg[0][4]
                    + (p[100] * state.plg[0][1]
                        + p[101] * state.plg[0][3]
                        + p[102] * state.plg[0][5])
                        * cd14
                        * flags.swc[5]
                    + (p[121] * state.plg[1][1]
                        + p[122] * state.plg[1][3]
                        + p[123] * state.plg[1][5])
                        * flags.swc[7]
                        * (hr * (tloc - p[124])).cos());
        }
    }
    if flags.sw[10] != 0.0 && input.g_long > -1000.0 {
        if flags.sw[11] != 0.0 {
            t[10] = (1.0 + p[80] * state.dfa * flags.swc[1])
                * ((p[64] * state.plg[1][2]
                    + p[65] * state.plg[1][4]
                    + p[66] * state.plg[1][6]
                    + p[103] * state.plg[1][1]
                    + p[104] * state.plg[1][3]
                    + p[105] * state.plg[1][5]
                    + flags.swc[5]
                        * (p[109] * state.plg[1][1]
                            + p[110] * state.plg[1][3]
                            + p[111] * state.plg[1][5])
                        * cd14)
                    * (dgtr * input.g_long).cos()
                    + (p[90] * state.plg[1][2]
                        + p[91] * state.plg[1][4]
                        + p[92] * state.plg[1][6]
                        + p[106] * state.plg[1][1]
                        + p[107] * state.plg[1][3]
                        + p[108] * state.plg[1][5]
                        + flags.swc[5]
                            * (p[112] * state.plg[1][1]
                                + p[113] * state.plg[1][3]
                                + p[114] * state.plg[1][5])
                            * cd14)
                        * (dgtr * input.g_long).sin());
        }
        if flags.sw[12] != 0.0 {
            t[11] = (1.0 + p[95] * state.plg[0][1])
                * (1.0 + p[81] * state.dfa * flags.swc[1])
                * (1.0 + p[119] * state.plg[0][1] * flags.swc[5] * cd14)
                * ((p[68] * state.plg[0][1] + p[69] * state.plg[0][3] + p[70] * state.plg[0][5])
                    * (sr * (input.sec - p[71])).cos());
            t[11] += flags.swc[11]
                * (p[76] * state.plg[2][3] + p[77] * state.plg[2][5] + p[78] * state.plg[2][7])
                * (sr * (input.sec - p[79]) + 2.0 * dgtr * input.g_long).cos()
                * (1.0 + p[137] * state.dfa * flags.swc[1]);
        }
        if flags.sw[13] != 0.0 {
            if flags.sw[9] as i32 == -1 {
                if let Some(ref _ap_a) = input.ap_a {
                    if p[51] != 0.0 {
                        t[12] = state.apt[0]
                            * flags.swc[11]
                            * (1.0 + p[132] * state.plg[0][1])
                            * ((p[52] * state.plg[1][2]
                                + p[98] * state.plg[1][4]
                                + p[67] * state.plg[1][6])
                                * (dgtr * (input.g_long - p[97])).cos())
                            + state.apt[0]
                                * flags.swc[11]
                                * flags.swc[5]
                                * (p[133] * state.plg[1][1]
                                    + p[134] * state.plg[1][3]
                                    + p[135] * state.plg[1][5])
                                * cd14
                                * (dgtr * (input.g_long - p[136])).cos()
                            + state.apt[0]
                                * flags.swc[12]
                                * (p[55] * state.plg[0][1]
                                    + p[56] * state.plg[0][3]
                                    + p[57] * state.plg[0][5])
                                * (sr * (input.sec - p[58])).cos();
                    }
                }
            } else {
                t[12] = state.apdf
                    * flags.swc[11]
                    * (1.0 + p[120] * state.plg[0][1])
                    * ((p[60] * state.plg[1][2]
                        + p[61] * state.plg[1][4]
                        + p[62] * state.plg[1][6])
                        * (dgtr * (input.g_long - p[63])).cos())
                    + state.apdf
                        * flags.swc[11]
                        * flags.swc[5]
                        * (p[115] * state.plg[1][1]
                            + p[116] * state.plg[1][3]
                            + p[117] * state.plg[1][5])
                        * cd14
                        * (dgtr * (input.g_long - p[118])).cos()
                    + state.apdf
                        * flags.swc[12]
                        * (p[83] * state.plg[0][1]
                            + p[84] * state.plg[0][3]
                            + p[85] * state.plg[0][5])
                        * (sr * (input.sec - p[75])).cos();
            }
        }
    }
    let mut tinf = p[30];
    for i in 0..14 {
        tinf += flags.sw[i + 1].abs() * t[i];
    }
    tinf
}

fn glob7s(p: &[f64], input: &NrlmsiseInput, flags: &NrlmsiseFlags, state: &NrlmsiseState) -> f64 {
    let dr = 1.72142E-2_f64;
    let dgtr = 1.74533E-2_f64;
    let _hr = 0.2618_f64;
    let mut t = [0.0_f64; 14];
    let cd32 = (dr * (input.doy as f64 - p[31])).cos();
    let cd18 = (2.0 * dr * (input.doy as f64 - p[17])).cos();
    let cd14 = (dr * (input.doy as f64 - p[13])).cos();
    let cd39 = (2.0 * dr * (input.doy as f64 - p[38])).cos();
    t[0] = p[21] * state.dfa;
    t[1] = p[1] * state.plg[0][2]
        + p[2] * state.plg[0][4]
        + p[22] * state.plg[0][6]
        + p[26] * state.plg[0][1]
        + p[14] * state.plg[0][3]
        + p[59] * state.plg[0][5];
    t[2] = (p[18] + p[47] * state.plg[0][2] + p[29] * state.plg[0][4]) * cd32;
    t[3] = (p[15] + p[16] * state.plg[0][2] + p[30] * state.plg[0][4]) * cd18;
    t[4] = (p[9] * state.plg[0][1] + p[10] * state.plg[0][3] + p[20] * state.plg[0][5]) * cd14;
    t[5] = p[37] * state.plg[0][1] * cd39;
    if flags.sw[7] != 0.0 {
        let t71 = p[11] * state.plg[1][2] * cd14 * flags.swc[5];
        let t72 = p[12] * state.plg[1][2] * cd14 * flags.swc[5];
        t[6] = (p[3] * state.plg[1][1] + p[4] * state.plg[1][3] + t71) * state.ctloc
            + (p[6] * state.plg[1][1] + p[7] * state.plg[1][3] + t72) * state.stloc;
    }
    if flags.sw[8] != 0.0 {
        let t81 = (p[23] * state.plg[2][3] + p[35] * state.plg[2][5]) * cd14 * flags.swc[5];
        let t82 = (p[33] * state.plg[2][3] + p[36] * state.plg[2][5]) * cd14 * flags.swc[5];
        t[7] = (p[5] * state.plg[2][2] + p[41] * state.plg[2][4] + t81) * state.c2tloc
            + (p[8] * state.plg[2][2] + p[42] * state.plg[2][4] + t82) * state.s2tloc;
    }
    if flags.sw[14] != 0.0 {
        t[13] = p[39] * state.plg[3][3] * state.s3tloc + p[40] * state.plg[3][3] * state.c3tloc;
    }
    if flags.sw[9] != 0.0 {
        if flags.sw[9] as i32 == 1 {
            t[8] = state.apdf * (p[32] + p[45] * state.plg[0][2] * flags.swc[2]);
        }
        if flags.sw[9] as i32 == -1 {
            t[8] = p[50] * state.apt[0] + p[96] * state.plg[0][2] * state.apt[0] * flags.swc[2];
        }
    }
    if flags.sw[10] != 0.0 && flags.sw[11] != 0.0 && input.g_long > -1000.0 {
        t[10] = (1.0
            + state.plg[0][1]
                * (p[80] * flags.swc[5] * (dr * (input.doy as f64 - p[81])).cos()
                    + p[85] * flags.swc[6] * (2.0 * dr * (input.doy as f64 - p[86])).cos())
            + p[83] * flags.swc[3] * (dr * (input.doy as f64 - p[84])).cos()
            + p[87] * flags.swc[4] * (2.0 * dr * (input.doy as f64 - p[88])).cos())
            * ((p[64] * state.plg[1][2]
                + p[65] * state.plg[1][4]
                + p[66] * state.plg[1][6]
                + p[74] * state.plg[1][1]
                + p[75] * state.plg[1][3]
                + p[76] * state.plg[1][5])
                * (dgtr * input.g_long).cos()
                + (p[90] * state.plg[1][2]
                    + p[91] * state.plg[1][4]
                    + p[92] * state.plg[1][6]
                    + p[77] * state.plg[1][1]
                    + p[78] * state.plg[1][3]
                    + p[79] * state.plg[1][5])
                    * (dgtr * input.g_long).sin());
    }
    let mut tt = 0.0_f64;
    for i in 0..14 {
        tt += flags.sw[i + 1].abs() * t[i];
    }
    tt
}

fn gts7(
    input: &mut NrlmsiseInput,
    flags: &NrlmsiseFlags,
    output: &mut NrlmsiseOutput,
    state: &mut NrlmsiseState,
) {
    let dgtr = 1.74533E-2_f64;
    let dr = 1.72142E-2_f64;
    let alpha = [-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0_f64];
    let altl = [200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0_f64];
    let za = PDL[1][15];
    let zn1 = [za, 110.0, 100.0, 90.0, 72.5];
    let mn1 = 5_usize;
    for j in 0..9 {
        output.d[j] = 0.0;
    }
    let tinf = if input.alt > zn1[0] {
        PTM[0] * PT[0] * (1.0 + flags.sw[16] * globe7(&PT, input, flags, state))
    } else {
        PTM[0] * PT[0]
    };
    output.t[0] = tinf;
    let g0 = if input.alt > zn1[4] {
        PTM[3] * PS[0] * (1.0 + flags.sw[19] * globe7(&PS, input, flags, state))
    } else {
        PTM[3] * PS[0]
    };
    let tlb = PTM[1] * (1.0 + flags.sw[17] * globe7(&PD[3], input, flags, state)) * PD[3][0];
    let s = g0 / (tinf - tlb);
    if input.alt < 300.0 {
        state.meso_tn1[1] =
            PTM[6] * PTL[0][0] / (1.0 - flags.sw[18] * glob7s(&PTL[0], input, flags, state));
        state.meso_tn1[2] =
            PTM[2] * PTL[1][0] / (1.0 - flags.sw[18] * glob7s(&PTL[1], input, flags, state));
        state.meso_tn1[3] =
            PTM[7] * PTL[2][0] / (1.0 - flags.sw[18] * glob7s(&PTL[2], input, flags, state));
        state.meso_tn1[4] = PTM[4] * PTL[3][0]
            / (1.0 - flags.sw[18] * flags.sw[20] * glob7s(&PTL[3], input, flags, state));
        state.meso_tgn1[1] = PTM[8]
            * PMA[8][0]
            * (1.0 + flags.sw[18] * flags.sw[20] * glob7s(&PMA[8], input, flags, state))
            * state.meso_tn1[4]
            * state.meso_tn1[4]
            / ((PTM[4] * PTL[3][0]) * (PTM[4] * PTL[3][0]));
    } else {
        state.meso_tn1[1] = PTM[6] * PTL[0][0];
        state.meso_tn1[2] = PTM[2] * PTL[1][0];
        state.meso_tn1[3] = PTM[7] * PTL[2][0];
        state.meso_tn1[4] = PTM[4] * PTL[3][0];
        state.meso_tgn1[1] = PTM[8] * PMA[8][0] * state.meso_tn1[4] * state.meso_tn1[4]
            / ((PTM[4] * PTL[3][0]) * (PTM[4] * PTL[3][0]));
    }
    let g28 = flags.sw[21] * globe7(&PD[2], input, flags, state);
    let zhf = PDL[1][24]
        * (1.0
            + flags.sw[5]
                * PDL[0][24]
                * (dgtr * input.g_lat).sin()
                * (dr * (input.doy as f64 - PT[13])).cos());
    output.t[0] = tinf;
    let xmm = PDM[2][4];
    let z = input.alt;
    let db28 = PDM[2][0] * g28.exp() * PD[2][0];
    output.d[2] = densu(
        z,
        db28,
        tinf,
        tlb,
        28.0,
        alpha[2],
        &mut output.t[1],
        PTM[5],
        s,
        mn1,
        &zn1,
        &mut state.meso_tn1,
        &mut state.meso_tgn1,
        state.gsurf,
        state.re,
    );
    state.dd = output.d[2];
    let zh28 = PDM[2][2] * zhf;
    let zhm28 = PDM[2][3] * PDL[1][5];
    let xmd = 28.0 - xmm;
    let b28 = densu(
        zh28,
        db28,
        tinf,
        tlb,
        xmd,
        alpha[2] - 1.0,
        &mut output.t[1],
        PTM[5],
        s,
        mn1,
        &zn1,
        &mut state.meso_tn1,
        &mut state.meso_tgn1,
        state.gsurf,
        state.re,
    );
    if flags.sw[15] != 0.0 && z <= altl[2] {
        state.dm28 = densu(
            z,
            b28,
            tinf,
            tlb,
            xmm,
            alpha[2],
            &mut output.t[1],
            PTM[5],
            s,
            mn1,
            &zn1,
            &mut state.meso_tn1,
            &mut state.meso_tgn1,
            state.gsurf,
            state.re,
        );
        output.d[2] = dnet(output.d[2], state.dm28, zhm28, xmm, 28.0);
    }
    let g4 = flags.sw[21] * globe7(&PD[0], input, flags, state);
    let db04 = PDM[0][0] * g4.exp() * PD[0][0];
    output.d[0] = densu(
        z,
        db04,
        tinf,
        tlb,
        4.0,
        alpha[0],
        &mut output.t[1],
        PTM[5],
        s,
        mn1,
        &zn1,
        &mut state.meso_tn1,
        &mut state.meso_tgn1,
        state.gsurf,
        state.re,
    );
    state.dd = output.d[0];
    if flags.sw[15] != 0.0 && z < altl[0] {
        let zh04 = PDM[0][2];
        let b04 = densu(
            zh04,
            db04,
            tinf,
            tlb,
            4.0 - xmm,
            alpha[0] - 1.0,
            &mut output.t[1],
            PTM[5],
            s,
            mn1,
            &zn1,
            &mut state.meso_tn1,
            &mut state.meso_tgn1,
            state.gsurf,
            state.re,
        );
        state.dm04 = densu(
            z,
            b04,
            tinf,
            tlb,
            xmm,
            0.0,
            &mut output.t[1],
            PTM[5],
            s,
            mn1,
            &zn1,
            &mut state.meso_tn1,
            &mut state.meso_tgn1,
            state.gsurf,
            state.re,
        );
        output.d[0] = dnet(output.d[0], state.dm04, zhm28, xmm, 4.0);
        let rl = (b28 * PDM[0][1] / b04).ln();
        let zc04 = PDM[0][4] * PDL[1][0];
        let hc04 = PDM[0][5] * PDL[1][1];
        output.d[0] *= ccor(z, rl, hc04, zc04);
    }
    let g16 = flags.sw[21] * globe7(&PD[1], input, flags, state);
    let db16 = PDM[1][0] * g16.exp() * PD[1][0];
    output.d[1] = densu(
        z,
        db16,
        tinf,
        tlb,
        16.0,
        alpha[1],
        &mut output.t[1],
        PTM[5],
        s,
        mn1,
        &zn1,
        &mut state.meso_tn1,
        &mut state.meso_tgn1,
        state.gsurf,
        state.re,
    );
    state.dd = output.d[1];
    if flags.sw[15] != 0.0 && z <= altl[1] {
        let zh16 = PDM[1][2];
        let b16 = densu(
            zh16,
            db16,
            tinf,
            tlb,
            16.0 - xmm,
            alpha[1] - 1.0,
            &mut output.t[1],
            PTM[5],
            s,
            mn1,
            &zn1,
            &mut state.meso_tn1,
            &mut state.meso_tgn1,
            state.gsurf,
            state.re,
        );
        state.dm16 = densu(
            z,
            b16,
            tinf,
            tlb,
            xmm,
            0.0,
            &mut output.t[1],
            PTM[5],
            s,
            mn1,
            &zn1,
            &mut state.meso_tn1,
            &mut state.meso_tgn1,
            state.gsurf,
            state.re,
        );
        output.d[1] = dnet(output.d[1], state.dm16, zhm28, xmm, 16.0);
        let rl = PDM[1][1] * PDL[1][16] * (1.0 + flags.sw[1] * PDL[0][23] * (input.f107a - 150.0));
        let hc16 = PDM[1][5] * PDL[1][3];
        let zc16 = PDM[1][4] * PDL[1][2];
        let hc216 = PDM[1][5] * PDL[1][4];
        output.d[1] *= ccor2(z, rl, hc16, zc16, hc216);
        let hcc16 = PDM[1][7] * PDL[1][13];
        let zcc16 = PDM[1][6] * PDL[1][12];
        let rc16 = PDM[1][3] * PDL[1][14];
        output.d[1] *= ccor(z, rc16, hcc16, zcc16);
    }
    let g32 = flags.sw[21] * globe7(&PD[4], input, flags, state);
    let db32 = PDM[3][0] * g32.exp() * PD[4][0];
    output.d[3] = densu(
        z,
        db32,
        tinf,
        tlb,
        32.0,
        alpha[3],
        &mut output.t[1],
        PTM[5],
        s,
        mn1,
        &zn1,
        &mut state.meso_tn1,
        &mut state.meso_tgn1,
        state.gsurf,
        state.re,
    );
    state.dd = output.d[3];
    if flags.sw[15] != 0.0 {
        if z <= altl[3] {
            let zh32 = PDM[3][2];
            let b32 = densu(
                zh32,
                db32,
                tinf,
                tlb,
                32.0 - xmm,
                alpha[3] - 1.0,
                &mut output.t[1],
                PTM[5],
                s,
                mn1,
                &zn1,
                &mut state.meso_tn1,
                &mut state.meso_tgn1,
                state.gsurf,
                state.re,
            );
            state.dm32 = densu(
                z,
                b32,
                tinf,
                tlb,
                xmm,
                0.0,
                &mut output.t[1],
                PTM[5],
                s,
                mn1,
                &zn1,
                &mut state.meso_tn1,
                &mut state.meso_tgn1,
                state.gsurf,
                state.re,
            );
            output.d[3] = dnet(output.d[3], state.dm32, zhm28, xmm, 32.0);
            let rl = (b28 * PDM[3][1] / b32).ln();
            let hc32 = PDM[3][5] * PDL[1][7];
            let zc32 = PDM[3][4] * PDL[1][6];
            output.d[3] *= ccor(z, rl, hc32, zc32);
        }
        let hcc32 = PDM[3][7] * PDL[1][22];
        let hcc232 = PDM[3][7] * PDL[0][22];
        let zcc32 = PDM[3][6] * PDL[1][21];
        let rc32 =
            PDM[3][3] * PDL[1][23] * (1.0 + flags.sw[1] * PDL[0][23] * (input.f107a - 150.0));
        output.d[3] *= ccor2(z, rc32, hcc32, zcc32, hcc232);
    }
    let g40 = flags.sw[21] * globe7(&PD[5], input, flags, state);
    let db40 = PDM[4][0] * g40.exp() * PD[5][0];
    output.d[4] = densu(
        z,
        db40,
        tinf,
        tlb,
        40.0,
        alpha[4],
        &mut output.t[1],
        PTM[5],
        s,
        mn1,
        &zn1,
        &mut state.meso_tn1,
        &mut state.meso_tgn1,
        state.gsurf,
        state.re,
    );
    state.dd = output.d[4];
    if flags.sw[15] != 0.0 && z <= altl[4] {
        let zh40 = PDM[4][2];
        let b40 = densu(
            zh40,
            db40,
            tinf,
            tlb,
            40.0 - xmm,
            alpha[4] - 1.0,
            &mut output.t[1],
            PTM[5],
            s,
            mn1,
            &zn1,
            &mut state.meso_tn1,
            &mut state.meso_tgn1,
            state.gsurf,
            state.re,
        );
        state.dm40 = densu(
            z,
            b40,
            tinf,
            tlb,
            xmm,
            0.0,
            &mut output.t[1],
            PTM[5],
            s,
            mn1,
            &zn1,
            &mut state.meso_tn1,
            &mut state.meso_tgn1,
            state.gsurf,
            state.re,
        );
        output.d[4] = dnet(output.d[4], state.dm40, zhm28, xmm, 40.0);
        let rl = (b28 * PDM[4][1] / b40).ln();
        let hc40 = PDM[4][5] * PDL[1][9];
        let zc40 = PDM[4][4] * PDL[1][8];
        output.d[4] *= ccor(z, rl, hc40, zc40);
    }
    let g1 = flags.sw[21] * globe7(&PD[6], input, flags, state);
    let db01 = PDM[5][0] * g1.exp() * PD[6][0];
    output.d[6] = densu(
        z,
        db01,
        tinf,
        tlb,
        1.0,
        alpha[6],
        &mut output.t[1],
        PTM[5],
        s,
        mn1,
        &zn1,
        &mut state.meso_tn1,
        &mut state.meso_tgn1,
        state.gsurf,
        state.re,
    );
    state.dd = output.d[6];
    if flags.sw[15] != 0.0 && z <= altl[6] {
        let zh01 = PDM[5][2];
        let b01 = densu(
            zh01,
            db01,
            tinf,
            tlb,
            1.0 - xmm,
            alpha[6] - 1.0,
            &mut output.t[1],
            PTM[5],
            s,
            mn1,
            &zn1,
            &mut state.meso_tn1,
            &mut state.meso_tgn1,
            state.gsurf,
            state.re,
        );
        state.dm01 = densu(
            z,
            b01,
            tinf,
            tlb,
            xmm,
            0.0,
            &mut output.t[1],
            PTM[5],
            s,
            mn1,
            &zn1,
            &mut state.meso_tn1,
            &mut state.meso_tgn1,
            state.gsurf,
            state.re,
        );
        output.d[6] = dnet(output.d[6], state.dm01, zhm28, xmm, 1.0);
        let rl = (b28 * PDM[5][1] * (PDL[1][17] * PDL[1][17]).sqrt() / b01).ln();
        let hc01 = PDM[5][5] * PDL[1][11];
        let zc01 = PDM[5][4] * PDL[1][10];
        output.d[6] *= ccor(z, rl, hc01, zc01);
        let hcc01 = PDM[5][7] * PDL[1][19];
        let zcc01 = PDM[5][6] * PDL[1][18];
        let rc01 = PDM[5][3] * PDL[1][20];
        output.d[6] *= ccor(z, rc01, hcc01, zcc01);
    }
    let g14 = flags.sw[21] * globe7(&PD[7], input, flags, state);
    let db14 = PDM[6][0] * g14.exp() * PD[7][0];
    output.d[7] = densu(
        z,
        db14,
        tinf,
        tlb,
        14.0,
        alpha[7],
        &mut output.t[1],
        PTM[5],
        s,
        mn1,
        &zn1,
        &mut state.meso_tn1,
        &mut state.meso_tgn1,
        state.gsurf,
        state.re,
    );
    state.dd = output.d[7];
    if flags.sw[15] != 0.0 && z <= altl[7] {
        let zh14 = PDM[6][2];
        let b14 = densu(
            zh14,
            db14,
            tinf,
            tlb,
            14.0 - xmm,
            alpha[7] - 1.0,
            &mut output.t[1],
            PTM[5],
            s,
            mn1,
            &zn1,
            &mut state.meso_tn1,
            &mut state.meso_tgn1,
            state.gsurf,
            state.re,
        );
        state.dm14 = densu(
            z,
            b14,
            tinf,
            tlb,
            xmm,
            0.0,
            &mut output.t[1],
            PTM[5],
            s,
            mn1,
            &zn1,
            &mut state.meso_tn1,
            &mut state.meso_tgn1,
            state.gsurf,
            state.re,
        );
        output.d[7] = dnet(output.d[7], state.dm14, zhm28, xmm, 14.0);
        let rl = (b28 * PDM[6][1] * (PDL[0][2] * PDL[0][2]).sqrt() / b14).ln();
        let hc14 = PDM[6][5] * PDL[0][1];
        let zc14 = PDM[6][4] * PDL[0][0];
        output.d[7] *= ccor(z, rl, hc14, zc14);
        let hcc14 = PDM[6][7] * PDL[0][4];
        let zcc14 = PDM[6][6] * PDL[0][3];
        let rc14 = PDM[6][3] * PDL[0][5];
        output.d[7] *= ccor(z, rc14, hcc14, zcc14);
    }
    let g16h = flags.sw[21] * globe7(&PD[8], input, flags, state);
    let db16h = PDM[7][0] * g16h.exp() * PD[8][0];
    let tho = PDM[7][9] * PDL[0][6];
    let dd_anom = densu(
        z,
        db16h,
        tho,
        tho,
        16.0,
        alpha[8],
        &mut output.t[1],
        PTM[5],
        s,
        mn1,
        &zn1,
        &mut state.meso_tn1,
        &mut state.meso_tgn1,
        state.gsurf,
        state.re,
    );
    let zsht = PDM[7][5];
    let zmho = PDM[7][4];
    let zsho = scalh(zmho, 16.0, tho, state.gsurf, state.re);
    output.d[8] = dd_anom * (-zsht / zsho * ((-(z - zmho) / zsht).exp() - 1.0)).exp();
    output.d[5] = 1.66E-24
        * (4.0 * output.d[0]
            + 16.0 * output.d[1]
            + 28.0 * output.d[2]
            + 32.0 * output.d[3]
            + 40.0 * output.d[4]
            + output.d[6]
            + 14.0 * output.d[7]);
    let z_abs = (input.alt * input.alt).sqrt();
    let _ddum = densu(
        z_abs,
        1.0,
        tinf,
        tlb,
        0.0,
        0.0,
        &mut output.t[1],
        PTM[5],
        s,
        mn1,
        &zn1,
        &mut state.meso_tn1,
        &mut state.meso_tgn1,
        state.gsurf,
        state.re,
    );
    if flags.sw[0] != 0.0 {
        for i in 0..9 {
            output.d[i] *= 1.0E6;
        }
        output.d[5] /= 1000.0;
    }
}

fn gtd7(
    input: &mut NrlmsiseInput,
    flags: &mut NrlmsiseFlags,
    output: &mut NrlmsiseOutput,
    state: &mut NrlmsiseState,
) {
    let mn3 = 5_usize;
    let zn3 = [32.5, 20.0, 15.0, 10.0, 0.0_f64];
    let mn2 = 4_usize;
    let zn2 = [72.5, 55.0, 45.0, 32.5_f64];
    let zmix = 62.5_f64;
    tselec(flags);
    let xlat = if flags.sw[2] == 0.0 {
        45.0
    } else {
        input.g_lat
    };
    glatf(xlat, &mut state.gsurf, &mut state.re);
    let xmm = PDM[2][4];
    let altt = if input.alt > zn2[0] {
        input.alt
    } else {
        zn2[0]
    };
    let tmp = input.alt;
    input.alt = altt;
    let mut soutput = NrlmsiseOutput {
        d: [0.0; 9],
        t: [0.0; 2],
    };
    gts7(input, flags, &mut soutput, state);
    input.alt = tmp;
    let dm28m = if flags.sw[0] != 0.0 {
        state.dm28 * 1.0E6
    } else {
        state.dm28
    };
    output.t[0] = soutput.t[0];
    output.t[1] = soutput.t[1];
    if input.alt >= zn2[0] {
        for i in 0..9 {
            output.d[i] = soutput.d[i];
        }
        return;
    }
    state.meso_tgn2[0] = state.meso_tgn1[1];
    state.meso_tn2[0] = state.meso_tn1[4];
    state.meso_tn2[1] =
        PMA[0][0] * PAVGM[0] / (1.0 - flags.sw[20] * glob7s(&PMA[0], input, flags, state));
    state.meso_tn2[2] =
        PMA[1][0] * PAVGM[1] / (1.0 - flags.sw[20] * glob7s(&PMA[1], input, flags, state));
    state.meso_tn2[3] = PMA[2][0] * PAVGM[2]
        / (1.0 - flags.sw[20] * flags.sw[22] * glob7s(&PMA[2], input, flags, state));
    state.meso_tgn2[1] = PAVGM[8]
        * PMA[9][0]
        * (1.0 + flags.sw[20] * flags.sw[22] * glob7s(&PMA[9], input, flags, state))
        * state.meso_tn2[3]
        * state.meso_tn2[3]
        / ((PMA[2][0] * PAVGM[2]) * (PMA[2][0] * PAVGM[2]));
    state.meso_tn3[0] = state.meso_tn2[3];
    if input.alt <= zn3[0] {
        state.meso_tgn3[0] = state.meso_tgn2[1];
        state.meso_tn3[1] =
            PMA[3][0] * PAVGM[3] / (1.0 - flags.sw[22] * glob7s(&PMA[3], input, flags, state));
        state.meso_tn3[2] =
            PMA[4][0] * PAVGM[4] / (1.0 - flags.sw[22] * glob7s(&PMA[4], input, flags, state));
        state.meso_tn3[3] =
            PMA[5][0] * PAVGM[5] / (1.0 - flags.sw[22] * glob7s(&PMA[5], input, flags, state));
        state.meso_tn3[4] =
            PMA[6][0] * PAVGM[6] / (1.0 - flags.sw[22] * glob7s(&PMA[6], input, flags, state));
        state.meso_tgn3[1] = PMA[7][0]
            * PAVGM[7]
            * (1.0 + flags.sw[22] * glob7s(&PMA[7], input, flags, state))
            * state.meso_tn3[4]
            * state.meso_tn3[4]
            / ((PMA[6][0] * PAVGM[6]) * (PMA[6][0] * PAVGM[6]));
    }
    let dmc = if input.alt > zmix {
        1.0 - (zn2[0] - input.alt) / (zn2[0] - zmix)
    } else {
        0.0
    };
    let dz28 = soutput.d[2];
    let mut tz = 0.0_f64;
    let dmr = soutput.d[2] / dm28m - 1.0;
    output.d[2] = densm(
        input.alt,
        dm28m,
        xmm,
        &mut tz,
        mn3,
        &zn3,
        &state.meso_tn3,
        &state.meso_tgn3,
        mn2,
        &zn2,
        &state.meso_tn2,
        &state.meso_tgn2,
        state.gsurf,
        state.re,
    );
    output.d[2] *= 1.0 + dmr * dmc;
    let dmr = soutput.d[0] / (dz28 * PDM[0][1]) - 1.0;
    output.d[0] = output.d[2] * PDM[0][1] * (1.0 + dmr * dmc);
    output.d[1] = 0.0;
    output.d[8] = 0.0;
    let dmr = soutput.d[3] / (dz28 * PDM[3][1]) - 1.0;
    output.d[3] = output.d[2] * PDM[3][1] * (1.0 + dmr * dmc);
    let dmr = soutput.d[4] / (dz28 * PDM[4][1]) - 1.0;
    output.d[4] = output.d[2] * PDM[4][1] * (1.0 + dmr * dmc);
    output.d[6] = 0.0;
    output.d[7] = 0.0;
    output.d[5] = 1.66E-24
        * (4.0 * output.d[0]
            + 16.0 * output.d[1]
            + 28.0 * output.d[2]
            + 32.0 * output.d[3]
            + 40.0 * output.d[4]
            + output.d[6]
            + 14.0 * output.d[7]);
    if flags.sw[0] != 0.0 {
        output.d[5] /= 1000.0;
    }
    let _dd = densm(
        input.alt,
        1.0,
        0.0,
        &mut tz,
        mn3,
        &zn3,
        &state.meso_tn3,
        &state.meso_tgn3,
        mn2,
        &zn2,
        &state.meso_tn2,
        &state.meso_tgn2,
        state.gsurf,
        state.re,
    );
    output.t[1] = tz;
}

fn gtd7d(
    input: &mut NrlmsiseInput,
    flags: &mut NrlmsiseFlags,
    output: &mut NrlmsiseOutput,
    state: &mut NrlmsiseState,
) {
    gtd7(input, flags, output, state);
    output.d[5] = 1.66E-24
        * (4.0 * output.d[0]
            + 16.0 * output.d[1]
            + 28.0 * output.d[2]
            + 32.0 * output.d[3]
            + 40.0 * output.d[4]
            + output.d[6]
            + 14.0 * output.d[7]
            + 16.0 * output.d[8]);
    if flags.sw[0] != 0.0 {
        output.d[5] /= 1000.0;
    }
}

/// NRL MSISE-00 model for atmosphere density
///
/// # Arguments
///
///   * `alt_km` -  Altitude in kilometers
///   * `lat_option` - Optional latitude in degrees (default: 0)
///   * `lon_option` - Optional longitude in degrees (default: 0)
///   * `time_option` -  Optional time, for when using space weather
///   * `use_spaceweather` -  Boolean to use space weather data
///
/// # Outputs
///
///   Tuple with two elements:
///      * Atmosphere mass density in kg / m^3
///      * Temperature in Kelvin
///
pub fn nrlmsise(
    alt_km: f64,
    lat_option: Option<f64>,
    lon_option: Option<f64>,
    time_option: Option<&impl TimeLike>,
    use_spaceweather: bool,
) -> (f64, f64) {
    let lat = lat_option.unwrap_or(0.0);
    let lon = lon_option.unwrap_or(0.0);
    let mut day_of_year: i32 = 1;
    let mut sec_of_day: f64 = 0.0;
    let mut f107a: f64 = 150.0;
    let mut f107: f64 = 150.0;
    let mut ap: f64 = 4.0;
    if let Some(time) = time_option {
        let time = time.as_instant();
        let (year, _mon, _day, dhour, dmin, dsec) = time.as_datetime();
        let fday: f64 = (time - Instant::from_date(year, 1, 1).unwrap()).as_days() + 1.0;
        day_of_year = fday.floor() as i32;
        sec_of_day = (dhour as f64).mul_add(3600.0, dmin as f64 * 60.0) + dsec;
        if use_spaceweather {
            let sw_time = time - Duration::from_days(1.0);
            if let Ok(r) = spaceweather::get(&sw_time) {
                f107a = r.f10p7_adj_c81;
                f107 = r.f10p7_adj;
                ap = r.ap_avg as f64;
            } else if let Some(predicted) = solar_cycle_forecast::get_predicted_f107(&time) {
                f107 = predicted;
                f107a = predicted;
            }
        }
    }
    let mut state = NrlmsiseState::new();
    let mut input = NrlmsiseInput {
        year: 2022,
        doy: day_of_year,
        sec: sec_of_day,
        alt: alt_km,
        g_lat: lat,
        g_long: lon,
        lst: sec_of_day / 3600.0 + lon / 15.0,
        f107a,
        f107,
        ap,
        ap_a: None,
    };
    let mut flags = NrlmsiseFlags {
        switches: [
            0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        ],
        sw: [0.0; 24],
        swc: [0.0; 24],
    };
    let mut output = NrlmsiseOutput {
        d: [0.0; 9],
        t: [0.0; 2],
    };
    gtd7d(&mut input, &mut flags, &mut output, &mut state);
    (output.d[5] * 1.0e3, output.t[1])
}

#[cfg(test)]
mod tests {
    use super::*;

    fn call_model(
        alt: f64,
        lat: f64,
        lon: f64,
        doy: i32,
        sec: f64,
        f107: f64,
        f107a: f64,
        ap: f64,
    ) -> NrlmsiseOutput {
        let mut state = NrlmsiseState::new();
        let mut input = NrlmsiseInput {
            year: 2000,
            doy,
            sec,
            alt,
            g_lat: lat,
            g_long: lon,
            lst: sec / 3600.0 + lon / 15.0,
            f107a,
            f107,
            ap,
            ap_a: None,
        };
        let mut flags = NrlmsiseFlags {
            switches: [
                0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
            ],
            sw: [0.0; 24],
            swc: [0.0; 24],
        };
        let mut output = NrlmsiseOutput {
            d: [0.0; 9],
            t: [0.0; 2],
        };
        gtd7d(&mut input, &mut flags, &mut output, &mut state);
        output
    }

    fn sanity_check(output: &NrlmsiseOutput, alt: f64) {
        assert!(
            output.d[5] > 0.0,
            "Total density must be positive at alt={alt}"
        );
        assert!(
            output.t[1] > 100.0,
            "Temperature must be > 100K at alt={alt}"
        );
        assert!(
            output.t[1] < 3000.0,
            "Temperature must be < 3000K at alt={alt}"
        );
    }

    #[test]
    fn test_nrlmsise_basic() {
        let tm: Instant = Instant::from_date(2010, 1, 1).unwrap()
            + Duration::from_days(171.0 + 29000.0 / 86400.0);
        let (density, temperature) = nrlmsise(400.0, Some(60.0), Some(-70.0), Some(&tm), true);
        assert!(density > 0.0, "Density must be positive");
        assert!(temperature > 0.0, "Temperature must be positive");
    }

    #[test]
    fn test_multiple_altitudes() {
        let altitudes = [100.0, 200.0, 400.0, 600.0, 800.0];
        let mut prev_density = f64::MAX;
        for &alt in &altitudes {
            let out = call_model(alt, 60.0, -70.0, 172, 29000.0, 150.0, 150.0, 4.0);
            sanity_check(&out, alt);
            if alt > 100.0 {
                assert!(
                    out.d[5] < prev_density,
                    "Density should decrease with altitude: alt={alt}"
                );
            }
            prev_density = out.d[5];
        }
    }

    #[test]
    fn test_different_latitudes() {
        for &lat in &[0.0, 45.0, 90.0] {
            let out = call_model(400.0, lat, 0.0, 172, 29000.0, 150.0, 150.0, 4.0);
            sanity_check(&out, 400.0);
        }
    }

    #[test]
    fn test_solar_activity_levels() {
        let mut prev_density = 0.0_f64;
        for &f107 in &[70.0, 150.0, 250.0] {
            let out = call_model(400.0, 45.0, 0.0, 172, 29000.0, f107, f107, 4.0);
            sanity_check(&out, 400.0);
            if f107 > 70.0 {
                assert!(
                    out.d[5] > prev_density,
                    "Higher F10.7 should increase density at 400km: f107={f107}"
                );
            }
            prev_density = out.d[5];
        }
    }

    #[test]
    fn test_no_spaceweather() {
        let (density, temperature) =
            nrlmsise(400.0, Some(45.0), Some(0.0), None::<&Instant>, false);
        assert!(density > 0.0);
        assert!(temperature > 0.0);
    }

    #[test]
    fn test_low_altitude() {
        let out = call_model(10.0, 45.0, 0.0, 172, 29000.0, 150.0, 150.0, 4.0);
        sanity_check(&out, 10.0);
        let out_high = call_model(400.0, 45.0, 0.0, 172, 29000.0, 150.0, 150.0, 4.0);
        assert!(
            out.d[5] > out_high.d[5] * 1e6,
            "10km density should be >> 400km density"
        );
    }

    #[test]
    fn test_exospheric_temperature() {
        let out_low = call_model(500.0, 0.0, 0.0, 172, 43200.0, 70.0, 70.0, 4.0);
        let out_high = call_model(500.0, 0.0, 0.0, 172, 43200.0, 250.0, 250.0, 4.0);
        assert!(
            out_high.t[0] > out_low.t[0],
            "Exospheric temp should increase with F10.7"
        );
    }

    #[test]
    fn test_deterministic() {
        let out1 = call_model(400.0, 45.0, -70.0, 172, 29000.0, 150.0, 150.0, 4.0);
        let out2 = call_model(400.0, 45.0, -70.0, 172, 29000.0, 150.0, 150.0, 4.0);
        for i in 0..9 {
            assert_eq!(out1.d[i], out2.d[i], "Density d[{i}] not deterministic");
        }
        assert_eq!(out1.t[0], out2.t[0]);
        assert_eq!(out1.t[1], out2.t[1]);
    }

    #[test]
    fn test_reference_values_400km() {
        let out = call_model(400.0, 60.0, -70.0, 172, 29000.0, 150.0, 150.0, 4.0);
        assert!(
            out.d[5] > 1.0e-16 && out.d[5] < 1.0e-10,
            "Total density at 400km out of expected range: {}",
            out.d[5]
        );
        assert!(
            out.t[1] > 700.0 && out.t[1] < 1500.0,
            "Temperature at 400km out of expected range: {}",
            out.t[1]
        );
    }

    #[test]
    fn test_species_densities_positive() {
        let out = call_model(400.0, 45.0, 0.0, 172, 29000.0, 150.0, 150.0, 4.0);
        assert!(out.d[0] > 0.0, "He density should be positive");
        assert!(out.d[1] > 0.0, "O density should be positive");
        assert!(out.d[2] > 0.0, "N2 density should be positive");
        assert!(out.d[3] > 0.0, "O2 density should be positive");
    }

    #[test]
    fn test_altitude_sweep_monotonic() {
        let alts = [100.0, 150.0, 200.0, 300.0, 400.0, 500.0, 600.0, 800.0];
        let mut prev = f64::MAX;
        for &alt in &alts {
            let out = call_model(alt, 0.0, 0.0, 172, 43200.0, 150.0, 150.0, 4.0);
            sanity_check(&out, alt);
            if alt >= 150.0 {
                assert!(out.d[5] < prev, "Density not decreasing at alt={alt}");
            }
            prev = out.d[5];
        }
    }

    #[test]
    fn test_against_c_reference() {
        // Reference values generated from the original C NRLMSISE-00 implementation
        // (Dominik Brodowski, release 20041227). 22 test cases covering altitudes
        // 50-1000km, latitudes -45 to 90, F10.7 70-250, Ap 4-40, all times of day.
        // Format: (doy, sec, alt_km, lat, lon, f107a, f107, ap, density_kg_m3, temp_K)
        let ref_values: &[(i32, f64, f64, f64, f64, f64, f64, f64, f64, f64)] = &[
            (
                172,
                29000.0,
                100.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                3.566_331_077_034_03e-7,
                213.4853662589559,
            ),
            (
                172,
                29000.0,
                200.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                2.405112682924217e-10,
                968.827_485_252_76,
            ),
            (
                172,
                29000.0,
                300.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                1.753041634659097e-11,
                1083.565143501857,
            ),
            (
                172,
                29000.0,
                400.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                2.400187349450909e-12,
                1_098.248_715_501_85,
            ),
            (
                172,
                29000.0,
                500.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                4.547669148745933e-13,
                1100.244572186091,
            ),
            (
                172,
                29000.0,
                600.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                1.032_110_876_054_66e-13,
                1100.531951866164,
            ),
            (
                172,
                29000.0,
                800.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                8.209468510644924e-15,
                1100.582691109543,
            ),
            (
                172,
                29000.0,
                1000.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                1.574632267239295e-15,
                1100.584084143845,
            ),
            (
                172,
                29000.0,
                400.0,
                0.0,
                0.0,
                150.0,
                150.0,
                4.0,
                2.751878914860656e-12,
                929.3974079388325,
            ),
            (
                172,
                29000.0,
                400.0,
                45.0,
                0.0,
                150.0,
                150.0,
                4.0,
                3.395091151031206e-12,
                1085.384248466074,
            ),
            (
                172,
                29000.0,
                400.0,
                90.0,
                0.0,
                150.0,
                150.0,
                4.0,
                3.648579006145896e-12,
                1186.441874023936,
            ),
            (
                172,
                29000.0,
                400.0,
                -45.0,
                0.0,
                150.0,
                150.0,
                4.0,
                2.207074681643701e-12,
                845.0483323181795,
            ),
            (
                172,
                29000.0,
                400.0,
                60.0,
                -70.0,
                70.0,
                70.0,
                4.0,
                4.566723510827348e-13,
                790.6239557119441,
            ),
            (
                172,
                29000.0,
                400.0,
                60.0,
                -70.0,
                250.0,
                250.0,
                4.0,
                6.957966569624449e-12,
                1349.462586622119,
            ),
            (
                172,
                29000.0,
                400.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                40.0,
                3.258252132785382e-12,
                1289.234960843807,
            ),
            (
                172,
                0.0,
                400.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                4.085403374565813e-12,
                1200.971191917643,
            ),
            (
                172,
                43200.0,
                400.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                3.230_763_252_933_96e-12,
                1165.060956208804,
            ),
            (
                1,
                29000.0,
                400.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                2.105972433040141e-12,
                843.1586610802336,
            ),
            (
                91,
                29000.0,
                400.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                3.180576942159127e-12,
                1030.982373989887,
            ),
            (
                274,
                29000.0,
                400.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                2.990266751727769e-12,
                988.6671231076829,
            ),
            (
                172,
                29000.0,
                50.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                1.295_541_393_193_92e-3,
                279.9138817220627,
            ),
            (
                172,
                29000.0,
                80.0,
                60.0,
                -70.0,
                150.0,
                150.0,
                4.0,
                2.583_834_857_191_45e-5,
                168.396_036_821_222,
            ),
        ];

        for (i, &(doy, sec, alt, lat, lon, f107a, f107, ap, ref_density, ref_temp)) in
            ref_values.iter().enumerate()
        {
            let out = call_model(alt, lat, lon, doy, sec, f107, f107a, ap);
            let density = out.d[5] * 1.0e3; // CGS g/cm3 -> kg/m3
            let temp = out.t[1];
            let density_rel_err = ((density - ref_density) / ref_density).abs();
            let temp_rel_err = ((temp - ref_temp) / ref_temp).abs();
            assert!(
                density_rel_err < 1.0e-10,
                "Case {i} (alt={alt}): density mismatch: got {density:.10e}, expected {ref_density:.10e}, rel_err={density_rel_err:.2e}"
            );
            assert!(
                temp_rel_err < 1.0e-10,
                "Case {i} (alt={alt}): temp mismatch: got {temp:.10}, expected {ref_temp:.10}, rel_err={temp_rel_err:.2e}"
            );
        }
    }
}