|
| 1 | +/* |
| 2 | + ISC License |
| 3 | +
|
| 4 | + Copyright (c) 2025, Department of Engineering Cybernetics, NTNU |
| 5 | +
|
| 6 | + Permission to use, copy, modify, and/or distribute this software for any |
| 7 | + purpose with or without fee is hereby granted, provided that the above |
| 8 | + copyright notice and this permission notice appear in all copies. |
| 9 | +
|
| 10 | + THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES |
| 11 | + WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF |
| 12 | + MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR |
| 13 | + ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES |
| 14 | + WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN |
| 15 | + ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF |
| 16 | + OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. |
| 17 | +
|
| 18 | + */ |
| 19 | + |
| 20 | +#include "architecture/utilities/ItuRefAtmosphere.h" |
| 21 | +#include <cmath> |
| 22 | +#include <algorithm> |
| 23 | + |
| 24 | +// Constructor |
| 25 | +ItuAtmosphere::ItuAtmosphere() { |
| 26 | +} |
| 27 | + |
| 28 | +// Destructor |
| 29 | +ItuAtmosphere::~ItuAtmosphere() { |
| 30 | +} |
| 31 | +const ItuAtmosphere& ItuAtmosphere::instance() { |
| 32 | + static ItuAtmosphere inst; |
| 33 | + return inst; |
| 34 | +} |
| 35 | + |
| 36 | +// Getters |
| 37 | +double ItuAtmosphere::getTempISA(double altitude_m) { |
| 38 | + return instance().calcTempISA(altitude_m); |
| 39 | +} |
| 40 | + |
| 41 | +double ItuAtmosphere::getPresISA(double altitude_m) { |
| 42 | + return instance().calcPresISA(altitude_m); |
| 43 | +} |
| 44 | + |
| 45 | +double ItuAtmosphere::getWaterVapDensityISA(double altitude_m) { |
| 46 | + return instance().calcWaterVapDensityISA(altitude_m); |
| 47 | +} |
| 48 | + |
| 49 | +size_t ItuAtmosphere::findLayerIndex(double altitude_m) const { |
| 50 | + size_t layer_idx = NUM_LAYERS - 1; |
| 51 | + for (size_t i = 0; i < NUM_LAYERS - 1; ++i) { |
| 52 | + if (altitude_m >= ATMOSPHERE_DATA[i].h && |
| 53 | + altitude_m < ATMOSPHERE_DATA[i + 1].h) |
| 54 | + { |
| 55 | + layer_idx = i; |
| 56 | + break; |
| 57 | + } |
| 58 | + } |
| 59 | + return layer_idx; |
| 60 | +} |
| 61 | + |
| 62 | +double ItuAtmosphere::geometricHeight(double altitude_m) const { |
| 63 | + // For altitudes above 84.852 km, convert geopotential height to geometric height |
| 64 | + return (R_EARTH_ITU * altitude_m * 1e-3) / (R_EARTH_ITU - altitude_m*1e-3); // ITU-R P.835-7 eq 1b |
| 65 | +} |
| 66 | + |
| 67 | +/* Calculate Temperature according to ITU-R P.835-7 |
| 68 | +*/ |
| 69 | +double ItuAtmosphere::calcTempISA(double altitude_m) const { |
| 70 | + // Clamp above model validity (85 km) |
| 71 | + if (altitude_m >= ATMOSPHERE_DATA[NUM_LAYERS - 1].h) { |
| 72 | + // if h above 84.852 km |
| 73 | + double Z = ItuAtmosphere::geometricHeight(altitude_m); // [km] Geometric height |
| 74 | + if (Z < 91.0) { |
| 75 | + return 186.8673; // [K] Constant temperature between 84.852 km and 91 km (ITU-R P.835-7 eq 4a) |
| 76 | + } else if (Z < 100.0) { |
| 77 | + return 263.1905 - 76.3232 * std::sqrt(1 - std::pow((Z - 91.0)/(19.9429),2)); // [K] Linear increase between 91 km and 100 km (ITU-R P.835-7 eq 4b) |
| 78 | + } else { |
| 79 | + return 195.08134; // [K] Above 100 km, assume constant temperature |
| 80 | + } |
| 81 | + } |
| 82 | + |
| 83 | + // Find which layer the altitude falls into (table in header file) |
| 84 | + size_t layer_idx = ItuAtmosphere::findLayerIndex(altitude_m); |
| 85 | + |
| 86 | + // Compute temperature at layer base: T_(layer_idx-1) |
| 87 | + double Ti = T0; |
| 88 | + for (size_t i = 0; i < layer_idx; ++i) { |
| 89 | + double h0 = ATMOSPHERE_DATA[i].h; // [m] Altitude at layer "i" base |
| 90 | + double h1 = ATMOSPHERE_DATA[i + 1].h; // [m] Altitude at layer "i+1" base |
| 91 | + double dT_dh_i = ATMOSPHERE_DATA[i].dT_dh; // [K/km] Temperature gradient in layer "i" |
| 92 | + Ti += dT_dh_i * (h1 - h0) * 1e-3; // [K] Temperature at top of layer "layer_idx-1" |
| 93 | + } |
| 94 | + |
| 95 | + // Linear variation within that layer |
| 96 | + double Hi = ATMOSPHERE_DATA[layer_idx].h; // [m] Altitude at layer "layer_idx" base |
| 97 | + double dT_dh = ATMOSPHERE_DATA[layer_idx].dT_dh; // [K/km] Temperature gradient in layer "layer_idx" |
| 98 | + |
| 99 | + return Ti + dT_dh * (altitude_m - Hi) * 1e-3; // [](ITU-R P.835-2 eq 2) |
| 100 | +} |
| 101 | + |
| 102 | +/* Calculate Pressure according to ITU-R P.835-7 |
| 103 | + * barometric formula for each atmospheric layer |
| 104 | + */ |
| 105 | +double ItuAtmosphere::calcPresISA(double altitude_m) const { |
| 106 | + // Find which layer the altitude falls into |
| 107 | + size_t layer_idx = ItuAtmosphere::findLayerIndex(altitude_m); |
| 108 | + double dT_dh = ATMOSPHERE_DATA[layer_idx].dT_dh; |
| 109 | + if (layer_idx == 0) { |
| 110 | + // 0 km < altitude < 11 km |
| 111 | + return 1013.25 * pow((288.15 / (288.15 - (dT_dh * altitude_m * 1e-3))), -itu_press_const/dT_dh); // ITU-R P.835-7 eq 3a |
| 112 | + } else if (layer_idx == 1) { |
| 113 | + // 11 km < altitude < 20 km |
| 114 | + return 226.3226 * exp(-0.157688 * (altitude_m * 1e-3 - 11.0)); // ITU-R P.835-7 eq 3b |
| 115 | + } else if (layer_idx == 2) { |
| 116 | + // 20 km < altitude < 32 km |
| 117 | + return 54.7498 * pow((216.65 / (216.65 + (altitude_m * 1e-3 - 20.0))), -itu_press_const); // ITU-R P.835-7 eq 3c |
| 118 | + } else if (layer_idx == 3) { |
| 119 | + // 32 km < altitude < 47 km |
| 120 | + return 8.680422 * pow((228.65 / (228.65 + (dT_dh * (altitude_m * 1e-3 - 32.0)))), -itu_press_const/dT_dh); // ITU-R P.835-7 eq 3d |
| 121 | + } else if (layer_idx == 4) { |
| 122 | + // 47 km < altitude < 51 km |
| 123 | + return 1.109106 * exp(-0.126226 * (altitude_m * 1e-3 - 47.0)); // ITU-R P.835-7 eq 3e |
| 124 | + } else if (layer_idx == 5) { |
| 125 | + // 51 km < altitude < 71 km |
| 126 | + return 0.6694167 * pow((270.65 / (270.65 - (dT_dh * (altitude_m * 1e-3 - 51.0)))), -itu_press_const/dT_dh); // ITU-R P.835-7 eq 3f |
| 127 | + } else if (layer_idx == 6) { |
| 128 | + // 71 km < altitude < 85 km |
| 129 | + return 0.03956649 * pow((214.65 / (214.65 - (dT_dh * (altitude_m * 1e-3 - 71.0)))), -itu_press_const/dT_dh); // ITU-R P.835-7 eq 3g |
| 130 | + } else { |
| 131 | + // altitude >= 85 km |
| 132 | + double Z = ItuAtmosphere::geometricHeight(altitude_m); // [km] Geometric height |
| 133 | + if (Z < 100.0e3) { |
| 134 | + return exp(95.571899 - 4.011801*Z + 6.424731*Z*Z*1e-2 - 4.789660*Z*Z*Z*1e-4 + 1.340543*Z*Z*Z*Z*1e-6); // ITU-R P.835-7 eq 5 |
| 135 | + } else { |
| 136 | + return 0.0; // Pressure above 100 km is assumed zero |
| 137 | + } |
| 138 | + } |
| 139 | +} |
| 140 | + |
| 141 | +/* Calculate Water Vapor Density according to ITU-R P.835-7 |
| 142 | + * Uses exponential decay model with different scale heights |
| 143 | + */ |
| 144 | +double ItuAtmosphere::calcWaterVapDensityISA(double altitude_m) const { |
| 145 | + // Calculate geometric height |
| 146 | + double Z = ItuAtmosphere::geometricHeight(altitude_m); // [km] Geometric height |
| 147 | + double rho = rho0 * std::exp(-Z / h0_water); // [g/m3] ITU-R P.835-2 eq 6 |
| 148 | + double T_Z = ItuAtmosphere::calcTempISA(altitude_m); // [K] |
| 149 | + double e_z = (rho*T_Z)/(216.7); // [hPa] Partial pressure at height Z |
| 150 | + double p_z = ItuAtmosphere::calcPresISA(altitude_m); // [hPa] Total pressure at height Z |
| 151 | + if (e_z/p_z < 2.0e-6) { |
| 152 | + // Use ITU-R 835-7 eq. 8 |
| 153 | + rho = 2.0e-6 * (p_z * 216.7)/(T_Z); // [g/m3] |
| 154 | + } |
| 155 | + return rho; |
| 156 | +} |
0 commit comments