INSTINCT Code Coverage Report

Directory: src/
File: Navigation/Atmosphere/Troposphere/MappingFunctions/GMF.cpp
Date: 2025-02-07 16:54:41
Exec Total Coverage
Lines: 87 87 100.0%
Functions: 3 3 100.0%
Branches: 73 128 57.0%
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 /// @file GMF.cpp
10 /// @brief Global Mapping Function (GMF)
11 /// @author T. Topp (topp@ins.uni-stuttgart.de)
12 /// @date 2024-04-21
13 /// @note See \cite Böhm2006a Böhm2006: Global Mapping Function (GMF): A new empirical mapping function based on numerical weather model data
14 /// @note See https://vmf.geo.tuwien.ac.at/codes/ for code sources in matlab.
15
16 #include "GMF.hpp"
17 #include "internal/GMFCoeffs.hpp"
18
19 namespace NAV::internal::GMF
20 {
21
22 namespace
23 {
24
25 // degree n and order m
26 constexpr int nmax = 9;
27
28 14220 std::array<Eigen::Matrix<double, nmax + 1, nmax + 1>, 2> calcLegendrePolynomials(const Eigen::Vector3d& lla_pos)
29 {
30 // unit vector
31
2/4
✓ Branch 1 taken 14220 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14220 times.
✗ Branch 5 not taken.
 14220 double x = std::cos(lla_pos(0)) * std::cos(lla_pos(1));
32
2/4
✓ Branch 1 taken 14220 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14220 times.
✗ Branch 5 not taken.
 14220 double y = std::cos(lla_pos(0)) * std::sin(lla_pos(1));
33
1/2
✓ Branch 1 taken 14220 times.
✗ Branch 2 not taken.
 14220 double z = std::sin(lla_pos(0));
34
35 // Legendre polynomials
36
1/2
✓ Branch 1 taken 14220 times.
✗ Branch 2 not taken.
 14220 Eigen::Matrix<double, nmax + 1, nmax + 1> V;
37
1/2
✓ Branch 1 taken 14220 times.
✗ Branch 2 not taken.
 14220 Eigen::Matrix<double, nmax + 1, nmax + 1> W;
38
39
1/2
✓ Branch 1 taken 14220 times.
✗ Branch 2 not taken.
 14220 V(0, 0) = 1;
40
1/2
✓ Branch 1 taken 14220 times.
✗ Branch 2 not taken.
 14220 W(0, 0) = 0;
41
2/4
✓ Branch 1 taken 14220 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14220 times.
✗ Branch 5 not taken.
 14220 V(1, 0) = z * V(0, 0);
42
1/2
✓ Branch 1 taken 14220 times.
✗ Branch 2 not taken.
 14220 W(1, 0) = 0;
43
44
2/2
✓ Branch 0 taken 113760 times.
✓ Branch 1 taken 14220 times.
 127980 for (int n = 2; n <= nmax; n++)
45 {
46 113760 auto dn = static_cast<double>(n);
47
3/6
✓ Branch 1 taken 113760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 113760 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 113760 times.
✗ Branch 8 not taken.
 113760 V(n, 0) = ((2 * dn - 1) * z * V(n - 1, 0) - (dn - 1) * V(n - 2, 0)) / dn;
48
1/2
✓ Branch 1 taken 113760 times.
✗ Branch 2 not taken.
 113760 W(n, 0) = 0;
49 }
50
2/2
✓ Branch 0 taken 127980 times.
✓ Branch 1 taken 14220 times.
 142200 for (int m = 1; m <= nmax; m++)
51 {
52 127980 auto dm = static_cast<double>(m);
53
3/6
✓ Branch 1 taken 127980 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 127980 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 127980 times.
✗ Branch 8 not taken.
 127980 V(m, m) = (2 * dm - 1) * (x * V(m - 1, m - 1) - y * W(m - 1, m - 1));
54
3/6
✓ Branch 1 taken 127980 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 127980 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 127980 times.
✗ Branch 8 not taken.
 127980 W(m, m) = (2 * dm - 1) * (x * W(m - 1, m - 1) + y * V(m - 1, m - 1));
55
2/2
✓ Branch 0 taken 113760 times.
✓ Branch 1 taken 14220 times.
 127980 if (m < nmax)
56 {
57
2/4
✓ Branch 1 taken 113760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 113760 times.
✗ Branch 5 not taken.
 113760 V(m + 1, m) = (2 * dm + 1) * z * V(m, m);
58
2/4
✓ Branch 1 taken 113760 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 113760 times.
✗ Branch 5 not taken.
 113760 W(m + 1, m) = (2 * dm + 1) * z * W(m, m);
59 }
60
2/2
✓ Branch 0 taken 398160 times.
✓ Branch 1 taken 127980 times.
 526140 for (int n = m + 2; n <= nmax; n++)
61 {
62 398160 auto dn = static_cast<double>(n);
63
3/6
✓ Branch 1 taken 398160 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 398160 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 398160 times.
✗ Branch 8 not taken.
 398160 V(n, m) = ((2 * dn - 1) * z * V(n - 1, m) - (dn + dm - 1) * V(n - 2, m)) / (dn - dm);
64
3/6
✓ Branch 1 taken 398160 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 398160 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 398160 times.
✗ Branch 8 not taken.
 398160 W(n, m) = ((2 * dn - 1) * z * W(n - 1, m) - (dn + dm - 1) * W(n - 2, m)) / (dn - dm);
65 }
66 }
67
68
2/4
✓ Branch 1 taken 14220 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14220 times.
✗ Branch 5 not taken.
 28440 return { V, W };
69 }
70
71 } // namespace
72
73 } // namespace NAV::internal::GMF
74
75 7110 double NAV::calcTropoMapFunc_GMFH(double mjd, const Eigen::Vector3d& lla_pos, double elevation)
76 {
77 using namespace internal::GMF; // NOLINT(google-build-using-namespace)
78
79 // reference day is 28 January
80 // this is taken from Niell (1996) to be consistent
81 7110 double doy = mjd - 44239.0 + 1 - 28;
82
83
1/2
✓ Branch 1 taken 7110 times.
✗ Branch 2 not taken.
 7110 auto [V, W] = calcLegendrePolynomials(lla_pos);
84
85 7110 double bh = 0.0029;
86 7110 double c0h = 0.062;
87 7110 double phh = 0.0;
88 7110 double c11h = 0.0;
89 7110 double c10h = 0.0;
90
3/4
✓ Branch 1 taken 7110 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 324 times.
✓ Branch 4 taken 6786 times.
 7110 if (lla_pos(0) < 0) // southern hemisphere
91 {
92 324 phh = M_PI;
93 324 c11h = 0.007;
94 324 c10h = 0.002;
95 }
96 else // northern hemisphere
97 {
98 6786 phh = 0;
99 6786 c11h = 0.005;
100 6786 c10h = 0.001;
101 }
102
1/2
✓ Branch 1 taken 7110 times.
✗ Branch 2 not taken.
 7110 double ch = c0h + ((std::cos(doy / 365.25 * 2 * M_PI + phh) + 1) * c11h / 2 + c10h) * (1 - std::cos(lla_pos(0)));
103
104 7110 double ahm = 0;
105 7110 double aha = 0;
106 7110 size_t i = 0;
107
2/2
✓ Branch 0 taken 71100 times.
✓ Branch 1 taken 7110 times.
 78210 for (int n = 0; n <= nmax; n++)
108 {
109
2/2
✓ Branch 0 taken 391050 times.
✓ Branch 1 taken 71100 times.
 462150 for (int m = 0; m <= n; m++)
110 {
111
4/8
✓ Branch 1 taken 391050 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 391050 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 391050 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 391050 times.
✗ Branch 11 not taken.
 391050 ahm = ahm + (ah_mean.at(i) * V(n, m) + bh_mean.at(i) * W(n, m));
112
4/8
✓ Branch 1 taken 391050 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 391050 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 391050 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 391050 times.
✗ Branch 11 not taken.
 391050 aha = aha + (ah_amp.at(i) * V(n, m) + bh_amp.at(i) * W(n, m));
113 391050 i = i + 1;
114 }
115 }
116 7110 double ah = (ahm + aha * std::cos(doy / 365.25 * 2 * M_PI)) * 1e-5;
117
118 7110 double sine = std::sin(elevation);
119 7110 double beta = bh / (sine + ch);
120 7110 double gamma = ah / (sine + beta);
121 7110 double topcon = (1 + ah / (1 + bh / (1 + ch)));
122 7110 double gmfh = topcon / (sine + gamma);
123
124 // height correction for hydrostatic mapping function from Niell (1996) in order to reduce the coefficients to sea level
125 7110 double a_ht = 2.53e-5;
126 7110 double b_ht = 5.49e-3;
127 7110 double c_ht = 1.14e-3;
128
1/2
✓ Branch 1 taken 7110 times.
✗ Branch 2 not taken.
 7110 double hs_km = lla_pos(2) / 1000;
129
130 7110 beta = b_ht / (sine + c_ht);
131 7110 gamma = a_ht / (sine + beta);
132 7110 topcon = (1 + a_ht / (1 + b_ht / (1 + c_ht)));
133 7110 double ht_corr_coef = 1 / sine - topcon / (sine + gamma);
134 7110 double ht_corr = ht_corr_coef * hs_km;
135 7110 gmfh += ht_corr;
136
137 7110 return gmfh;
138 }
139
140 7110 double NAV::calcTropoMapFunc_GMFW(double mjd, const Eigen::Vector3d& lla_pos, double elevation)
141 {
142 using namespace internal::GMF; // NOLINT(google-build-using-namespace)
143
144 // reference day is 28 January
145 // this is taken from Niell (1996) to be consistent
146 7110 double doy = mjd - 44239.0 + 1 - 28;
147
148
1/2
✓ Branch 1 taken 7110 times.
✗ Branch 2 not taken.
 7110 auto [V, W] = calcLegendrePolynomials(lla_pos);
149
150 7110 double bw = 0.00146;
151 7110 double cw = 0.04391;
152
153 7110 double awm = 0.0;
154 7110 double awa = 0.0;
155 7110 size_t i = 0;
156
2/2
✓ Branch 0 taken 71100 times.
✓ Branch 1 taken 7110 times.
 78210 for (int n = 0; n <= nmax; n++)
157 {
158
2/2
✓ Branch 0 taken 391050 times.
✓ Branch 1 taken 71100 times.
 462150 for (int m = 0; m <= n; m++)
159 {
160
4/8
✓ Branch 1 taken 391050 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 391050 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 391050 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 391050 times.
✗ Branch 11 not taken.
 391050 awm = awm + (aw_mean.at(i) * V(n, m) + bw_mean.at(i) * W(n, m));
161
4/8
✓ Branch 1 taken 391050 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 391050 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 391050 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 391050 times.
✗ Branch 11 not taken.
 391050 awa = awa + (aw_amp.at(i) * V(n, m) + bw_amp.at(i) * W(n, m));
162 391050 i = i + 1;
163 }
164 }
165 7110 double aw = (awm + awa * std::cos(doy / 365.25 * 2 * M_PI)) * 1e-5;
166
167 7110 double sine = std::sin(elevation);
168 7110 double beta = bw / (sine + cw);
169 7110 double gamma = aw / (sine + beta);
170 7110 double topcon = (1 + aw / (1 + bw / (1 + cw)));
171 7110 double gmfw = topcon / (sine + gamma);
172
173 7110 return gmfw;
174 }
175