INSTINCT Code Coverage Report

Directory: src/
File: Navigation/Atmosphere/Troposphere/MappingFunctions/GMF.cpp
Date: 2025-11-25 23:34:18
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 125908 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 125908 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 125908 times.
✗ Branch 5 not taken.
 125908 double x = std::cos(lla_pos(0)) * std::cos(lla_pos(1));
32
2/4
✓ Branch 1 taken 125908 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 125908 times.
✗ Branch 5 not taken.
 125908 double y = std::cos(lla_pos(0)) * std::sin(lla_pos(1));
33
1/2
✓ Branch 1 taken 125908 times.
✗ Branch 2 not taken.
 125908 double z = std::sin(lla_pos(0));
34
35 // Legendre polynomials
36
1/2
✓ Branch 1 taken 125908 times.
✗ Branch 2 not taken.
 125908 Eigen::Matrix<double, nmax + 1, nmax + 1> V;
37
1/2
✓ Branch 1 taken 125908 times.
✗ Branch 2 not taken.
 125908 Eigen::Matrix<double, nmax + 1, nmax + 1> W;
38
39
1/2
✓ Branch 1 taken 125908 times.
✗ Branch 2 not taken.
 125908 V(0, 0) = 1;
40
1/2
✓ Branch 1 taken 125908 times.
✗ Branch 2 not taken.
 125908 W(0, 0) = 0;
41
2/4
✓ Branch 1 taken 125908 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 125908 times.
✗ Branch 5 not taken.
 125908 V(1, 0) = z * V(0, 0);
42
1/2
✓ Branch 1 taken 125908 times.
✗ Branch 2 not taken.
 125908 W(1, 0) = 0;
43
44
2/2
✓ Branch 0 taken 1007264 times.
✓ Branch 1 taken 125908 times.
 1133172 for (int n = 2; n <= nmax; n++)
45 {
46 1007264 auto dn = static_cast<double>(n);
47
3/6
✓ Branch 1 taken 1007264 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1007264 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1007264 times.
✗ Branch 8 not taken.
 1007264 V(n, 0) = ((2 * dn - 1) * z * V(n - 1, 0) - (dn - 1) * V(n - 2, 0)) / dn;
48
1/2
✓ Branch 1 taken 1007264 times.
✗ Branch 2 not taken.
 1007264 W(n, 0) = 0;
49 }
50
2/2
✓ Branch 0 taken 1133172 times.
✓ Branch 1 taken 125908 times.
 1259080 for (int m = 1; m <= nmax; m++)
51 {
52 1133172 auto dm = static_cast<double>(m);
53
3/6
✓ Branch 1 taken 1133172 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1133172 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1133172 times.
✗ Branch 8 not taken.
 1133172 V(m, m) = (2 * dm - 1) * (x * V(m - 1, m - 1) - y * W(m - 1, m - 1));
54
3/6
✓ Branch 1 taken 1133172 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1133172 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1133172 times.
✗ Branch 8 not taken.
 1133172 W(m, m) = (2 * dm - 1) * (x * W(m - 1, m - 1) + y * V(m - 1, m - 1));
55
2/2
✓ Branch 0 taken 1007264 times.
✓ Branch 1 taken 125908 times.
 1133172 if (m < nmax)
56 {
57
2/4
✓ Branch 1 taken 1007264 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1007264 times.
✗ Branch 5 not taken.
 1007264 V(m + 1, m) = (2 * dm + 1) * z * V(m, m);
58
2/4
✓ Branch 1 taken 1007264 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1007264 times.
✗ Branch 5 not taken.
 1007264 W(m + 1, m) = (2 * dm + 1) * z * W(m, m);
59 }
60
2/2
✓ Branch 0 taken 3525424 times.
✓ Branch 1 taken 1133172 times.
 4658596 for (int n = m + 2; n <= nmax; n++)
61 {
62 3525424 auto dn = static_cast<double>(n);
63
3/6
✓ Branch 1 taken 3525424 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3525424 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3525424 times.
✗ Branch 8 not taken.
 3525424 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 3525424 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3525424 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3525424 times.
✗ Branch 8 not taken.
 3525424 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 125908 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 125908 times.
✗ Branch 5 not taken.
 251816 return { V, W };
69 }
70
71 } // namespace
72
73 } // namespace NAV::internal::GMF
74
75 62954 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 62954 double doy = mjd - 44239.0 + 1 - 28;
82
83
1/2
✓ Branch 1 taken 62954 times.
✗ Branch 2 not taken.
 62954 auto [V, W] = calcLegendrePolynomials(lla_pos);
84
85 62954 double bh = 0.0029;
86 62954 double c0h = 0.062;
87 62954 double phh = 0.0;
88 62954 double c11h = 0.0;
89 62954 double c10h = 0.0;
90
3/4
✓ Branch 1 taken 62954 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 324 times.
✓ Branch 4 taken 62630 times.
 62954 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 62630 phh = 0;
99 62630 c11h = 0.005;
100 62630 c10h = 0.001;
101 }
102
1/2
✓ Branch 1 taken 62954 times.
✗ Branch 2 not taken.
 62954 double ch = c0h + ((std::cos(doy / 365.25 * 2 * M_PI + phh) + 1) * c11h / 2 + c10h) * (1 - std::cos(lla_pos(0)));
103
104 62954 double ahm = 0;
105 62954 double aha = 0;
106 62954 size_t i = 0;
107
2/2
✓ Branch 0 taken 629540 times.
✓ Branch 1 taken 62954 times.
 692494 for (int n = 0; n <= nmax; n++)
108 {
109
2/2
✓ Branch 0 taken 3462470 times.
✓ Branch 1 taken 629540 times.
 4092010 for (int m = 0; m <= n; m++)
110 {
111
4/8
✓ Branch 1 taken 3462470 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3462470 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3462470 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3462470 times.
✗ Branch 11 not taken.
 3462470 ahm = ahm + (ah_mean.at(i) * V(n, m) + bh_mean.at(i) * W(n, m));
112
4/8
✓ Branch 1 taken 3462470 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3462470 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3462470 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3462470 times.
✗ Branch 11 not taken.
 3462470 aha = aha + (ah_amp.at(i) * V(n, m) + bh_amp.at(i) * W(n, m));
113 3462470 i = i + 1;
114 }
115 }
116 62954 double ah = (ahm + aha * std::cos(doy / 365.25 * 2 * M_PI)) * 1e-5;
117
118 62954 double sine = std::sin(elevation);
119 62954 double beta = bh / (sine + ch);
120 62954 double gamma = ah / (sine + beta);
121 62954 double topcon = (1 + ah / (1 + bh / (1 + ch)));
122 62954 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 62954 double a_ht = 2.53e-5;
126 62954 double b_ht = 5.49e-3;
127 62954 double c_ht = 1.14e-3;
128
1/2
✓ Branch 1 taken 62954 times.
✗ Branch 2 not taken.
 62954 double hs_km = lla_pos(2) / 1000;
129
130 62954 beta = b_ht / (sine + c_ht);
131 62954 gamma = a_ht / (sine + beta);
132 62954 topcon = (1 + a_ht / (1 + b_ht / (1 + c_ht)));
133 62954 double ht_corr_coef = 1 / sine - topcon / (sine + gamma);
134 62954 double ht_corr = ht_corr_coef * hs_km;
135 62954 gmfh += ht_corr;
136
137 62954 return gmfh;
138 }
139
140 62954 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 62954 double doy = mjd - 44239.0 + 1 - 28;
147
148
1/2
✓ Branch 1 taken 62954 times.
✗ Branch 2 not taken.
 62954 auto [V, W] = calcLegendrePolynomials(lla_pos);
149
150 62954 double bw = 0.00146;
151 62954 double cw = 0.04391;
152
153 62954 double awm = 0.0;
154 62954 double awa = 0.0;
155 62954 size_t i = 0;
156
2/2
✓ Branch 0 taken 629540 times.
✓ Branch 1 taken 62954 times.
 692494 for (int n = 0; n <= nmax; n++)
157 {
158
2/2
✓ Branch 0 taken 3462470 times.
✓ Branch 1 taken 629540 times.
 4092010 for (int m = 0; m <= n; m++)
159 {
160
4/8
✓ Branch 1 taken 3462470 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3462470 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3462470 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3462470 times.
✗ Branch 11 not taken.
 3462470 awm = awm + (aw_mean.at(i) * V(n, m) + bw_mean.at(i) * W(n, m));
161
4/8
✓ Branch 1 taken 3462470 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3462470 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3462470 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3462470 times.
✗ Branch 11 not taken.
 3462470 awa = awa + (aw_amp.at(i) * V(n, m) + bw_amp.at(i) * W(n, m));
162 3462470 i = i + 1;
163 }
164 }
165 62954 double aw = (awm + awa * std::cos(doy / 365.25 * 2 * M_PI)) * 1e-5;
166
167 62954 double sine = std::sin(elevation);
168 62954 double beta = bw / (sine + cw);
169 62954 double gamma = aw / (sine + beta);
170 62954 double topcon = (1 + aw / (1 + bw / (1 + cw)));
171 62954 double gmfw = topcon / (sine + gamma);
172
173 62954 return gmfw;
174 }
175