//! 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}"
);
}
}
}