| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // This file is part of INSTINCT, the INS Toolkit for Integrated | ||
| 2 | // Navigation Concepts and Training by the Institute of Navigation of | ||
| 3 | // the University of Stuttgart, Germany. | ||
| 4 | // | ||
| 5 | // This Source Code Form is subject to the terms of the Mozilla Public | ||
| 6 | // License, v. 2.0. If a copy of the MPL was not distributed with this | ||
| 7 | // file, You can obtain one at https://mozilla.org/MPL/2.0/. | ||
| 8 | |||
| 9 | #include "GPT.hpp" | ||
| 10 | |||
| 11 | #include <array> | ||
| 12 | #include <cmath> | ||
| 13 | |||
| 14 | #include "internal/GPT2Coeffs.hpp" | ||
| 15 | #include "internal/GPT3Coeffs.hpp" | ||
| 16 | |||
| 17 | #include "Navigation/Constants.hpp" | ||
| 18 | #include "Navigation/Math/Math.hpp" | ||
| 19 | |||
| 20 | namespace NAV | ||
| 21 | { | ||
| 22 | 725 | GPT2output GPT2_param(const double& mjd, const Eigen::Vector3d& lla_pos) | |
| 23 | { | ||
| 24 | using internal::GPT2_grid; | ||
| 25 | |||
| 26 | 725 | GPT2output gpt2outputs; | |
| 27 | |||
| 28 | // change the reference epoch to January 1 2000 | ||
| 29 | 725 | double mjd1 = mjd - 51544.5; | |
| 30 | // GPT2 with time variation | ||
| 31 | 725 | double cosfy = cos(mjd1 / 365.25 * 2.0 * M_PI); | |
| 32 | 725 | double coshy = cos(mjd1 / 365.25 * 4.0 * M_PI); | |
| 33 | 725 | double sinfy = sin(mjd1 / 365.25 * 2.0 * M_PI); | |
| 34 | 725 | double sinhy = sin(mjd1 / 365.25 * 4.0 * M_PI); | |
| 35 | |||
| 36 |
1/2✓ Branch 1 taken 725 times.
✗ Branch 2 not taken.
|
725 | double dlat = lla_pos(0); |
| 37 |
1/2✓ Branch 1 taken 725 times.
✗ Branch 2 not taken.
|
725 | double dlon = lla_pos(1); |
| 38 |
1/2✓ Branch 1 taken 725 times.
✗ Branch 2 not taken.
|
725 | double hell = lla_pos(2); |
| 39 | |||
| 40 | // only positive longitude in degrees | ||
| 41 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 725 times.
|
725 | double plon = dlon < 0.0 ? (dlon + 2.0 * M_PI) * 180.0 / M_PI : dlon * 180.0 / M_PI; |
| 42 | |||
| 43 | // transform to polar distance in degrees | ||
| 44 | 725 | double ppod = (-dlat + M_PI / 2.0) * 180.0 / M_PI; | |
| 45 | |||
| 46 | // find the index (line in the grid file) of the nearest point | ||
| 47 | 725 | auto ipod = static_cast<int>(floor(ppod + 1.0)); | |
| 48 | 725 | auto ilon = static_cast<int>(floor(plon + 1.0)); | |
| 49 | |||
| 50 | // normalized (to one) differences, can be positive or negative | ||
| 51 | 725 | double diffpod = ppod - (ipod - 0.5); | |
| 52 | 725 | double difflon = plon - (ilon - 0.5); | |
| 53 | |||
| 54 | // added by HCY | ||
| 55 |
2/2✓ Branch 0 taken 36 times.
✓ Branch 1 taken 689 times.
|
725 | if (ipod == 181) { ipod = 180; } |
| 56 | // added by GP | ||
| 57 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 725 times.
|
725 | if (ilon == 361) { ilon = 1; } |
| 58 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 725 times.
|
725 | if (ilon == 0) { ilon = 360; } |
| 59 | |||
| 60 | 725 | constexpr size_t N = 4; | |
| 61 | 725 | std::array<size_t, N> indx{}; | |
| 62 | // get the number of the corresponding line | ||
| 63 | 725 | indx[0] = static_cast<size_t>((ipod - 1) * 360 + ilon) - 1; | |
| 64 | |||
| 65 | // near the poles: nearest neighbour interpolation | ||
| 66 | // otherwise: bilinear, with the 1 degree grid the limits are lower and upper (GP) | ||
| 67 |
4/4✓ Branch 0 taken 689 times.
✓ Branch 1 taken 36 times.
✓ Branch 2 taken 36 times.
✓ Branch 3 taken 653 times.
|
725 | if (ppod <= 0.5 || ppod >= 179.5) // case of nearest neighborhood |
| 68 | { | ||
| 69 | 72 | auto ix = indx[0]; | |
| 70 | // transforming ellipsoidal height to orthometric height | ||
| 71 |
1/2✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
|
72 | gpt2outputs.undu = GPT2_grid.at(ix).undulation; // undu: geoid undulation in m |
| 72 | 72 | double hgt = hell - gpt2outputs.undu; | |
| 73 | |||
| 74 | // pressure, temperature at the height of the grid | ||
| 75 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | double T0 = GPT2_grid.at(ix).T_A0 + GPT2_grid.at(ix).T_A1 * cosfy + GPT2_grid.at(ix).T_B1 * sinfy + GPT2_grid.at(ix).T_A2 * coshy + GPT2_grid.at(ix).T_B2 * sinhy; |
| 76 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | double p0 = GPT2_grid.at(ix).p_A0 + GPT2_grid.at(ix).p_A1 * cosfy + GPT2_grid.at(ix).p_B1 * sinfy + GPT2_grid.at(ix).p_A2 * coshy + GPT2_grid.at(ix).p_B2 * sinhy; |
| 77 | |||
| 78 | // specific humidity | ||
| 79 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | double Q = (GPT2_grid.at(ix).Q_A0 + GPT2_grid.at(ix).Q_A1 * cosfy + GPT2_grid.at(ix).Q_B1 * sinfy + GPT2_grid.at(ix).Q_A2 * coshy + GPT2_grid.at(ix).Q_B2 * sinhy) |
| 80 | 72 | / 1e3; | |
| 81 | |||
| 82 | // lapse rate of the temperature | ||
| 83 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | gpt2outputs.dT = (GPT2_grid.at(ix).dT_A0 + GPT2_grid.at(ix).dT_A1 * cosfy + GPT2_grid.at(ix).dT_B1 * sinfy + GPT2_grid.at(ix).dT_A2 * coshy + GPT2_grid.at(ix).dT_B2 * sinhy) |
| 84 | 72 | / 1e3; | |
| 85 | |||
| 86 | // station height - grid height | ||
| 87 |
1/2✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
|
72 | double redh = hgt - GPT2_grid.at(ix).Hs; // Hs: orthometric grid height in m |
| 88 | |||
| 89 | // temperature at station height in Celsius | ||
| 90 | 72 | gpt2outputs.T = T0 + gpt2outputs.dT * redh - 273.150; | |
| 91 | |||
| 92 | // temperature lapse rate in degrees / km | ||
| 93 | 72 | gpt2outputs.dT = gpt2outputs.dT * 1e3; | |
| 94 | |||
| 95 | // virtual temperature in Kelvin | ||
| 96 | 72 | double Tv = T0 * (1 + 0.6077 * Q); | |
| 97 | |||
| 98 | 72 | double c = InsConst::G_NORM * InsConst::dMtr / (InsConst::Rg * Tv); | |
| 99 | |||
| 100 | // pressure in hPa | ||
| 101 | 72 | gpt2outputs.p = (p0 * exp(-c * redh)) / 100.0; | |
| 102 | |||
| 103 | // hydrostatic coefficient ah | ||
| 104 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | gpt2outputs.ah = (GPT2_grid.at(ix).h_A0 + GPT2_grid.at(ix).h_A1 * cosfy + GPT2_grid.at(ix).h_B1 * sinfy + GPT2_grid.at(ix).h_A2 * coshy + GPT2_grid.at(ix).h_B2 * sinhy) |
| 105 | 72 | / 1e3; | |
| 106 | |||
| 107 | // wet coefficient aw | ||
| 108 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | gpt2outputs.aw = (GPT2_grid.at(ix).w_A0 + GPT2_grid.at(ix).w_A1 * cosfy + GPT2_grid.at(ix).w_B1 * sinfy + GPT2_grid.at(ix).w_A2 * coshy + GPT2_grid.at(ix).w_B2 * sinhy) |
| 109 | 72 | / 1e3; | |
| 110 | |||
| 111 | // water vapour decrease factor la - added by GP | ||
| 112 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | gpt2outputs.la = GPT2_grid.at(ix).la_A0 + GPT2_grid.at(ix).la_A1 * cosfy + GPT2_grid.at(ix).la_B1 * sinfy + GPT2_grid.at(ix).la_A2 * coshy + GPT2_grid.at(ix).la_B2 * sinhy; |
| 113 | |||
| 114 | // mean temperature of the water vapor Tm - added by GP | ||
| 115 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | gpt2outputs.Tm = GPT2_grid.at(ix).Tm_A0 + GPT2_grid.at(ix).Tm_A1 * cosfy + GPT2_grid.at(ix).Tm_B1 * sinfy + GPT2_grid.at(ix).Tm_A2 * coshy + GPT2_grid.at(ix).Tm_B2 * sinhy; |
| 116 | |||
| 117 | // water vapor pressure in hPa - changed by GP | ||
| 118 | 72 | double e0 = Q * p0 / (0.622 + 0.378 * Q) / 100.0; // on the grid | |
| 119 | |||
| 120 | 72 | gpt2outputs.e = e0 * pow((100.0 * gpt2outputs.p / p0), (gpt2outputs.la + 1)); // on the station height - (14) Askne and Nordius, 1987 | |
| 121 | 72 | } | |
| 122 | else // bilinear interpolation | ||
| 123 | { | ||
| 124 | 653 | double ipod1 = ipod + int(math::sign(1.0, diffpod)); | |
| 125 | 653 | double ilon1 = ilon + int(math::sign(1.0, difflon)); | |
| 126 | |||
| 127 | // changed for the 1 degree grid (GP) | ||
| 128 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 653 times.
|
653 | if (ilon1 == 361) { ilon1 = 1; } |
| 129 |
2/2✓ Branch 0 taken 17 times.
✓ Branch 1 taken 636 times.
|
653 | if (ilon1 == 0) { ilon1 = 360; } |
| 130 | // get the number of the line | ||
| 131 | 653 | indx[1] = static_cast<size_t>((ipod1 - 1) * 360 + ilon) - 1; // along same longitude | |
| 132 | 653 | indx[2] = static_cast<size_t>((ipod - 1) * 360 + ilon1) - 1; // along same polar distance | |
| 133 | 653 | indx[3] = static_cast<size_t>((ipod1 - 1) * 360 + ilon1) - 1; // diagonal | |
| 134 | |||
| 135 | 653 | std::array<double, N> undul{}; | |
| 136 | 653 | std::array<double, N> Ql{}; | |
| 137 | 653 | std::array<double, N> dTl{}; | |
| 138 | 653 | std::array<double, N> Tl{}; | |
| 139 | 653 | std::array<double, N> pl{}; | |
| 140 | 653 | std::array<double, N> ahl{}; | |
| 141 | 653 | std::array<double, N> awl{}; | |
| 142 | 653 | std::array<double, N> lal{}; | |
| 143 | 653 | std::array<double, N> Tml{}; | |
| 144 | 653 | std::array<double, N> el{}; | |
| 145 | |||
| 146 |
2/2✓ Branch 0 taken 2612 times.
✓ Branch 1 taken 653 times.
|
3265 | for (size_t l = 0; l < N; l++) |
| 147 | { | ||
| 148 | // transforming ellipsoidal height to orthometric height: | ||
| 149 | // Hortho = -N + Hell | ||
| 150 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | auto ix = indx.at(l); |
| 151 | |||
| 152 |
2/4✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
|
2612 | undul.at(l) = GPT2_grid.at(ix).undulation; // undu: geoid undulation in m |
| 153 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | double hgt = hell - undul.at(l); |
| 154 | |||
| 155 | // pressure, temperature at the height of the grid | ||
| 156 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | double T0 = GPT2_grid.at(ix).T_A0 + GPT2_grid.at(ix).T_A1 * cosfy + GPT2_grid.at(ix).T_B1 * sinfy + GPT2_grid.at(ix).T_A2 * coshy + GPT2_grid.at(ix).T_B2 * sinhy; |
| 157 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | double p0 = GPT2_grid.at(ix).p_A0 + GPT2_grid.at(ix).p_A1 * cosfy + GPT2_grid.at(ix).p_B1 * sinfy + GPT2_grid.at(ix).p_A2 * coshy + GPT2_grid.at(ix).p_B2 * sinhy; |
| 158 | |||
| 159 | // humidity | ||
| 160 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | Ql.at(l) = (GPT2_grid.at(ix).Q_A0 + GPT2_grid.at(ix).Q_A1 * cosfy + GPT2_grid.at(ix).Q_B1 * sinfy + GPT2_grid.at(ix).Q_A2 * coshy + GPT2_grid.at(ix).Q_B2 * sinhy) |
| 161 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | / 1e3; |
| 162 | |||
| 163 | // reduction = stationheight - gridheight | ||
| 164 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | double redh = hgt - GPT2_grid.at(ix).Hs; // Hs: orthometric grid height in m |
| 165 | |||
| 166 | // lapse rate of the temperature in degree / m | ||
| 167 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | dTl.at(l) = (GPT2_grid.at(ix).dT_A0 + GPT2_grid.at(ix).dT_A1 * cosfy + GPT2_grid.at(ix).dT_B1 * sinfy + GPT2_grid.at(ix).dT_A2 * coshy + GPT2_grid.at(ix).dT_B2 * sinhy) |
| 168 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | / 1e3; |
| 169 | |||
| 170 | // temperature reduction to station height | ||
| 171 |
2/4✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
|
2612 | Tl.at(l) = T0 + dTl.at(l) * redh - 273.150; |
| 172 | |||
| 173 | // virtual temperature | ||
| 174 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | double Tv = T0 * (1 + 0.6077 * Ql.at(l)); |
| 175 | 2612 | double c = InsConst::G_NORM * InsConst::dMtr / (InsConst::Rg * Tv); | |
| 176 | |||
| 177 | // pressure in hPa | ||
| 178 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | pl.at(l) = (p0 * exp(-c * redh)) / 100.0; |
| 179 | |||
| 180 | // hydrostatic coefficient ah | ||
| 181 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | ahl.at(l) = (GPT2_grid.at(ix).h_A0 + GPT2_grid.at(ix).h_A1 * cosfy + GPT2_grid.at(ix).h_B1 * sinfy + GPT2_grid.at(ix).h_A2 * coshy + GPT2_grid.at(ix).h_B2 * sinhy) |
| 182 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | / 1e3; |
| 183 | |||
| 184 | // wet coefficient aw | ||
| 185 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | awl.at(l) = (GPT2_grid.at(ix).w_A0 + GPT2_grid.at(ix).w_A1 * cosfy + GPT2_grid.at(ix).w_B1 * sinfy + GPT2_grid.at(ix).w_A2 * coshy + GPT2_grid.at(ix).w_B2 * sinhy) |
| 186 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | / 1e3; |
| 187 | |||
| 188 | // water vapor decrease factor la - added by GP | ||
| 189 |
6/12✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2612 times.
✗ Branch 17 not taken.
|
2612 | lal.at(l) = GPT2_grid.at(ix).la_A0 + GPT2_grid.at(ix).la_A1 * cosfy + GPT2_grid.at(ix).la_B1 * sinfy + GPT2_grid.at(ix).la_A2 * coshy + GPT2_grid.at(ix).la_B2 * sinhy; |
| 190 | |||
| 191 | // mean temperature of the water vapor Tm - added by GP | ||
| 192 |
6/12✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2612 times.
✗ Branch 17 not taken.
|
2612 | Tml.at(l) = GPT2_grid.at(ix).Tm_A0 + GPT2_grid.at(ix).Tm_A1 * cosfy + GPT2_grid.at(ix).Tm_B1 * sinfy + GPT2_grid.at(ix).Tm_A2 * coshy + GPT2_grid.at(ix).Tm_B2 * sinhy; |
| 193 | |||
| 194 | // water vapor pressure in hPa - changed by GP | ||
| 195 |
2/4✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
|
2612 | double e0 = Ql.at(l) * p0 / (0.622 + 0.378 * Ql.at(l)) / 100.0; // on the grid |
| 196 | |||
| 197 |
3/6✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
|
2612 | el.at(l) = e0 * pow((100.0 * pl.at(l) / p0), (lal.at(l) + 1.0)); // on the station height (14) Askne and Nordius, 1987 |
| 198 | } | ||
| 199 | |||
| 200 | 653 | double dnpod1 = fabs(diffpod); // distance nearer point | |
| 201 | 653 | double dnpod2 = 1.0 - dnpod1; // distance to distant point | |
| 202 | 653 | double dnlon1 = fabs(difflon); | |
| 203 | 653 | double dnlon2 = 1.0 - dnlon1; | |
| 204 | |||
| 205 | // pressure | ||
| 206 | 653 | double R1 = dnpod2 * pl[0] + dnpod1 * pl[1]; | |
| 207 | 653 | double R2 = dnpod2 * pl[2] + dnpod1 * pl[3]; | |
| 208 | 653 | gpt2outputs.p = dnlon2 * R1 + dnlon1 * R2; | |
| 209 | |||
| 210 | // temperature | ||
| 211 | 653 | R1 = dnpod2 * Tl[0] + dnpod1 * Tl[1]; | |
| 212 | 653 | R2 = dnpod2 * Tl[2] + dnpod1 * Tl[3]; | |
| 213 | 653 | gpt2outputs.T = dnlon2 * R1 + dnlon1 * R2; | |
| 214 | |||
| 215 | // temperature in degree per km | ||
| 216 | 653 | R1 = dnpod2 * dTl[0] + dnpod1 * dTl[1]; | |
| 217 | 653 | R2 = dnpod2 * dTl[2] + dnpod1 * dTl[3]; | |
| 218 | 653 | gpt2outputs.dT = (dnlon2 * R1 + dnlon1 * R2) * 1e3; | |
| 219 | |||
| 220 | // water vapor pressure in hPa - changed by GP | ||
| 221 | 653 | R1 = dnpod2 * el[0] + dnpod1 * el[1]; | |
| 222 | 653 | R2 = dnpod2 * el[2] + dnpod1 * el[3]; | |
| 223 | 653 | gpt2outputs.e = dnlon2 * R1 + dnlon1 * R2; | |
| 224 | |||
| 225 | // hydrostatic | ||
| 226 | 653 | R1 = dnpod2 * ahl[0] + dnpod1 * ahl[1]; | |
| 227 | 653 | R2 = dnpod2 * ahl[2] + dnpod1 * ahl[3]; | |
| 228 | 653 | gpt2outputs.ah = dnlon2 * R1 + dnlon1 * R2; | |
| 229 | |||
| 230 | // wet | ||
| 231 | 653 | R1 = dnpod2 * awl[0] + dnpod1 * awl[1]; | |
| 232 | 653 | R2 = dnpod2 * awl[2] + dnpod1 * awl[3]; | |
| 233 | 653 | gpt2outputs.aw = dnlon2 * R1 + dnlon1 * R2; | |
| 234 | |||
| 235 | // undulation | ||
| 236 | 653 | R1 = dnpod2 * undul[0] + dnpod1 * undul[1]; | |
| 237 | 653 | R2 = dnpod2 * undul[2] + dnpod1 * undul[3]; | |
| 238 | 653 | gpt2outputs.undu = dnlon2 * R1 + dnlon1 * R2; | |
| 239 | |||
| 240 | // water vapor decrease factor la - added by GP | ||
| 241 | 653 | R1 = dnpod2 * lal[0] + dnpod1 * lal[1]; | |
| 242 | 653 | R2 = dnpod2 * lal[2] + dnpod1 * lal[3]; | |
| 243 | 653 | gpt2outputs.la = dnlon2 * R1 + dnlon1 * R2; | |
| 244 | |||
| 245 | // mean temperature of the water vapor - added by GP | ||
| 246 | 653 | R1 = dnpod2 * Tml[0] + dnpod1 * Tml[1]; | |
| 247 | 653 | R2 = dnpod2 * Tml[2] + dnpod1 * Tml[3]; | |
| 248 | 653 | gpt2outputs.Tm = dnlon2 * R1 + dnlon1 * R2; | |
| 249 | } | ||
| 250 | |||
| 251 | 1450 | return gpt2outputs; | |
| 252 | } | ||
| 253 | |||
| 254 | 725 | GPT3output GPT3_param(const double& mjd, const Eigen::Vector3d& lla_pos) | |
| 255 | { | ||
| 256 | using internal::GPT3_grid; | ||
| 257 | |||
| 258 | 725 | GPT3output gpt3outputs; | |
| 259 | |||
| 260 | // convert mjd to doy | ||
| 261 |
1/2✓ Branch 1 taken 725 times.
✗ Branch 2 not taken.
|
725 | double doy = mjd2doy(mjd); |
| 262 | // factors for amplitudes | ||
| 263 | 725 | double cosfy = cos(doy / 365.25 * 2.0 * M_PI); | |
| 264 | 725 | double coshy = cos(doy / 365.25 * 4.0 * M_PI); | |
| 265 | 725 | double sinfy = sin(doy / 365.25 * 2.0 * M_PI); | |
| 266 | 725 | double sinhy = sin(doy / 365.25 * 4.0 * M_PI); | |
| 267 | |||
| 268 |
1/2✓ Branch 1 taken 725 times.
✗ Branch 2 not taken.
|
725 | double dlat = lla_pos(0); |
| 269 |
1/2✓ Branch 1 taken 725 times.
✗ Branch 2 not taken.
|
725 | double dlon = lla_pos(1); |
| 270 |
1/2✓ Branch 1 taken 725 times.
✗ Branch 2 not taken.
|
725 | double hell = lla_pos(2); |
| 271 | |||
| 272 | // only positive longitude in degrees | ||
| 273 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 725 times.
|
725 | double plon = dlon < 0.0 ? (dlon + 2.0 * M_PI) * 180.0 / M_PI : dlon * 180.0 / M_PI; |
| 274 | |||
| 275 | // transform to polar distance in degrees | ||
| 276 | 725 | double ppod = (-dlat + M_PI / 2.0) * 180.0 / M_PI; | |
| 277 | |||
| 278 | // find the index (line in the grid file) of the nearest point | ||
| 279 | 725 | auto ipod = static_cast<int>(floor(ppod + 1.0)); | |
| 280 | 725 | auto ilon = static_cast<int>(floor(plon + 1.0)); | |
| 281 | |||
| 282 | // normalized (to one) differences, can be positive or negative | ||
| 283 | 725 | double diffpod = ppod - (ipod - 0.5); | |
| 284 | 725 | double difflon = plon - (ilon - 0.5); | |
| 285 | |||
| 286 | // added by HCY | ||
| 287 |
2/2✓ Branch 0 taken 36 times.
✓ Branch 1 taken 689 times.
|
725 | if (ipod == 181) { ipod = 180; } |
| 288 | // added by GP | ||
| 289 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 725 times.
|
725 | if (ilon == 361) { ilon = 1; } |
| 290 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 725 times.
|
725 | if (ilon == 0) { ilon = 360; } |
| 291 | |||
| 292 | 725 | constexpr size_t N = 4; | |
| 293 | 725 | std::array<size_t, N> indx{}; | |
| 294 | // get the number of the corresponding line | ||
| 295 | 725 | indx[0] = static_cast<size_t>((ipod - 1) * 360 + ilon) - 1; | |
| 296 | |||
| 297 | // near the poles: nearest neighbour interpolation | ||
| 298 | // otherwise: bilinear, with the 1 degree grid the limits are lower and upper (GP) | ||
| 299 |
4/4✓ Branch 0 taken 689 times.
✓ Branch 1 taken 36 times.
✓ Branch 2 taken 36 times.
✓ Branch 3 taken 653 times.
|
725 | if (ppod <= 0.5 || ppod >= 179.5) // case of nearest neighborhood |
| 300 | { | ||
| 301 | 72 | auto ix = indx[0]; | |
| 302 | // transforming ellipsoidal height to orthometric height | ||
| 303 |
1/2✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
|
72 | gpt3outputs.undu = GPT3_grid.at(ix).undulation; // undu: geoid undulation in m |
| 304 | 72 | double hgt = hell - gpt3outputs.undu; | |
| 305 | |||
| 306 | // pressure, temperature at the height of the grid | ||
| 307 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | double T0 = GPT3_grid.at(ix).T_A0 + GPT3_grid.at(ix).T_A1 * cosfy + GPT3_grid.at(ix).T_B1 * sinfy + GPT3_grid.at(ix).T_A2 * coshy + GPT3_grid.at(ix).T_B2 * sinhy; |
| 308 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | double p0 = GPT3_grid.at(ix).p_A0 + GPT3_grid.at(ix).p_A1 * cosfy + GPT3_grid.at(ix).p_B1 * sinfy + GPT3_grid.at(ix).p_A2 * coshy + GPT3_grid.at(ix).p_B2 * sinhy; |
| 309 | |||
| 310 | // specific humidity | ||
| 311 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | double Q = (GPT3_grid.at(ix).Q_A0 + GPT3_grid.at(ix).Q_A1 * cosfy + GPT3_grid.at(ix).Q_B1 * sinfy + GPT3_grid.at(ix).Q_A2 * coshy + GPT3_grid.at(ix).Q_B2 * sinhy) |
| 312 | 72 | / 1e3; | |
| 313 | |||
| 314 | // lapse rate of the temperature | ||
| 315 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | gpt3outputs.dT = (GPT3_grid.at(ix).dT_A0 + GPT3_grid.at(ix).dT_A1 * cosfy + GPT3_grid.at(ix).dT_B1 * sinfy + GPT3_grid.at(ix).dT_A2 * coshy + GPT3_grid.at(ix).dT_B2 * sinhy) |
| 316 | 72 | / 1e3; | |
| 317 | |||
| 318 | // station height - grid height | ||
| 319 |
1/2✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
|
72 | double redh = hgt - GPT3_grid.at(ix).Hs; // Hs: orthometric grid height in m |
| 320 | |||
| 321 | // temperature at station height in Celsius | ||
| 322 | 72 | gpt3outputs.T = T0 + gpt3outputs.dT * redh - 273.150; | |
| 323 | |||
| 324 | // temperature lapse rate in degrees / km | ||
| 325 | 72 | gpt3outputs.dT = gpt3outputs.dT * 1e3; | |
| 326 | |||
| 327 | // virtual temperature in Kelvin | ||
| 328 | 72 | double Tv = T0 * (1 + 0.6077 * Q); | |
| 329 | |||
| 330 | 72 | double c = InsConst::G_NORM * InsConst::dMtr / (InsConst::Rg * Tv); | |
| 331 | |||
| 332 | // pressure in hPa | ||
| 333 | 72 | gpt3outputs.p = (p0 * exp(-c * redh)) / 100.0; | |
| 334 | |||
| 335 | // hydrostatic coefficient ah | ||
| 336 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | gpt3outputs.ah = (GPT3_grid.at(ix).h_A0 + GPT3_grid.at(ix).h_A1 * cosfy + GPT3_grid.at(ix).h_B1 * sinfy + GPT3_grid.at(ix).h_A2 * coshy + GPT3_grid.at(ix).h_B2 * sinhy) |
| 337 | 72 | / 1e3; | |
| 338 | |||
| 339 | // wet coefficient aw | ||
| 340 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | gpt3outputs.aw = (GPT3_grid.at(ix).w_A0 + GPT3_grid.at(ix).w_A1 * cosfy + GPT3_grid.at(ix).w_B1 * sinfy + GPT3_grid.at(ix).w_A2 * coshy + GPT3_grid.at(ix).w_B2 * sinhy) |
| 341 | 72 | / 1e3; | |
| 342 | |||
| 343 | // water vapour decrease factor la - added by GP | ||
| 344 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | gpt3outputs.la = GPT3_grid.at(ix).la_A0 + GPT3_grid.at(ix).la_A1 * cosfy + GPT3_grid.at(ix).la_B1 * sinfy + GPT3_grid.at(ix).la_A2 * coshy + GPT3_grid.at(ix).la_B2 * sinhy; |
| 345 | |||
| 346 | // mean temperature of the water vapor Tm - added by GP | ||
| 347 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | gpt3outputs.Tm = GPT3_grid.at(ix).Tm_A0 + GPT3_grid.at(ix).Tm_A1 * cosfy + GPT3_grid.at(ix).Tm_B1 * sinfy + GPT3_grid.at(ix).Tm_A2 * coshy + GPT3_grid.at(ix).Tm_B2 * sinhy; |
| 348 | |||
| 349 | // water vapor pressure in hPa - changed by GP | ||
| 350 | 72 | double e0 = Q * p0 / (0.622 + 0.378 * Q) / 100.0; // on the grid | |
| 351 | |||
| 352 | 72 | gpt3outputs.e = e0 * pow((100.0 * gpt3outputs.p / p0), (gpt3outputs.la + 1)); // on the station height - (14) Askne and Nordius, 1987 | |
| 353 | |||
| 354 | // north and east gradients (total, hydrostatic and wet) | ||
| 355 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | gpt3outputs.Gn_h = (GPT3_grid.at(ix).Gnh_A0 + GPT3_grid.at(ix).Gnh_A1 * cosfy + GPT3_grid.at(ix).Gnh_B1 * sinfy + GPT3_grid.at(ix).Gnh_A2 * coshy + GPT3_grid.at(ix).Gnh_B2 * sinhy) |
| 356 | 72 | / 1e5; | |
| 357 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | gpt3outputs.Ge_h = (GPT3_grid.at(ix).Geh_A0 + GPT3_grid.at(ix).Geh_A1 * cosfy + GPT3_grid.at(ix).Geh_B1 * sinfy + GPT3_grid.at(ix).Geh_A2 * coshy + GPT3_grid.at(ix).Geh_B2 * sinhy) |
| 358 | 72 | / 1e5; | |
| 359 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | gpt3outputs.Gn_w = (GPT3_grid.at(ix).Gnw_A0 + GPT3_grid.at(ix).Gnw_A1 * cosfy + GPT3_grid.at(ix).Gnw_B1 * sinfy + GPT3_grid.at(ix).Gnw_A2 * coshy + GPT3_grid.at(ix).Gnw_B2 * sinhy) |
| 360 | 72 | / 1e5; | |
| 361 |
5/10✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 72 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 72 times.
✗ Branch 14 not taken.
|
72 | gpt3outputs.Ge_w = (GPT3_grid.at(ix).Gew_A0 + GPT3_grid.at(ix).Gew_A1 * cosfy + GPT3_grid.at(ix).Gew_B1 * sinfy + GPT3_grid.at(ix).Gew_A2 * coshy + GPT3_grid.at(ix).Gew_B2 * sinhy) |
| 362 | 72 | / 1e5; | |
| 363 | 72 | } | |
| 364 | else // bilinear interpolation | ||
| 365 | { | ||
| 366 | 653 | double ipod1 = ipod + int(math::sign(1.0, diffpod)); | |
| 367 | 653 | double ilon1 = ilon + int(math::sign(1.0, difflon)); | |
| 368 | |||
| 369 | // changed for the 1 degree grid (GP) | ||
| 370 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 653 times.
|
653 | if (ilon1 == 361) { ilon1 = 1; } |
| 371 |
2/2✓ Branch 0 taken 17 times.
✓ Branch 1 taken 636 times.
|
653 | if (ilon1 == 0) { ilon1 = 360; } |
| 372 | // get the number of the line | ||
| 373 | 653 | indx[1] = static_cast<size_t>((ipod1 - 1) * 360 + ilon) - 1; // along same longitude | |
| 374 | 653 | indx[2] = static_cast<size_t>((ipod - 1) * 360 + ilon1) - 1; // along same polar distance | |
| 375 | 653 | indx[3] = static_cast<size_t>((ipod1 - 1) * 360 + ilon1) - 1; // diagonal | |
| 376 | |||
| 377 | 653 | std::array<double, N> undul{}; | |
| 378 | 653 | std::array<double, N> Ql{}; | |
| 379 | 653 | std::array<double, N> dTl{}; | |
| 380 | 653 | std::array<double, N> Tl{}; | |
| 381 | 653 | std::array<double, N> pl{}; | |
| 382 | 653 | std::array<double, N> ahl{}; | |
| 383 | 653 | std::array<double, N> awl{}; | |
| 384 | 653 | std::array<double, N> lal{}; | |
| 385 | 653 | std::array<double, N> Tml{}; | |
| 386 | 653 | std::array<double, N> el{}; | |
| 387 | 653 | std::array<double, N> Gn_hl{}; | |
| 388 | 653 | std::array<double, N> Ge_hl{}; | |
| 389 | 653 | std::array<double, N> Gn_wl{}; | |
| 390 | 653 | std::array<double, N> Ge_wl{}; | |
| 391 | |||
| 392 |
2/2✓ Branch 0 taken 2612 times.
✓ Branch 1 taken 653 times.
|
3265 | for (size_t l = 0; l < N; l++) |
| 393 | { | ||
| 394 | // transforming ellipsoidal height to orthometric height: | ||
| 395 | // Hortho = -N + Hell | ||
| 396 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | auto ix = indx.at(l); |
| 397 | |||
| 398 |
2/4✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
|
2612 | undul.at(l) = GPT3_grid.at(ix).undulation; // undu: geoid undulation in m |
| 399 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | double hgt = hell - undul.at(l); |
| 400 | |||
| 401 | // pressure, temperature at the height of the grid | ||
| 402 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | double T0 = GPT3_grid.at(ix).T_A0 + GPT3_grid.at(ix).T_A1 * cosfy + GPT3_grid.at(ix).T_B1 * sinfy + GPT3_grid.at(ix).T_A2 * coshy + GPT3_grid.at(ix).T_B2 * sinhy; |
| 403 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | double p0 = GPT3_grid.at(ix).p_A0 + GPT3_grid.at(ix).p_A1 * cosfy + GPT3_grid.at(ix).p_B1 * sinfy + GPT3_grid.at(ix).p_A2 * coshy + GPT3_grid.at(ix).p_B2 * sinhy; |
| 404 | |||
| 405 | // humidity | ||
| 406 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | Ql.at(l) = (GPT3_grid.at(ix).Q_A0 + GPT3_grid.at(ix).Q_A1 * cosfy + GPT3_grid.at(ix).Q_B1 * sinfy + GPT3_grid.at(ix).Q_A2 * coshy + GPT3_grid.at(ix).Q_B2 * sinhy) |
| 407 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | / 1e3; |
| 408 | |||
| 409 | // reduction = stationheight - gridheight | ||
| 410 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | double redh = hgt - GPT3_grid.at(ix).Hs; // Hs: orthometric grid height in m |
| 411 | |||
| 412 | // lapse rate of the temperature in degree / m | ||
| 413 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | dTl.at(l) = (GPT3_grid.at(ix).dT_A0 + GPT3_grid.at(ix).dT_A1 * cosfy + GPT3_grid.at(ix).dT_B1 * sinfy + GPT3_grid.at(ix).dT_A2 * coshy + GPT3_grid.at(ix).dT_B2 * sinhy) |
| 414 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | / 1e3; |
| 415 | |||
| 416 | // temperature reduction to station height | ||
| 417 |
2/4✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
|
2612 | Tl.at(l) = T0 + dTl.at(l) * redh - 273.150; |
| 418 | |||
| 419 | // virtual temperature | ||
| 420 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | double Tv = T0 * (1 + 0.6077 * Ql.at(l)); |
| 421 | 2612 | double c = InsConst::G_NORM * InsConst::dMtr / (InsConst::Rg * Tv); | |
| 422 | |||
| 423 | // pressure in hPa | ||
| 424 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | pl.at(l) = (p0 * exp(-c * redh)) / 100.0; |
| 425 | |||
| 426 | // hydrostatic coefficient ah | ||
| 427 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | ahl.at(l) = (GPT3_grid.at(ix).h_A0 + GPT3_grid.at(ix).h_A1 * cosfy + GPT3_grid.at(ix).h_B1 * sinfy + GPT3_grid.at(ix).h_A2 * coshy + GPT3_grid.at(ix).h_B2 * sinhy) |
| 428 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | / 1e3; |
| 429 | |||
| 430 | // wet coefficient aw | ||
| 431 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | awl.at(l) = (GPT3_grid.at(ix).w_A0 + GPT3_grid.at(ix).w_A1 * cosfy + GPT3_grid.at(ix).w_B1 * sinfy + GPT3_grid.at(ix).w_A2 * coshy + GPT3_grid.at(ix).w_B2 * sinhy) |
| 432 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | / 1e3; |
| 433 | |||
| 434 | // water vapor decrease factor la - added by GP | ||
| 435 |
6/12✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2612 times.
✗ Branch 17 not taken.
|
2612 | lal.at(l) = GPT3_grid.at(ix).la_A0 + GPT3_grid.at(ix).la_A1 * cosfy + GPT3_grid.at(ix).la_B1 * sinfy + GPT3_grid.at(ix).la_A2 * coshy + GPT3_grid.at(ix).la_B2 * sinhy; |
| 436 | |||
| 437 | // mean temperature of the water vapor Tm - added by GP | ||
| 438 |
6/12✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2612 times.
✗ Branch 17 not taken.
|
2612 | Tml.at(l) = GPT3_grid.at(ix).Tm_A0 + GPT3_grid.at(ix).Tm_A1 * cosfy + GPT3_grid.at(ix).Tm_B1 * sinfy + GPT3_grid.at(ix).Tm_A2 * coshy + GPT3_grid.at(ix).Tm_B2 * sinhy; |
| 439 | |||
| 440 | // water vapor pressure in hPa - changed by GP | ||
| 441 |
2/4✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
|
2612 | double e0 = Ql.at(l) * p0 / (0.622 + 0.378 * Ql.at(l)) / 100.0; // on the grid |
| 442 | |||
| 443 |
3/6✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
|
2612 | el.at(l) = e0 * pow((100.0 * pl.at(l) / p0), (lal.at(l) + 1.0)); // on the station height (14) Askne and Nordius, 1987 |
| 444 | |||
| 445 | // north and east gradients (total, hydrostatic and wet) | ||
| 446 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | Gn_hl.at(l) = (GPT3_grid.at(ix).Gnh_A0 + GPT3_grid.at(ix).Gnh_A1 * cosfy + GPT3_grid.at(ix).Gnh_B1 * sinfy + GPT3_grid.at(ix).Gnh_A2 * coshy + GPT3_grid.at(ix).Gnh_B2 * sinhy) |
| 447 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | / 1e5; |
| 448 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | Ge_hl.at(l) = (GPT3_grid.at(ix).Geh_A0 + GPT3_grid.at(ix).Geh_A1 * cosfy + GPT3_grid.at(ix).Geh_B1 * sinfy + GPT3_grid.at(ix).Geh_A2 * coshy + GPT3_grid.at(ix).Geh_B2 * sinhy) |
| 449 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | / 1e5; |
| 450 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
2612 | Gn_wl.at(l) = (GPT3_grid.at(ix).Gnw_A0 + GPT3_grid.at(ix).Gnw_A1 * cosfy + GPT3_grid.at(ix).Gnw_B1 * sinfy + GPT3_grid.at(ix).Gnw_A2 * coshy + GPT3_grid.at(ix).Gnw_B2 * sinhy) |
| 451 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | / 1e5; |
| 452 |
5/10✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2612 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2612 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2612 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2612 times.
✗ Branch 14 not taken.
|
5224 | Ge_wl.at(l) = (GPT3_grid.at(ix).Gew_A0 + GPT3_grid.at(ix).Gew_A1 * cosfy + GPT3_grid.at(ix).Gew_B1 * sinfy + GPT3_grid.at(ix).Gew_A2 * coshy + GPT3_grid.at(ix).Gew_B2 * sinhy) |
| 453 |
1/2✓ Branch 1 taken 2612 times.
✗ Branch 2 not taken.
|
2612 | / 1e5; |
| 454 | } | ||
| 455 | |||
| 456 | 653 | double dnpod1 = fabs(diffpod); // distance nearer point | |
| 457 | 653 | double dnpod2 = 1.0 - dnpod1; // distance to distant point | |
| 458 | 653 | double dnlon1 = fabs(difflon); | |
| 459 | 653 | double dnlon2 = 1.0 - dnlon1; | |
| 460 | |||
| 461 | // pressure | ||
| 462 | 653 | double R1 = dnpod2 * pl[0] + dnpod1 * pl[1]; | |
| 463 | 653 | double R2 = dnpod2 * pl[2] + dnpod1 * pl[3]; | |
| 464 | 653 | gpt3outputs.p = dnlon2 * R1 + dnlon1 * R2; | |
| 465 | |||
| 466 | // temperature | ||
| 467 | 653 | R1 = dnpod2 * Tl[0] + dnpod1 * Tl[1]; | |
| 468 | 653 | R2 = dnpod2 * Tl[2] + dnpod1 * Tl[3]; | |
| 469 | 653 | gpt3outputs.T = dnlon2 * R1 + dnlon1 * R2; | |
| 470 | |||
| 471 | // temperature in degree per km | ||
| 472 | 653 | R1 = dnpod2 * dTl[0] + dnpod1 * dTl[1]; | |
| 473 | 653 | R2 = dnpod2 * dTl[2] + dnpod1 * dTl[3]; | |
| 474 | 653 | gpt3outputs.dT = (dnlon2 * R1 + dnlon1 * R2) * 1e3; | |
| 475 | |||
| 476 | // water vapor pressure in hPa - changed by GP | ||
| 477 | 653 | R1 = dnpod2 * el[0] + dnpod1 * el[1]; | |
| 478 | 653 | R2 = dnpod2 * el[2] + dnpod1 * el[3]; | |
| 479 | 653 | gpt3outputs.e = dnlon2 * R1 + dnlon1 * R2; | |
| 480 | |||
| 481 | // hydrostatic | ||
| 482 | 653 | R1 = dnpod2 * ahl[0] + dnpod1 * ahl[1]; | |
| 483 | 653 | R2 = dnpod2 * ahl[2] + dnpod1 * ahl[3]; | |
| 484 | 653 | gpt3outputs.ah = dnlon2 * R1 + dnlon1 * R2; | |
| 485 | |||
| 486 | // wet | ||
| 487 | 653 | R1 = dnpod2 * awl[0] + dnpod1 * awl[1]; | |
| 488 | 653 | R2 = dnpod2 * awl[2] + dnpod1 * awl[3]; | |
| 489 | 653 | gpt3outputs.aw = dnlon2 * R1 + dnlon1 * R2; | |
| 490 | |||
| 491 | // undulation | ||
| 492 | 653 | R1 = dnpod2 * undul[0] + dnpod1 * undul[1]; | |
| 493 | 653 | R2 = dnpod2 * undul[2] + dnpod1 * undul[3]; | |
| 494 | 653 | gpt3outputs.undu = dnlon2 * R1 + dnlon1 * R2; | |
| 495 | |||
| 496 | // water vapor decrease factor la - added by GP | ||
| 497 | 653 | R1 = dnpod2 * lal[0] + dnpod1 * lal[1]; | |
| 498 | 653 | R2 = dnpod2 * lal[2] + dnpod1 * lal[3]; | |
| 499 | 653 | gpt3outputs.la = dnlon2 * R1 + dnlon1 * R2; | |
| 500 | |||
| 501 | // mean temperature of the water vapor - added by GP | ||
| 502 | 653 | R1 = dnpod2 * Tml[0] + dnpod1 * Tml[1]; | |
| 503 | 653 | R2 = dnpod2 * Tml[2] + dnpod1 * Tml[3]; | |
| 504 | 653 | gpt3outputs.Tm = dnlon2 * R1 + dnlon1 * R2; | |
| 505 | |||
| 506 | // gradients | ||
| 507 | 653 | R1 = dnpod2 * Gn_hl[0] + dnpod1 * Gn_hl[1]; | |
| 508 | 653 | R2 = dnpod2 * Gn_hl[2] + dnpod1 * Gn_hl[3]; | |
| 509 | 653 | gpt3outputs.Gn_h = dnlon2 * R1 + dnlon1 * R2; | |
| 510 | |||
| 511 | 653 | R1 = dnpod2 * Ge_hl[0] + dnpod1 * Ge_hl[1]; | |
| 512 | 653 | R2 = dnpod2 * Ge_hl[2] + dnpod1 * Ge_hl[3]; | |
| 513 | 653 | gpt3outputs.Ge_h = dnlon2 * R1 + dnlon1 * R2; | |
| 514 | |||
| 515 | 653 | R1 = dnpod2 * Gn_wl[0] + dnpod1 * Gn_wl[1]; | |
| 516 | 653 | R2 = dnpod2 * Gn_wl[2] + dnpod1 * Gn_wl[3]; | |
| 517 | 653 | gpt3outputs.Gn_w = dnlon2 * R1 + dnlon1 * R2; | |
| 518 | |||
| 519 | 653 | R1 = dnpod2 * Ge_wl[0] + dnpod1 * Ge_wl[1]; | |
| 520 | 653 | R2 = dnpod2 * Ge_wl[2] + dnpod1 * Ge_wl[3]; | |
| 521 | 653 | gpt3outputs.Ge_w = dnlon2 * R1 + dnlon1 * R2; | |
| 522 | } | ||
| 523 | |||
| 524 | 1450 | return gpt3outputs; | |
| 525 | } | ||
| 526 | |||
| 527 | 725 | double mjd2doy(const double& mjd) | |
| 528 | { | ||
| 529 | // convert mjd to doy | ||
| 530 | 725 | double hour = floor((mjd - floor(mjd)) * 24); // get hours | |
| 531 | 725 | double minu = floor((((mjd - floor(mjd)) * 24) - hour) * 60); // get minutes | |
| 532 | 725 | double sec = (((((mjd - floor(mjd)) * 24) - hour) * 60) - minu) * 60; // get seconds | |
| 533 | |||
| 534 | // change secs, min hour whose sec==60 and days, whose hour==24 | ||
| 535 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 725 times.
|
725 | if (sec == 60.0) { minu = minu + 1; } |
| 536 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 725 times.
|
725 | if (minu == 60.0) { hour = hour + 1; } |
| 537 | |||
| 538 | // calc jd (yet wrong for hour==24) | ||
| 539 | 725 | double jd = mjd + 2400000.5; | |
| 540 | |||
| 541 | // if hour==24, correct jd and set hour==0 | ||
| 542 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 725 times.
|
725 | if (hour == 24.0) { jd = jd + 1; } |
| 543 | |||
| 544 | // integer Julian date | ||
| 545 | 725 | double jd_int = floor(jd + 0.5); | |
| 546 | |||
| 547 | 725 | double aa = jd_int + 32044; | |
| 548 | 725 | double bb = floor((4 * aa + 3) / 146097); | |
| 549 | 725 | double cc = aa - floor((bb * 146097) / 4); | |
| 550 | 725 | double dd = floor((4 * cc + 3) / 1461); | |
| 551 | 725 | double ee = cc - floor((1461 * dd) / 4); | |
| 552 | 725 | double mm = floor((5 * ee + 2) / 153); | |
| 553 | |||
| 554 | 725 | auto day = static_cast<int>(ee - floor((153 * mm + 2) / 5) + 1); | |
| 555 | 725 | auto month = static_cast<int>(mm + 3 - 12 * floor(mm / 10)); | |
| 556 | 725 | auto year = static_cast<int>(bb * 100 + dd - 4800 + floor(mm / 10)); | |
| 557 | |||
| 558 | // first check if the specified year is leap year or not (logical output) | ||
| 559 | 725 | double leapYear = 0.0; | |
| 560 |
2/6✓ Branch 0 taken 725 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 725 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
725 | if ((year % 4 == 0 && year % 100 != 0) || (year % 400 == 0)) |
| 561 | { | ||
| 562 | 725 | leapYear = 1.0; | |
| 563 | } | ||
| 564 | else | ||
| 565 | { | ||
| 566 | ✗ | leapYear = 0.0; | |
| 567 | } | ||
| 568 |
2/4✓ Branch 1 taken 725 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 725 times.
✗ Branch 7 not taken.
|
2900 | auto days = [inits = { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 }] { return std::vector<int>{ std::begin(inits), std::end(inits) }; }(); |
| 569 | 725 | int sum = 0; | |
| 570 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 725 times.
|
725 | for (size_t i = 0; i < static_cast<size_t>(month) - 1; i++) |
| 571 | { | ||
| 572 | ✗ | sum = sum + days[i]; | |
| 573 | } | ||
| 574 | 725 | double doy = sum + day; | |
| 575 | |||
| 576 |
2/4✓ Branch 0 taken 725 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 725 times.
|
725 | if (leapYear == 1 && month > 2) |
| 577 | { | ||
| 578 | ✗ | doy = doy + 1; | |
| 579 | } | ||
| 580 | |||
| 581 | // Add decimal places | ||
| 582 | 725 | doy = doy + mjd - floor(mjd); | |
| 583 | |||
| 584 | 725 | return doy; | |
| 585 | 725 | } | |
| 586 | |||
| 587 | ✗ | double asknewet(const double& e, const double& Tm, const double& la) | |
| 588 | { | ||
| 589 | // coefficients | ||
| 590 | ✗ | double k1 = 77.604; // K/hPa | |
| 591 | ✗ | double k2 = 64.79; // K/hPa | |
| 592 | ✗ | double k2p = k2 - k1 * 18.0152 / 28.9644; // K/hPa | |
| 593 | ✗ | double k3 = 377600.0; // KK/hPa | |
| 594 | |||
| 595 | // specific gas constant for dry consituents | ||
| 596 | ✗ | constexpr double Rd = InsConst::Rg / InsConst::dMtr; | |
| 597 | |||
| 598 | ✗ | return 1.0e-6 * (k2p + k3 / Tm) * Rd / (la + 1.0) / InsConst::G_NORM * e; | |
| 599 | } | ||
| 600 | |||
| 601 | } // namespace NAV | ||
| 602 |