0.3.0
Loading...
Searching...
No Matches
GPT.cpp
Go to the documentation of this file.
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
20namespace NAV
21{
22GPT2output GPT2_param(const double& mjd, const Eigen::Vector3d& lla_pos)
23{
24 using internal::GPT2_grid;
25
26 GPT2output gpt2outputs;
27
28 // change the reference epoch to January 1 2000
29 double mjd1 = mjd - 51544.5;
30 // GPT2 with time variation
31 double cosfy = cos(mjd1 / 365.25 * 2.0 * M_PI);
32 double coshy = cos(mjd1 / 365.25 * 4.0 * M_PI);
33 double sinfy = sin(mjd1 / 365.25 * 2.0 * M_PI);
34 double sinhy = sin(mjd1 / 365.25 * 4.0 * M_PI);
35
36 double dlat = lla_pos(0);
37 double dlon = lla_pos(1);
38 double hell = lla_pos(2);
39
40 // only positive longitude in degrees
41 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 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 auto ipod = static_cast<int>(floor(ppod + 1.0));
48 auto ilon = static_cast<int>(floor(plon + 1.0));
49
50 // normalized (to one) differences, can be positive or negative
51 double diffpod = ppod - (ipod - 0.5);
52 double difflon = plon - (ilon - 0.5);
53
54 // added by HCY
55 if (ipod == 181) { ipod = 180; }
56 // added by GP
57 if (ilon == 361) { ilon = 1; }
58 if (ilon == 0) { ilon = 360; }
59
60 constexpr size_t N = 4;
61 std::array<size_t, N> indx{};
62 // get the number of the corresponding line
63 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 if (ppod <= 0.5 || ppod >= 179.5) // case of nearest neighborhood
68 {
69 auto ix = indx[0];
70 // transforming ellipsoidal height to orthometric height
71 gpt2outputs.undu = GPT2_grid.at(ix).undulation; // undu: geoid undulation in m
72 double hgt = hell - gpt2outputs.undu;
73
74 // pressure, temperature at the height of the grid
75 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 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 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 / 1e3;
81
82 // lapse rate of the temperature
83 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 / 1e3;
85
86 // station height - grid height
87 double redh = hgt - GPT2_grid.at(ix).Hs; // Hs: orthometric grid height in m
88
89 // temperature at station height in Celsius
90 gpt2outputs.T = T0 + gpt2outputs.dT * redh - 273.150;
91
92 // temperature lapse rate in degrees / km
93 gpt2outputs.dT = gpt2outputs.dT * 1e3;
94
95 // virtual temperature in Kelvin
96 double Tv = T0 * (1 + 0.6077 * Q);
97
98 double c = InsConst::G_NORM * InsConst::dMtr / (InsConst::Rg * Tv);
99
100 // pressure in hPa
101 gpt2outputs.p = (p0 * exp(-c * redh)) / 100.0;
102
103 // hydrostatic coefficient ah
104 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 / 1e3;
106
107 // wet coefficient aw
108 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 / 1e3;
110
111 // water vapour decrease factor la - added by GP
112 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 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 double e0 = Q * p0 / (0.622 + 0.378 * Q) / 100.0; // on the grid
119
120 gpt2outputs.e = e0 * pow((100.0 * gpt2outputs.p / p0), (gpt2outputs.la + 1)); // on the station height - (14) Askne and Nordius, 1987
121 }
122 else // bilinear interpolation
123 {
124 double ipod1 = ipod + int(math::sign(1.0, diffpod));
125 double ilon1 = ilon + int(math::sign(1.0, difflon));
126
127 // changed for the 1 degree grid (GP)
128 if (ilon1 == 361) { ilon1 = 1; }
129 if (ilon1 == 0) { ilon1 = 360; }
130 // get the number of the line
131 indx[1] = static_cast<size_t>((ipod1 - 1) * 360 + ilon) - 1; // along same longitude
132 indx[2] = static_cast<size_t>((ipod - 1) * 360 + ilon1) - 1; // along same polar distance
133 indx[3] = static_cast<size_t>((ipod1 - 1) * 360 + ilon1) - 1; // diagonal
134
135 std::array<double, N> undul{};
136 std::array<double, N> Ql{};
137 std::array<double, N> dTl{};
138 std::array<double, N> Tl{};
139 std::array<double, N> pl{};
140 std::array<double, N> ahl{};
141 std::array<double, N> awl{};
142 std::array<double, N> lal{};
143 std::array<double, N> Tml{};
144 std::array<double, N> el{};
145
146 for (size_t l = 0; l < N; l++)
147 {
148 // transforming ellipsoidal height to orthometric height:
149 // Hortho = -N + Hell
150 auto ix = indx.at(l);
151
152 undul.at(l) = GPT2_grid.at(ix).undulation; // undu: geoid undulation in m
153 double hgt = hell - undul.at(l);
154
155 // pressure, temperature at the height of the grid
156 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 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 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 / 1e3;
162
163 // reduction = stationheight - gridheight
164 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 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 / 1e3;
169
170 // temperature reduction to station height
171 Tl.at(l) = T0 + dTl.at(l) * redh - 273.150;
172
173 // virtual temperature
174 double Tv = T0 * (1 + 0.6077 * Ql.at(l));
175 double c = InsConst::G_NORM * InsConst::dMtr / (InsConst::Rg * Tv);
176
177 // pressure in hPa
178 pl.at(l) = (p0 * exp(-c * redh)) / 100.0;
179
180 // hydrostatic coefficient ah
181 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 / 1e3;
183
184 // wet coefficient aw
185 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 / 1e3;
187
188 // water vapor decrease factor la - added by GP
189 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 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 double e0 = Ql.at(l) * p0 / (0.622 + 0.378 * Ql.at(l)) / 100.0; // on the grid
196
197 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 double dnpod1 = fabs(diffpod); // distance nearer point
201 double dnpod2 = 1.0 - dnpod1; // distance to distant point
202 double dnlon1 = fabs(difflon);
203 double dnlon2 = 1.0 - dnlon1;
204
205 // pressure
206 double R1 = dnpod2 * pl[0] + dnpod1 * pl[1];
207 double R2 = dnpod2 * pl[2] + dnpod1 * pl[3];
208 gpt2outputs.p = dnlon2 * R1 + dnlon1 * R2;
209
210 // temperature
211 R1 = dnpod2 * Tl[0] + dnpod1 * Tl[1];
212 R2 = dnpod2 * Tl[2] + dnpod1 * Tl[3];
213 gpt2outputs.T = dnlon2 * R1 + dnlon1 * R2;
214
215 // temperature in degree per km
216 R1 = dnpod2 * dTl[0] + dnpod1 * dTl[1];
217 R2 = dnpod2 * dTl[2] + dnpod1 * dTl[3];
218 gpt2outputs.dT = (dnlon2 * R1 + dnlon1 * R2) * 1e3;
219
220 // water vapor pressure in hPa - changed by GP
221 R1 = dnpod2 * el[0] + dnpod1 * el[1];
222 R2 = dnpod2 * el[2] + dnpod1 * el[3];
223 gpt2outputs.e = dnlon2 * R1 + dnlon1 * R2;
224
225 // hydrostatic
226 R1 = dnpod2 * ahl[0] + dnpod1 * ahl[1];
227 R2 = dnpod2 * ahl[2] + dnpod1 * ahl[3];
228 gpt2outputs.ah = dnlon2 * R1 + dnlon1 * R2;
229
230 // wet
231 R1 = dnpod2 * awl[0] + dnpod1 * awl[1];
232 R2 = dnpod2 * awl[2] + dnpod1 * awl[3];
233 gpt2outputs.aw = dnlon2 * R1 + dnlon1 * R2;
234
235 // undulation
236 R1 = dnpod2 * undul[0] + dnpod1 * undul[1];
237 R2 = dnpod2 * undul[2] + dnpod1 * undul[3];
238 gpt2outputs.undu = dnlon2 * R1 + dnlon1 * R2;
239
240 // water vapor decrease factor la - added by GP
241 R1 = dnpod2 * lal[0] + dnpod1 * lal[1];
242 R2 = dnpod2 * lal[2] + dnpod1 * lal[3];
243 gpt2outputs.la = dnlon2 * R1 + dnlon1 * R2;
244
245 // mean temperature of the water vapor - added by GP
246 R1 = dnpod2 * Tml[0] + dnpod1 * Tml[1];
247 R2 = dnpod2 * Tml[2] + dnpod1 * Tml[3];
248 gpt2outputs.Tm = dnlon2 * R1 + dnlon1 * R2;
249 }
250
251 return gpt2outputs;
252}
253
254GPT3output GPT3_param(const double& mjd, const Eigen::Vector3d& lla_pos)
255{
256 using internal::GPT3_grid;
257
258 GPT3output gpt3outputs;
259
260 // convert mjd to doy
261 double doy = mjd2doy(mjd);
262 // factors for amplitudes
263 double cosfy = cos(doy / 365.25 * 2.0 * M_PI);
264 double coshy = cos(doy / 365.25 * 4.0 * M_PI);
265 double sinfy = sin(doy / 365.25 * 2.0 * M_PI);
266 double sinhy = sin(doy / 365.25 * 4.0 * M_PI);
267
268 double dlat = lla_pos(0);
269 double dlon = lla_pos(1);
270 double hell = lla_pos(2);
271
272 // only positive longitude in degrees
273 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 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 auto ipod = static_cast<int>(floor(ppod + 1.0));
280 auto ilon = static_cast<int>(floor(plon + 1.0));
281
282 // normalized (to one) differences, can be positive or negative
283 double diffpod = ppod - (ipod - 0.5);
284 double difflon = plon - (ilon - 0.5);
285
286 // added by HCY
287 if (ipod == 181) { ipod = 180; }
288 // added by GP
289 if (ilon == 361) { ilon = 1; }
290 if (ilon == 0) { ilon = 360; }
291
292 constexpr size_t N = 4;
293 std::array<size_t, N> indx{};
294 // get the number of the corresponding line
295 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 if (ppod <= 0.5 || ppod >= 179.5) // case of nearest neighborhood
300 {
301 auto ix = indx[0];
302 // transforming ellipsoidal height to orthometric height
303 gpt3outputs.undu = GPT3_grid.at(ix).undulation; // undu: geoid undulation in m
304 double hgt = hell - gpt3outputs.undu;
305
306 // pressure, temperature at the height of the grid
307 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 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 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 / 1e3;
313
314 // lapse rate of the temperature
315 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 / 1e3;
317
318 // station height - grid height
319 double redh = hgt - GPT3_grid.at(ix).Hs; // Hs: orthometric grid height in m
320
321 // temperature at station height in Celsius
322 gpt3outputs.T = T0 + gpt3outputs.dT * redh - 273.150;
323
324 // temperature lapse rate in degrees / km
325 gpt3outputs.dT = gpt3outputs.dT * 1e3;
326
327 // virtual temperature in Kelvin
328 double Tv = T0 * (1 + 0.6077 * Q);
329
330 double c = InsConst::G_NORM * InsConst::dMtr / (InsConst::Rg * Tv);
331
332 // pressure in hPa
333 gpt3outputs.p = (p0 * exp(-c * redh)) / 100.0;
334
335 // hydrostatic coefficient ah
336 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 / 1e3;
338
339 // wet coefficient aw
340 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 / 1e3;
342
343 // water vapour decrease factor la - added by GP
344 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 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 double e0 = Q * p0 / (0.622 + 0.378 * Q) / 100.0; // on the grid
351
352 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 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 / 1e5;
357 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 / 1e5;
359 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 / 1e5;
361 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 / 1e5;
363 }
364 else // bilinear interpolation
365 {
366 double ipod1 = ipod + int(math::sign(1.0, diffpod));
367 double ilon1 = ilon + int(math::sign(1.0, difflon));
368
369 // changed for the 1 degree grid (GP)
370 if (ilon1 == 361) { ilon1 = 1; }
371 if (ilon1 == 0) { ilon1 = 360; }
372 // get the number of the line
373 indx[1] = static_cast<size_t>((ipod1 - 1) * 360 + ilon) - 1; // along same longitude
374 indx[2] = static_cast<size_t>((ipod - 1) * 360 + ilon1) - 1; // along same polar distance
375 indx[3] = static_cast<size_t>((ipod1 - 1) * 360 + ilon1) - 1; // diagonal
376
377 std::array<double, N> undul{};
378 std::array<double, N> Ql{};
379 std::array<double, N> dTl{};
380 std::array<double, N> Tl{};
381 std::array<double, N> pl{};
382 std::array<double, N> ahl{};
383 std::array<double, N> awl{};
384 std::array<double, N> lal{};
385 std::array<double, N> Tml{};
386 std::array<double, N> el{};
387 std::array<double, N> Gn_hl{};
388 std::array<double, N> Ge_hl{};
389 std::array<double, N> Gn_wl{};
390 std::array<double, N> Ge_wl{};
391
392 for (size_t l = 0; l < N; l++)
393 {
394 // transforming ellipsoidal height to orthometric height:
395 // Hortho = -N + Hell
396 auto ix = indx.at(l);
397
398 undul.at(l) = GPT3_grid.at(ix).undulation; // undu: geoid undulation in m
399 double hgt = hell - undul.at(l);
400
401 // pressure, temperature at the height of the grid
402 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 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 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 / 1e3;
408
409 // reduction = stationheight - gridheight
410 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 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 / 1e3;
415
416 // temperature reduction to station height
417 Tl.at(l) = T0 + dTl.at(l) * redh - 273.150;
418
419 // virtual temperature
420 double Tv = T0 * (1 + 0.6077 * Ql.at(l));
421 double c = InsConst::G_NORM * InsConst::dMtr / (InsConst::Rg * Tv);
422
423 // pressure in hPa
424 pl.at(l) = (p0 * exp(-c * redh)) / 100.0;
425
426 // hydrostatic coefficient ah
427 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 / 1e3;
429
430 // wet coefficient aw
431 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 / 1e3;
433
434 // water vapor decrease factor la - added by GP
435 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 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 double e0 = Ql.at(l) * p0 / (0.622 + 0.378 * Ql.at(l)) / 100.0; // on the grid
442
443 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 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 / 1e5;
448 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 / 1e5;
450 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 / 1e5;
452 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 / 1e5;
454 }
455
456 double dnpod1 = fabs(diffpod); // distance nearer point
457 double dnpod2 = 1.0 - dnpod1; // distance to distant point
458 double dnlon1 = fabs(difflon);
459 double dnlon2 = 1.0 - dnlon1;
460
461 // pressure
462 double R1 = dnpod2 * pl[0] + dnpod1 * pl[1];
463 double R2 = dnpod2 * pl[2] + dnpod1 * pl[3];
464 gpt3outputs.p = dnlon2 * R1 + dnlon1 * R2;
465
466 // temperature
467 R1 = dnpod2 * Tl[0] + dnpod1 * Tl[1];
468 R2 = dnpod2 * Tl[2] + dnpod1 * Tl[3];
469 gpt3outputs.T = dnlon2 * R1 + dnlon1 * R2;
470
471 // temperature in degree per km
472 R1 = dnpod2 * dTl[0] + dnpod1 * dTl[1];
473 R2 = dnpod2 * dTl[2] + dnpod1 * dTl[3];
474 gpt3outputs.dT = (dnlon2 * R1 + dnlon1 * R2) * 1e3;
475
476 // water vapor pressure in hPa - changed by GP
477 R1 = dnpod2 * el[0] + dnpod1 * el[1];
478 R2 = dnpod2 * el[2] + dnpod1 * el[3];
479 gpt3outputs.e = dnlon2 * R1 + dnlon1 * R2;
480
481 // hydrostatic
482 R1 = dnpod2 * ahl[0] + dnpod1 * ahl[1];
483 R2 = dnpod2 * ahl[2] + dnpod1 * ahl[3];
484 gpt3outputs.ah = dnlon2 * R1 + dnlon1 * R2;
485
486 // wet
487 R1 = dnpod2 * awl[0] + dnpod1 * awl[1];
488 R2 = dnpod2 * awl[2] + dnpod1 * awl[3];
489 gpt3outputs.aw = dnlon2 * R1 + dnlon1 * R2;
490
491 // undulation
492 R1 = dnpod2 * undul[0] + dnpod1 * undul[1];
493 R2 = dnpod2 * undul[2] + dnpod1 * undul[3];
494 gpt3outputs.undu = dnlon2 * R1 + dnlon1 * R2;
495
496 // water vapor decrease factor la - added by GP
497 R1 = dnpod2 * lal[0] + dnpod1 * lal[1];
498 R2 = dnpod2 * lal[2] + dnpod1 * lal[3];
499 gpt3outputs.la = dnlon2 * R1 + dnlon1 * R2;
500
501 // mean temperature of the water vapor - added by GP
502 R1 = dnpod2 * Tml[0] + dnpod1 * Tml[1];
503 R2 = dnpod2 * Tml[2] + dnpod1 * Tml[3];
504 gpt3outputs.Tm = dnlon2 * R1 + dnlon1 * R2;
505
506 // gradients
507 R1 = dnpod2 * Gn_hl[0] + dnpod1 * Gn_hl[1];
508 R2 = dnpod2 * Gn_hl[2] + dnpod1 * Gn_hl[3];
509 gpt3outputs.Gn_h = dnlon2 * R1 + dnlon1 * R2;
510
511 R1 = dnpod2 * Ge_hl[0] + dnpod1 * Ge_hl[1];
512 R2 = dnpod2 * Ge_hl[2] + dnpod1 * Ge_hl[3];
513 gpt3outputs.Ge_h = dnlon2 * R1 + dnlon1 * R2;
514
515 R1 = dnpod2 * Gn_wl[0] + dnpod1 * Gn_wl[1];
516 R2 = dnpod2 * Gn_wl[2] + dnpod1 * Gn_wl[3];
517 gpt3outputs.Gn_w = dnlon2 * R1 + dnlon1 * R2;
518
519 R1 = dnpod2 * Ge_wl[0] + dnpod1 * Ge_wl[1];
520 R2 = dnpod2 * Ge_wl[2] + dnpod1 * Ge_wl[3];
521 gpt3outputs.Ge_w = dnlon2 * R1 + dnlon1 * R2;
522 }
523
524 return gpt3outputs;
525}
526
527double mjd2doy(const double& mjd)
528{
529 // convert mjd to doy
530 double hour = floor((mjd - floor(mjd)) * 24); // get hours
531 double minu = floor((((mjd - floor(mjd)) * 24) - hour) * 60); // get minutes
532 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 if (sec == 60.0) { minu = minu + 1; }
536 if (minu == 60.0) { hour = hour + 1; }
537
538 // calc jd (yet wrong for hour==24)
539 double jd = mjd + 2400000.5;
540
541 // if hour==24, correct jd and set hour==0
542 if (hour == 24.0) { jd = jd + 1; }
543
544 // integer Julian date
545 double jd_int = floor(jd + 0.5);
546
547 double aa = jd_int + 32044;
548 double bb = floor((4 * aa + 3) / 146097);
549 double cc = aa - floor((bb * 146097) / 4);
550 double dd = floor((4 * cc + 3) / 1461);
551 double ee = cc - floor((1461 * dd) / 4);
552 double mm = floor((5 * ee + 2) / 153);
553
554 auto day = static_cast<int>(ee - floor((153 * mm + 2) / 5) + 1);
555 auto month = static_cast<int>(mm + 3 - 12 * floor(mm / 10));
556 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 double leapYear = 0.0;
560 if ((year % 4 == 0 && year % 100 != 0) || (year % 400 == 0))
561 {
562 leapYear = 1.0;
563 }
564 else
565 {
566 leapYear = 0.0;
567 }
568 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 int sum = 0;
570 for (size_t i = 0; i < static_cast<size_t>(month) - 1; i++)
571 {
572 sum = sum + days[i];
573 }
574 double doy = sum + day;
575
576 if (leapYear == 1 && month > 2)
577 {
578 doy = doy + 1;
579 }
580
581 // Add decimal places
582 doy = doy + mjd - floor(mjd);
583
584 return doy;
585}
586
587double 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
Constants.hpp
Holds all Constants.
GPT2Coeffs.hpp
GPT2 Coefficients.
GPT3Coeffs.hpp
GPT3 Coefficients.
GPT.hpp
GPT2/3 (Global Pressure and Temperature) models.
Math.hpp
Simple Math functions.
NAV::InsConst::Rg
static constexpr double Rg
Universal gas constant in [J/K/mol].
Definition Constants.hpp:230
NAV::InsConst::dMtr
static constexpr double dMtr
Molar mass of dry air in [kg/mol].
Definition Constants.hpp:233
NAV::InsConst::G_NORM
static constexpr double G_NORM
Standard gravity in [m / s^2].
Definition Constants.hpp:40
NAV::internal::GPT3_grid
std::array< GPT3Data, 64800 > GPT3_grid
GPT3 grid.
Definition GPT3Coeffs.cpp:19
NAV::internal::GPT2_grid
std::array< GPT2Data, 64800 > GPT2_grid
GPT2 grid.
Definition GPT2Coeffs.cpp:19
NAV::math::sign
T sign(const T &x, const T &y)
Change the sign of x according to the value of y.
Definition Math.hpp:335
NAV::asknewet
double asknewet(const double &e, const double &Tm, const double &la)
This subroutine determines the zenith wet delay.
Definition GPT.cpp:587
NAV::mjd2doy
double mjd2doy(const double &mjd)
To calculate the day of year.
Definition GPT.cpp:527
NAV::GPT2_param
GPT2output GPT2_param(const double &mjd, const Eigen::Vector3d &lla_pos)
Determine pressure, temperature, temperature lapse rate, mean temperature of the water vapor,...
Definition GPT.cpp:22
NAV::GPT3_param
GPT3output GPT3_param(const double &mjd, const Eigen::Vector3d &lla_pos)
Determine pressure, temperature, temperature lapse rate, mean temperature of the water vapor,...
Definition GPT.cpp:254
NAV::GPT2output
GPT2 output parameters.
Definition GPT.hpp:24
NAV::GPT2output::ah
double ah
hydrostatic mapping function coefficient at zero height (VMF1) (vector of length nstat)
Definition GPT.hpp:30
NAV::GPT2output::e
double e
water vapor pressure in [hPa] (vector of length nstat)
Definition GPT.hpp:29
NAV::GPT2output::Tm
double Tm
mean temperature of the water vapor in [degrees Kelvin] (vector of length nstat)
Definition GPT.hpp:28
NAV::GPT2output::la
double la
water vapor decrease factor (vector of length nstat)
Definition GPT.hpp:32
NAV::GPT2output::undu
double undu
geoid undulation in [m] (vector of length nstat)
Definition GPT.hpp:33
NAV::GPT2output::T
double T
temperature in [degrees Celsius] (vector of length nstat)
Definition GPT.hpp:26
NAV::GPT2output::p
double p
p pressure in [hPa] (vector of length nstat)
Definition GPT.hpp:25
NAV::GPT2output::aw
double aw
wet mapping function coefficient (VMF1) (vector of length nstat)
Definition GPT.hpp:31
NAV::GPT2output::dT
double dT
temperature lapse rate in [degrees per km] (vector of length nstat)
Definition GPT.hpp:27
NAV::GPT3output
GPT3 output parameters.
Definition GPT.hpp:38
NAV::GPT3output::Ge_w
double Ge_w
wet east gradient in [m]
Definition GPT.hpp:42
NAV::GPT3output::Ge_h
double Ge_h
hydrostatic east gradient in [m]
Definition GPT.hpp:40
NAV::GPT3output::Gn_w
double Gn_w
wet north gradient in [m]
Definition GPT.hpp:41
NAV::GPT3output::Gn_h
double Gn_h
hydrostatic north gradient in [m]
Definition GPT.hpp:39