GlobalTropModel.cpp
Go to the documentation of this file.
1 //==============================================================================
2 //
3 // This file is part of GNSSTk, the ARL:UT GNSS Toolkit.
4 //
5 // The GNSSTk is free software; you can redistribute it and/or modify
6 // it under the terms of the GNU Lesser General Public License as published
7 // by the Free Software Foundation; either version 3.0 of the License, or
8 // any later version.
9 //
10 // The GNSSTk is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU Lesser General Public License for more details.
14 //
15 // You should have received a copy of the GNU Lesser General Public
16 // License along with GNSSTk; if not, write to the Free Software Foundation,
17 // Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA
18 //
19 // This software was developed by Applied Research Laboratories at the
20 // University of Texas at Austin.
21 // Copyright 2004-2022, The Board of Regents of The University of Texas System
22 //
23 //==============================================================================
24 
25 //==============================================================================
26 //
27 // This software was developed by Applied Research Laboratories at the
28 // University of Texas at Austin, under contract to an agency or agencies
29 // within the U.S. Department of Defense. The U.S. Government retains all
30 // rights to use, duplicate, distribute, disclose, or release this software.
31 //
32 // Pursuant to DoD Directive 523024
33 //
34 // DISTRIBUTION STATEMENT A: This software has been approved for public
35 // release, distribution is unlimited.
36 //
37 //==============================================================================
38 
39 #include "GlobalTropModel.hpp"
40 #include "MJD.hpp"
41 
42 namespace gnsstk
43 {
44 
45  // Constants for Global mapping functions
46  const double GlobalTropModel::ADryMean[55] = {
47  +1.2517e+02, +8.503e-01, +6.936e-02, -6.760e+00, +1.771e-01,
48  +1.130e-02, +5.963e-01, +1.808e-02, +2.801e-03, -1.414e-03,
49  -1.212e+00, +9.300e-02, +3.683e-03, +1.095e-03, +4.671e-05,
50  +3.959e-01, -3.867e-02, +5.413e-03, -5.289e-04, +3.229e-04,
51  +2.067e-05, +3.000e-01, +2.031e-02, +5.900e-03, +4.573e-04,
52  -7.619e-05, +2.327e-06, +3.845e-06, +1.182e-01, +1.158e-02,
53  +5.445e-03, +6.219e-05, +4.204e-06, -2.093e-06, +1.540e-07,
54  -4.280e-08, -4.751e-01, -3.490e-02, +1.758e-03, +4.019e-04,
55  -2.799e-06, -1.287e-06, +5.468e-07, +7.580e-08, -6.300e-09,
56  -1.160e-01, +8.301e-03, +8.771e-04, +9.955e-05, -1.718e-06,
57  -2.012e-06, +1.170e-08, +1.790e-08, -1.300e-09, +1.000e-10 };
58 
59  const double GlobalTropModel::BDryMean[55] = {
60  +0.000e+00, +0.000e+00, +3.249e-02, +0.000e+00, +3.324e-02,
61  +1.850e-02, +0.000e+00, -1.115e-01, +2.519e-02, +4.923e-03,
62  +0.000e+00, +2.737e-02, +1.595e-02, -7.332e-04, +1.933e-04,
63  +0.000e+00, -4.796e-02, +6.381e-03, -1.599e-04, -3.685e-04,
64  +1.815e-05, +0.000e+00, +7.033e-02, +2.426e-03, -1.111e-03,
65  -1.357e-04, -7.828e-06, +2.547e-06, +0.000e+00, +5.779e-03,
66  +3.133e-03, -5.312e-04, -2.028e-05, +2.323e-07, -9.100e-08,
67  -1.650e-08, +0.000e+00, +3.688e-02, -8.638e-04, -8.514e-05,
68  -2.828e-05, +5.403e-07, +4.390e-07, +1.350e-08, +1.800e-09,
69  +0.000e+00, -2.736e-02, -2.977e-04, +8.113e-05, +2.329e-07,
70  +8.451e-07, +4.490e-08, -8.100e-09, -1.500e-09, +2.000e-10 };
71 
72  const double GlobalTropModel::ADryAmp[55] = {
73  -2.738e-01, -2.837e+00, +1.298e-02, -3.588e-01, +2.413e-02,
74  +3.427e-02, -7.624e-01, +7.272e-02, +2.160e-02, -3.385e-03,
75  +4.424e-01, +3.722e-02, +2.195e-02, -1.503e-03, +2.426e-04,
76  +3.013e-01, +5.762e-02, +1.019e-02, -4.476e-04, +6.790e-05,
77  +3.227e-05, +3.123e-01, -3.535e-02, +4.840e-03, +3.025e-06,
78  -4.363e-05, +2.854e-07, -1.286e-06, -6.725e-01, -3.730e-02,
79  +8.964e-04, +1.399e-04, -3.990e-06, +7.431e-06, -2.796e-07,
80  -1.601e-07, +4.068e-02, -1.352e-02, +7.282e-04, +9.594e-05,
81  +2.070e-06, -9.620e-08, -2.742e-07, -6.370e-08, -6.300e-09,
82  +8.625e-02, -5.971e-03, +4.705e-04, +2.335e-05, +4.226e-06,
83  +2.475e-07, -8.850e-08, -3.600e-08, -2.900e-09, +0.000e+00 };
84 
85  const double GlobalTropModel::BDryAmp[55] = {
86  +0.000e+00, +0.000e+00, -1.136e-01, +0.000e+00, -1.868e-01,
87  -1.399e-02, +0.000e+00, -1.043e-01, +1.175e-02, -2.240e-03,
88  +0.000e+00, -3.222e-02, +1.333e-02, -2.647e-03, -2.316e-05,
89  +0.000e+00, +5.339e-02, +1.107e-02, -3.116e-03, -1.079e-04,
90  -1.299e-05, +0.000e+00, +4.861e-03, +8.891e-03, -6.448e-04,
91  -1.279e-05, +6.358e-06, -1.417e-07, +0.000e+00, +3.041e-02,
92  +1.150e-03, -8.743e-04, -2.781e-05, +6.367e-07, -1.140e-08,
93  -4.200e-08, +0.000e+00, -2.982e-02, -3.000e-03, +1.394e-05,
94  -3.290e-05, -1.705e-07, +7.440e-08, +2.720e-08, -6.600e-09,
95  +0.000e+00, +1.236e-02, -9.981e-04, -3.792e-05, -1.355e-05,
96  +1.162e-06, -1.789e-07, +1.470e-08, -2.400e-09, -4.000e-10 };
97 
98  const double GlobalTropModel::AWetMean[55] = {
99  +5.640e+01, +1.555e+00, -1.011e+00, -3.975e+00, +3.171e-02,
100  +1.065e-01, +6.175e-01, +1.376e-01, +4.229e-02, +3.028e-03,
101  +1.688e+00, -1.692e-01, +5.478e-02, +2.473e-02, +6.059e-04,
102  +2.278e+00, +6.614e-03, -3.505e-04, -6.697e-03, +8.402e-04,
103  +7.033e-04, -3.236e+00, +2.184e-01, -4.611e-02, -1.613e-02,
104  -1.604e-03, +5.420e-05, +7.922e-05, -2.711e-01, -4.406e-01,
105  -3.376e-02, -2.801e-03, -4.090e-04, -2.056e-05, +6.894e-06,
106  +2.317e-06, +1.941e+00, -2.562e-01, +1.598e-02, +5.449e-03,
107  +3.544e-04, +1.148e-05, +7.503e-06, -5.667e-07, -3.660e-08,
108  +8.683e-01, -5.931e-02, -1.864e-03, -1.277e-04, +2.029e-04,
109  +1.269e-05, +1.629e-06, +9.660e-08, -1.015e-07, -5.000e-10 };
110 
111  const double GlobalTropModel::BWetMean[55] = {
112  +0.000e+00, +0.000e+00, +2.592e-01, +0.000e+00, +2.974e-02,
113  -5.471e-01, +0.000e+00, -5.926e-01, -1.030e-01, -1.567e-02,
114  +0.000e+00, +1.710e-01, +9.025e-02, +2.689e-02, +2.243e-03,
115  +0.000e+00, +3.439e-01, +2.402e-02, +5.410e-03, +1.601e-03,
116  +9.669e-05, +0.000e+00, +9.502e-02, -3.063e-02, -1.055e-03,
117  -1.067e-04, -1.130e-04, +2.124e-05, +0.000e+00, -3.129e-01,
118  +8.463e-03, +2.253e-04, +7.413e-05, -9.376e-05, -1.606e-06,
119  +2.060e-06, +0.000e+00, +2.739e-01, +1.167e-03, -2.246e-05,
120  -1.287e-04, -2.438e-05, -7.561e-07, +1.158e-06, +4.950e-08,
121  +0.000e+00, -1.344e-01, +5.342e-03, +3.775e-04, -6.756e-05,
122  -1.686e-06, -1.184e-06, +2.768e-07, +2.730e-08, +5.700e-09 };
123 
124  const double GlobalTropModel::AWetAmp[55] = {
125  +1.023e-01, -2.695e+00, +3.417e-01, -1.405e-01, +3.175e-01,
126  +2.116e-01, +3.536e+00, -1.505e-01, -1.660e-02, +2.967e-02,
127  +3.819e-01, -1.695e-01, -7.444e-02, +7.409e-03, -6.262e-03,
128  -1.836e+00, -1.759e-02, -6.256e-02, -2.371e-03, +7.947e-04,
129  +1.501e-04, -8.603e-01, -1.360e-01, -3.629e-02, -3.706e-03,
130  -2.976e-04, +1.857e-05, +3.021e-05, +2.248e+00, -1.178e-01,
131  +1.255e-02, +1.134e-03, -2.161e-04, -5.817e-06, +8.836e-07,
132  -1.769e-07, +7.313e-01, -1.188e-01, +1.145e-02, +1.011e-03,
133  +1.083e-04, +2.570e-06, -2.140e-06, -5.710e-08, +2.000e-08,
134  -1.632e+00, -6.948e-03, -3.893e-03, +8.592e-04, +7.577e-05,
135  +4.539e-06, -3.852e-07, -2.213e-07, -1.370e-08, +5.800e-09 };
136 
137  const double GlobalTropModel::BWetAmp[55] = {
138  +0.000e+00, +0.000e+00, -8.865e-02, +0.000e+00, -4.309e-01,
139  +6.340e-02, +0.000e+00, +1.162e-01, +6.176e-02, -4.234e-03,
140  +0.000e+00, +2.530e-01, +4.017e-02, -6.204e-03, +4.977e-03,
141  +0.000e+00, -1.737e-01, -5.638e-03, +1.488e-04, +4.857e-04,
142  -1.809e-04, +0.000e+00, -1.514e-01, -1.685e-02, +5.333e-03,
143  -7.611e-05, +2.394e-05, +8.195e-06, +0.000e+00, +9.326e-02,
144  -1.275e-02, -3.071e-04, +5.374e-05, -3.391e-05, -7.436e-06,
145  +6.747e-07, +0.000e+00, -8.637e-02, -3.807e-03, -6.833e-04,
146  -3.861e-05, -2.268e-05, +1.454e-06, +3.860e-07, -1.068e-07,
147  +0.000e+00, -2.658e-02, -1.947e-03, +7.131e-04, -3.506e-05,
148  +1.885e-07, +5.792e-07, +3.990e-08, +2.000e-08, -5.700e-09 };
149 
150  const double GlobalTropModel::Ageoid[55] = {
151  -5.6195e-01,-6.0794e-02,-2.0125e-01,-6.4180e-02,-3.6997e-02,
152  +1.0098e+01,+1.6436e+01,+1.4065e+01,+1.9881e+00,+6.4414e-01,
153  -4.7482e+00,-3.2290e+00,+5.0652e-01,+3.8279e-01,-2.6646e-02,
154  +1.7224e+00,-2.7970e-01,+6.8177e-01,-9.6658e-02,-1.5113e-02,
155  +2.9206e-03,-3.4621e+00,-3.8198e-01,+3.2306e-02,+6.9915e-03,
156  -2.3068e-03,-1.3548e-03,+4.7324e-06,+2.3527e+00,+1.2985e+00,
157  +2.1232e-01,+2.2571e-02,-3.7855e-03,+2.9449e-05,-1.6265e-04,
158  +1.1711e-07,+1.6732e+00,+1.9858e-01,+2.3975e-02,-9.0013e-04,
159  -2.2475e-03,-3.3095e-05,-1.2040e-05,+2.2010e-06,-1.0083e-06,
160  +8.6297e-01,+5.8231e-01,+2.0545e-02,-7.8110e-03,-1.4085e-04,
161  -8.8459e-06,+5.7256e-06,-1.5068e-06,+4.0095e-07,-2.4185e-08 };
162 
163  const double GlobalTropModel::Bgeoid[55] = {
164  +0.0000e+00,+0.0000e+00,-6.5993e-02,+0.0000e+00,+6.5364e-02,
165  -5.8320e+00,+0.0000e+00,+1.6961e+00,-1.3557e+00,+1.2694e+00,
166  +0.0000e+00,-2.9310e+00,+9.4805e-01,-7.6243e-02,+4.1076e-02,
167  +0.0000e+00,-5.1808e-01,-3.4583e-01,-4.3632e-02,+2.2101e-03,
168  -1.0663e-02,+0.0000e+00,+1.0927e-01,-2.9463e-01,+1.4371e-03,
169  -1.1452e-02,-2.8156e-03,-3.5330e-04,+0.0000e+00,+4.4049e-01,
170  +5.5653e-02,-2.0396e-02,-1.7312e-03,+3.5805e-05,+7.2682e-05,
171  +2.2535e-06,+0.0000e+00,+1.9502e-02,+2.7919e-02,-8.1812e-03,
172  +4.4540e-04,+8.8663e-05,+5.5596e-05,+2.4826e-06,+1.0279e-06,
173  +0.0000e+00,+6.0529e-02,-3.5824e-02,-5.1367e-03,+3.0119e-05,
174  -2.9911e-05,+1.9844e-05,-1.2349e-06,-7.6756e-09,+5.0100e-08 };
175 
176  const double GlobalTropModel::APressMean[55] = {
177  +1.0108e+03,+8.4886e+00,+1.4799e+00,-1.3897e+01,+3.7516e-03,
178  -1.4936e-01,+1.2232e+01,-7.6615e-01,-6.7699e-02,+8.1002e-03,
179  -1.5874e+01,+3.6614e-01,-6.7807e-02,-3.6309e-03,+5.9966e-04,
180  +4.8163e+00,-3.7363e-01,-7.2071e-02,+1.9998e-03,-6.2385e-04,
181  -3.7916e-04,+4.7609e+00,-3.9534e-01,+8.6667e-03,+1.1569e-02,
182  +1.1441e-03,-1.4193e-04,-8.5723e-05,+6.5008e-01,-5.0889e-01,
183  -1.5754e-02,-2.8305e-03,+5.7458e-04,+3.2577e-05,-9.6052e-06,
184  -2.7974e-06,+1.3530e+00,-2.7271e-01,-3.0276e-04,+3.6286e-03,
185  -2.0398e-04,+1.5846e-05,-7.7787e-06,+1.1210e-06,+9.9020e-08,
186  +5.5046e-01,-2.7312e-01,+3.2532e-03,-2.4277e-03,+1.1596e-04,
187  +2.6421e-07,-1.3263e-06,+2.7322e-07,+1.4058e-07,+4.9414e-09 };
188 
189  const double GlobalTropModel::BPressMean[55] = {
190  +0.0000e+00,+0.0000e+00,-1.2878e+00,+0.0000e+00,+7.0444e-01,
191  +3.3222e-01,+0.0000e+00,-2.9636e-01,+7.2248e-03,+7.9655e-03,
192  +0.0000e+00,+1.0854e+00,+1.1145e-02,-3.6513e-02,+3.1527e-03,
193  +0.0000e+00,-4.8434e-01,+5.2023e-02,-1.3091e-02,+1.8515e-03,
194  +1.5422e-04,+0.0000e+00,+6.8298e-01,+2.5261e-03,-9.9703e-04,
195  -1.0829e-03,+1.7688e-04,-3.1418e-05,+0.0000e+00,-3.7018e-01,
196  +4.3234e-02,+7.2559e-03,+3.1516e-04,+2.0024e-05,-8.0581e-06,
197  -2.3653e-06,+0.0000e+00,+1.0298e-01,-1.5086e-02,+5.6186e-03,
198  +3.2613e-05,+4.0567e-05,-1.3925e-06,-3.6219e-07,-2.0176e-08,
199  +0.0000e+00,-1.8364e-01,+1.8508e-02,+7.5016e-04,-9.6139e-05,
200  -3.1995e-06,+1.3868e-07,-1.9486e-07,+3.0165e-10,-6.4376e-10 };
201 
202  const double GlobalTropModel::APressAmp[55] = {
203  -1.0444e-01,+1.6618e-01,-6.3974e-02,+1.0922e+00,+5.7472e-01,
204  -3.0277e-01,-3.5087e+00,+7.1264e-03,-1.4030e-01,+3.7050e-02,
205  +4.0208e-01,-3.0431e-01,-1.3292e-01,+4.6746e-03,-1.5902e-04,
206  +2.8624e+00,-3.9315e-01,-6.4371e-02,+1.6444e-02,-2.3403e-03,
207  +4.2127e-05,+1.9945e+00,-6.0907e-01,-3.5386e-02,-1.0910e-03,
208  -1.2799e-04,+4.0970e-05,+2.2131e-05,-5.3292e-01,-2.9765e-01,
209  -3.2877e-02,+1.7691e-03,+5.9692e-05,+3.1725e-05,+2.0741e-05,
210  -3.7622e-07,+2.6372e+00,-3.1165e-01,+1.6439e-02,+2.1633e-04,
211  +1.7485e-04,+2.1587e-05,+6.1064e-06,-1.3755e-08,-7.8748e-08,
212  -5.9152e-01,-1.7676e-01,+8.1807e-03,+1.0445e-03,+2.3432e-04,
213  +9.3421e-06,+2.8104e-06,-1.5788e-07,-3.0648e-08,+2.6421e-10 };
214 
215  const double GlobalTropModel::BPressAmp[55] = {
216  +0.0000e+00,+0.0000e+00,+9.3340e-01,+0.0000e+00,+8.2346e-01,
217  +2.2082e-01,+0.0000e+00,+9.6177e-01,-1.5650e-02,+1.2708e-03,
218  +0.0000e+00,-3.9913e-01,+2.8020e-02,+2.8334e-02,+8.5980e-04,
219  +0.0000e+00,+3.0545e-01,-2.1691e-02,+6.4067e-04,-3.6528e-05,
220  -1.1166e-04,+0.0000e+00,-7.6974e-02,-1.8986e-02,+5.6896e-03,
221  -2.4159e-04,-2.3033e-04,-9.6783e-06,+0.0000e+00,-1.0218e-01,
222  -1.3916e-02,-4.1025e-03,-5.1340e-05,-7.0114e-05,-3.3152e-07,
223  +1.6901e-06,+0.0000e+00,-1.2422e-02,+2.5072e-03,+1.1205e-03,
224  -1.3034e-04,-2.3971e-05,-2.6622e-06,+5.7852e-07,+4.5847e-08,
225  +0.0000e+00,+4.4777e-02,-3.0421e-03,+2.6062e-05,-7.2421e-05,
226  +1.9119e-06,+3.9236e-07,+2.2390e-07,+2.9765e-09,-4.6452e-09 };
227 
228  const double GlobalTropModel::ATempMean[55] = {
229  +1.6257e+01,+2.1224e+00,+9.2569e-01,-2.5974e+01,+1.4510e+00,
230  +9.2468e-02,-5.3192e-01,+2.1094e-01,-6.9210e-02,-3.4060e-02,
231  -4.6569e+00,+2.6385e-01,-3.6093e-02,+1.0198e-02,-1.8783e-03,
232  +7.4983e-01,+1.1741e-01,+3.9940e-02,+5.1348e-03,+5.9111e-03,
233  +8.6133e-06,+6.3057e-01,+1.5203e-01,+3.9702e-02,+4.6334e-03,
234  +2.4406e-04,+1.5189e-04,+1.9581e-07,+5.4414e-01,+3.5722e-01,
235  +5.2763e-02,+4.1147e-03,-2.7239e-04,-5.9957e-05,+1.6394e-06,
236  -7.3045e-07,-2.9394e+00,+5.5579e-02,+1.8852e-02,+3.4272e-03,
237  -2.3193e-05,-2.9349e-05,+3.6397e-07,+2.0490e-06,-6.4719e-08,
238  -5.2225e-01,+2.0799e-01,+1.3477e-03,+3.1613e-04,-2.2285e-04,
239  -1.8137e-05,-1.5177e-07,+6.1343e-07,+7.8566e-08,+1.0749e-09 };
240 
241  const double GlobalTropModel::BTempMean[55] = {
242  +0.0000e+00,+0.0000e+00,+1.0210e+00,+0.0000e+00,+6.0194e-01,
243  +1.2292e-01,+0.0000e+00,-4.2184e-01,+1.8230e-01,+4.2329e-02,
244  +0.0000e+00,+9.3312e-02,+9.5346e-02,-1.9724e-03,+5.8776e-03,
245  +0.0000e+00,-2.0940e-01,+3.4199e-02,-5.7672e-03,-2.1590e-03,
246  +5.6815e-04,+0.0000e+00,+2.2858e-01,+1.2283e-02,-9.3679e-03,
247  -1.4233e-03,-1.5962e-04,+4.0160e-05,+0.0000e+00,+3.6353e-02,
248  -9.4263e-04,-3.6762e-03,+5.8608e-05,-2.6391e-05,+3.2095e-06,
249  -1.1605e-06,+0.0000e+00,+1.6306e-01,+1.3293e-02,-1.1395e-03,
250  +5.1097e-05,+3.3977e-05,+7.6449e-06,-1.7602e-07,-7.6558e-08,
251  +0.0000e+00,-4.5415e-02,-1.8027e-02,+3.6561e-04,-1.1274e-04,
252  +1.3047e-05,+2.0001e-06,-1.5152e-07,-2.7807e-08,+7.7491e-09 };
253 
254  const double GlobalTropModel::ATempAmp[55] = {
255  -1.8654e+00,-9.0041e+00,-1.2974e-01,-3.6053e+00,+2.0284e-02,
256  +2.1872e-01,-1.3015e+00,+4.0355e-01,+2.2216e-01,-4.0605e-03,
257  +1.9623e+00,+4.2887e-01,+2.1437e-01,-1.0061e-02,-1.1368e-03,
258  -6.9235e-02,+5.6758e-01,+1.1917e-01,-7.0765e-03,+3.0017e-04,
259  +3.0601e-04,+1.6559e+00,+2.0722e-01,+6.0013e-02,+1.7023e-04,
260  -9.2424e-04,+1.1269e-05,-6.9911e-06,-2.0886e+00,-6.7879e-02,
261  -8.5922e-04,-1.6087e-03,-4.5549e-05,+3.3178e-05,-6.1715e-06,
262  -1.4446e-06,-3.7210e-01,+1.5775e-01,-1.7827e-03,-4.4396e-04,
263  +2.2844e-04,-1.1215e-05,-2.1120e-06,-9.6421e-07,-1.4170e-08,
264  +7.8720e-01,-4.4238e-02,-1.5120e-03,-9.4119e-04,+4.0645e-06,
265  -4.9253e-06,-1.8656e-06,-4.0736e-07,-4.9594e-08,+1.6134e-09 };
266 
267  const double GlobalTropModel::BTempAmp[55] = {
268  +0.0000e+00,+0.0000e+00,-8.9895e-01,+0.0000e+00,-1.0790e+00,
269  -1.2699e-01,+0.0000e+00,-5.9033e-01,+3.4865e-02,-3.2614e-02,
270  +0.0000e+00,-2.4310e-02,+1.5607e-02,-2.9833e-02,-5.9048e-03,
271  +0.0000e+00,+2.8383e-01,+4.0509e-02,-1.8834e-02,-1.2654e-03,
272  -1.3794e-04,+0.0000e+00,+1.3306e-01,+3.4960e-02,-3.6799e-03,
273  -3.5626e-04,+1.4814e-04,+3.7932e-06,+0.0000e+00,+2.0801e-01,
274  +6.5640e-03,-3.4893e-03,-2.7395e-04,+7.4296e-05,-7.9927e-06,
275  -1.0277e-06,+0.0000e+00,+3.6515e-02,-7.4319e-03,-6.2873e-04,
276  -8.2461e-05,+3.1095e-05,-5.3860e-07,-1.2055e-07,-1.1517e-07,
277  +0.0000e+00,+3.1404e-02,+1.5580e-02,-1.1428e-03,+3.3529e-05,
278  +1.0387e-05,-1.9378e-06,-2.7327e-07,+7.5833e-09,-9.2323e-09 };
279 
280  const double GlobalTropModel::Factorial[19] = {
281  1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800,
282  479001600, 6227020800, 87178291200, 1307674368000, 20922789888000,
283  355687428096000, 6402373705728000 };
284 
285  const double GlobalTropModel::HEIGHT_LIMIT = 44243.;
286 
287  GlobalTropModel :: GlobalTropModel()
288  : validCoeff(false), validHeight(false), validLat(false),
289  validLon(false), validDay(false), height(0.0), latitude(0.0),
290  longitude(0.0), dayfactor(0.0), undul(0.0)
291  {
292  // yes setting everything to 0 is the same as IEEE 0.0
293  memset(P, 0, sizeof(P));
294  memset(aP, 0, sizeof(aP));
295  memset(bP, 0, sizeof(bP));
296  TropModel::humid = 50.0;
297  valid = false;
298  }
299 
300 
301  double GlobalTropModel::correction(double elevation) const
302  {
303  try { testValidity(); }
304  catch(InvalidTropModel& e) { GNSSTK_RETHROW(e); }
305 
306  // Global mapping functions good down to 3 degrees of elevation
307  if(elevation < 3.0) { return 0.0; }
308 
309  double map_dry(GlobalTropModel::dry_mapping_function(elevation));
310  double map_wet(GlobalTropModel::wet_mapping_function(elevation));
311 
312  // Compute total tropospheric delay
313  double tropDelay((GlobalTropModel::dry_zenith_delay() * map_dry) +
314  (GlobalTropModel::wet_zenith_delay() * map_wet));
315 
316  return tropDelay;
317 
318  } // end GlobalTropModel::correction(elevation)
319 
320 
321  double GlobalTropModel::correction(const Position& RX, const Position& SV)
322  {
323  try {
324  double p;
325  p = RX.getAltitude(); if(p != height) setReceiverHeight(p);
326  p = RX.getGeodeticLatitude(); if(p != latitude) setReceiverLatitude(p);
327  p = RX.getLongitude(); if(p != longitude) setReceiverLongitude(p);
328  }
329  catch(GeometryException& e) {
330  validHeight = validLat = valid = false;
331  GNSSTK_RETHROW(e);
332  }
333 
334  try { testValidity(); }
335  catch(InvalidTropModel& e) { GNSSTK_RETHROW(e); }
336 
337  double c;
338  try {
339  c = GlobalTropModel::correction(RX.elevationGeodetic(SV));
340  }
341  catch(InvalidTropModel& e) { GNSSTK_RETHROW(e); }
342 
343  return c;
344 
345  } // end GlobalTropModel::correction(RX,SV)
346 
347 
348  double GlobalTropModel::dry_zenith_delay() const
349  {
350  try { testValidity(); } catch(InvalidTropModel& e) { GNSSTK_RETHROW(e); }
351  return SaasDryDelay(press,latitude,height);
352 
353  } // end GlobalTropModel::dry_zenith_delay()
354 
355 
356  double GlobalTropModel::wet_zenith_delay() const
357  {
358  double T = temp + CELSIUS_TO_KELVIN;
359  double pwv = 0.01 * humid * ::exp(-37.2465 + (0.213166-0.000256908*T)*T);
360  return (0.0122 + 0.00943 * pwv);
361  }
362 
363 
364  double GlobalTropModel::dry_mapping_function(double elevation) const
365  {
366  try { testValidity(); } catch(InvalidTropModel& e) { GNSSTK_RETHROW(e); }
367  if(elevation < 3.0) { return 0.0; }
368 
369  double clat = ::cos(latitude*DEG_TO_RAD);
370  double phh, c11h, c10h;
371 
372  static const double bh = 0.0029;
373  static const double c0h = 0.062;
374  if(latitude < 0) {
375  phh = PI;
376  c11h = 0.007;
377  c10h = 0.002;
378  }
379  else {
380  phh = 0.0;
381  c11h = 0.005;
382  c10h = 0.001;
383  }
384  double ch = c0h + ((::cos(dayfactor + phh)+1.0)*c11h/2.0 + c10h)*(1.0-clat);
385 
386  double am(0.0), aa(0.0);
387  for(int i=0; i<55; i++) {
388  am += (ADryMean[i]*aP[i] + BDryMean[i]*bP[i]) * 1.0e-5;
389  aa += (ADryAmp[i]*aP[i] + BDryAmp[i]*bP[i]) * 1.0e-5;
390  }
391  double ah = am + aa*::cos(dayfactor);
392 
393  double sine = ::sin(elevation*DEG_TO_RAD);
394  //std::cout << "sine " << std::fixed << std::setprecision(16) << sine
395  // << std::endl;
396  //std::cout << "ah bh ch " << std::fixed << std::setprecision(16)
397  // << ah << " " << bh << " " << ch << std::endl;
398  double map = (1.0 + ah/(1.0 + bh/(1.0 + ch)))
399  / (sine + ah/(sine + bh/(sine + ch)));
400  //std::cout << "map0 " << std::fixed << std::setprecision(16) << map
401  // << std::endl;
402 
403  // height correction
404  static const double a_ht = 2.53e-5;
405  static const double b_ht = 5.49e-3;
406  static const double c_ht = 1.14e-3;
407 
408  map += ( (1.0/sine) - (1.0 + a_ht/(1.0 + b_ht/(1.0 + c_ht)))
409  / (sine + a_ht/(sine + b_ht/(sine + c_ht)))
410  ) * (height/1000.0);
411 
412  return map;
413 
414  } // end GlobalTropModel::dry_mapping_function()
415 
416  /* Computing the derivative - do it numerically instead
417  * note that sin[elev] = cos[zenith]
418  * f(1) a
419  * map = ------------ where f(x) = x + ---------
420  * f(sin[elev]) b
421  * x + -----
422  * x + c
423  *
424  * a (x+c) x(x^2 + xc + b) + a(x+c)
425  * f(x) = x + ------------ = ------------------------
426  * x^2 + xc + b x^2 + xc + b
427  *
428  * x^3 + x^2 c + x(a+b) + ac
429  * = -------------------------
430  * x^2 + xc + b
431  *
432  * so -f(1)f'(x) -map(x)f'(x) N'D-D'N
433  * map'(x) = ---------- = -------------, where f'(x) = -------
434  * f^2(x) f(x) D^2
435  *
436  * [3x^2+2xc+a+b][x^2+xc+b]-[2x+c][x^3+x^2c+x(a+b)+ac]
437  * f'(x) = x' --------------------------------------------------------------
438  * [x^2 + xc + b]^2
439  *
440  * [3x^4 + x^3(5c) + x^2(a+4b+2c^2) + x(a+b+2bc) + ab+b^2]
441  * + [-2x^4 - x^3(3c) - x^2(2a+2b+c^2) - x(c(3a+b)) - ac^2]
442  * = x' ----------------------------------------------------------
443  * x^4 + x^3(2c) + x^2(2b+c^2) + x(2bc) + b^2
444  *
445  * x^4 + x^3(2c) + x^2(-a+2b+c^2) + x(a+b+3bc+3ac) + ab+b^2-ac^2
446  * = x' -------------------------------------------------------------
447  * x^4 + x^3(2c) + x^2(2b+c^2) + x(2bc) + b^2
448  *
449  */
450  double GlobalTropModel::wet_mapping_function(double elevation) const
451  {
452  try { testValidity(); }
453  catch(InvalidTropModel& e) { GNSSTK_RETHROW(e); }
454 
455  if(elevation < 3.0) { return 0.0; }
456 
457  static const double bw = 0.00146;
458  static const double cw = 0.04391;
459 
460  double am(0.0), aa(0.0);
461  for(int i=0; i<55; i++) {
462  am += (AWetMean[i]*aP[i] + BWetMean[i]*bP[i]) * 1.0e-5;
463  aa += (AWetAmp[i]*aP[i] + BWetAmp[i]*bP[i]) * 1.0e-5;
464  }
465  double aw = am + aa*::cos(dayfactor);
466 
467  double sine = ::sin(elevation*DEG_TO_RAD);
468  //std::cout << "sine " << std::fixed << std::setprecision(16) << sine
469  // << std::endl;
470  //std::cout << "aw bw cw " << std::fixed << std::setprecision(16)
471  // << aw << " " << bw << " " << cw << std::endl;
472  double f1(1.0 + aw/(1.0 + bw/(1.0 + cw)));
473  double f(sine + aw/(sine + bw/(sine + cw)));
474  double map = f1/f;
475 
488  return map;
489 
490  } // end GlobalTropModel::wet_mapping_function()
491 
492 
493  void GlobalTropModel::getGPT(double& P, double& T, double& U)
494  {
495  try { testValidity(); }
496  catch(InvalidTropModel& e) { GNSSTK_RETHROW(e); }
497 
498  int i;
499 
500  // undulation and orthometric height
501  U = 0.0;
502  for(i=0; i<55; i++) U += (Ageoid[i]*aP[i] + Bgeoid[i]*bP[i]);
503  double orthoht(height - U);
504  if(orthoht > HEIGHT_LIMIT)
505  {
506  InvalidTropModel exc("Invalid Global trop model: Rx Height exceeds limit");
507  GNSSTK_THROW(exc);
508  }
509 
510  // press at geoid
511  double am(0.0),aa(0.0),v0;
512  for(i=0; i<55; i++) {
513  am += (APressMean[i]*aP[i] + BPressMean[i]*bP[i]);
514  aa += (APressAmp[i]*aP[i] + BPressAmp[i]*bP[i]);
515  }
516  v0 = am + aa * ::cos(dayfactor);
517 
518  // pressure at height
519  // @note this implies any orthoht > 1/2.26e-5 == 44247.78m is invalid!
520  P = v0 * ::pow(1.0-2.26e-5*orthoht,5.225);
521 
522  // temper on geoid
523  am = aa = 0.0;
524  for(i=0; i<55; i++) {
525  am += (ATempMean[i]*aP[i] + BTempMean[i]*bP[i]);
526  aa += (ATempAmp[i]*aP[i] + BTempAmp[i]*bP[i]);
527  }
528  v0 = am + aa * ::cos(dayfactor);
529 
530  // temp at height
531  T = v0 - 6.5e-3 * orthoht;
532  }
533 
534 
535  void GlobalTropModel::setReceiverHeight(const double& ht)
536  {
537  if(height != ht) {
538  height = ht;
539  validHeight = true;
540  validCoeff = false;
541  setValid(); // calls updateGTMCoeff()
542  }
543  }
544 
545 
546  void GlobalTropModel::setReceiverLatitude(const double& lat)
547  {
548  if(latitude != lat) {
549  latitude = lat;
550  validLat = true;
551  validCoeff = false;
552  setValid(); // calls updateGTMCoeff()
553  }
554  }
555 
556 
557  void GlobalTropModel::setReceiverLongitude(const double& lon)
558  {
559  if(longitude != lon) {
560  longitude = lon;
561  validLon = true;
562  validCoeff = false;
563  setValid(); // calls updateGTMCoeff()
564  }
565  }
566 
567 
568  void GlobalTropModel::setTime(const double& mjd)
569  {
570  double df(TWO_PI*(mjd - 44266.0)/365.25); // -44239 + 1 - 28
571  if(df != dayfactor) {
572  dayfactor = df;
573  validDay = true;
574  validCoeff = false;
575  setValid();
576  }
577  }
578 
579 
580  void GlobalTropModel::setTime(const CommonTime& time)
581  {
582  double mjd = static_cast<MJD>(time).mjd;
583  setTime(mjd);
584  }
585 
586 
587  void GlobalTropModel::setDayOfYear(const int& doy)
588  {
589  double mjd = 44266.0 + (double)doy; // year doesn't matter
590  setTime(mjd);
591  }
592 
593 
594  void GlobalTropModel::setParameters(const CommonTime& time, const Position& rxPos)
595  {
596  validDay = validHeight = validLat = validLon = validCoeff = false;
597  setTime(time);
598  setReceiverHeight(rxPos.getHeight());
599  setReceiverLatitude(rxPos.getGeodeticLatitude());
600  setReceiverLongitude(rxPos.getLongitude());
601 
602  setValid(); // calls updateGTMCoeff()
603  }
604 
605 
606  void GlobalTropModel::updateGTMCoeff()
607  {
608  if(!validLon || !validLat) return;
609 
610  // compute Legendre functions and spherical harmonics
611  int i,j,k;
612  double sinlat(::sin(latitude*DEG_TO_RAD));
613  for(i=0; i<=9; i++) {
614  for(j=0; j<=i; j++) {
615  int ir((i-j)/2);
616  double sum(0.0),term;
617  int sign(1);
618  for(k=0; k<=ir; k++) {
619  term = sign*(((Factorial[2*i-2*k] / Factorial[k]) / Factorial[i-k])
620  / Factorial[i-j-2*k]) * ::pow(sinlat,i-j-2*k);
621  sum += term;
622  sign = -sign;
623  }
624  P[i][j] = (1.0/(::pow(2.0,i)) * ::sqrt(::pow(1.0-sinlat*sinlat,j)) * sum);
625  }
626  }
627 
628  // spherical harmonics
629  double rlon(longitude*DEG_TO_RAD);
630  i = 0;
631  for(j=0; j<=9; j++) {
632  for(k=0; k<=j; k++) {
633  aP[i] = P[j][k] * ::cos(k*rlon);
634  bP[i] = P[j][k] * ::sin(k*rlon);
635  i++;
636  }
637  }
638 
639  }
640 
641 
642  void GlobalTropModel::testValidity() const
643  {
644  if(!valid) {
645  if(!validLat)
646  GNSSTK_THROW(InvalidTropModel("Invalid Global trop model: Rx Latitude"));
647  if(!validHeight)
648  GNSSTK_THROW(InvalidTropModel("Invalid Global trop model: Rx Height"));
649  if(!validDay)
650  GNSSTK_THROW(InvalidTropModel("Invalid Global trop model: day of year"));
651 
652  GNSSTK_THROW(InvalidTropModel("Valid flag corrupted in Global trop model"));
653  }
654  }
655 
656 } // end namespace gnsstk
gnsstk::GlobalTropModel::latitude
double latitude
Definition: GlobalTropModel.hpp:311
gnsstk::GlobalTropModel::testValidity
void testValidity() const
Definition: GlobalTropModel.cpp:642
gnsstk::GlobalTropModel::setReceiverHeight
virtual void setReceiverHeight(const double &ht)
Definition: GlobalTropModel.cpp:535
gnsstk::GlobalTropModel::P
double P[10][10]
Definition: GlobalTropModel.hpp:312
example6.mjd
mjd
Definition: example6.py:102
gnsstk::GlobalTropModel::setDayOfYear
virtual void setDayOfYear(const int &doy)
Definition: GlobalTropModel.cpp:587
gnsstk::GlobalTropModel::setParameters
virtual void setParameters(const CommonTime &time, const Position &rxPos)
Definition: GlobalTropModel.cpp:594
gnsstk::GlobalTropModel::validCoeff
bool validCoeff
Definition: GlobalTropModel.hpp:313
gnsstk::GlobalTropModel::BWetMean
static const double BWetMean[55]
Definition: GlobalTropModel.hpp:290
gnsstk::GlobalTropModel::BPressAmp
static const double BPressAmp[55]
Definition: GlobalTropModel.hpp:299
gnsstk::GlobalTropModel::setReceiverLatitude
virtual void setReceiverLatitude(const double &lat)
Definition: GlobalTropModel.cpp:546
gnsstk::GlobalTropModel::APressMean
static const double APressMean[55]
Definition: GlobalTropModel.hpp:296
gnsstk::GlobalTropModel::dayfactor
double dayfactor
Definition: GlobalTropModel.hpp:311
gnsstk::sum
T sum(const ConstVectorBase< T, BaseClass > &l)
Definition: VectorBaseOperators.hpp:84
gnsstk::Position::getAltitude
double getAltitude() const noexcept
return height above ellipsoid (meters)
Definition: Position.hpp:469
gnsstk::TropModel::SaasDryDelay
double SaasDryDelay(const double pr, const double lat, const double ht) const
Definition: TropModel.hpp:262
gnsstk::GlobalTropModel::AWetAmp
static const double AWetAmp[55]
Definition: GlobalTropModel.hpp:291
gnsstk::TWO_PI
const double TWO_PI
GPS value of PI*2.
Definition: GNSSconstants.hpp:68
gnsstk::GlobalTropModel::BTempAmp
static const double BTempAmp[55]
Definition: GlobalTropModel.hpp:303
gnsstk::PI
const double PI
GPS value of PI; also specified by GAL.
Definition: GNSSconstants.hpp:62
gnsstk::GlobalTropModel::setValid
void setValid()
Definition: GlobalTropModel.hpp:327
gnsstk::GlobalTropModel::validDay
bool validDay
Definition: GlobalTropModel.hpp:313
gnsstk::GlobalTropModel::validHeight
bool validHeight
Definition: GlobalTropModel.hpp:313
gnsstk::GlobalTropModel::HEIGHT_LIMIT
static const GNSSTK_EXPORT double HEIGHT_LIMIT
Model is limited in height, at this value, in m.
Definition: GlobalTropModel.hpp:309
gnsstk
For Sinex::InputHistory.
Definition: BasicFramework.cpp:50
MJD.hpp
gnsstk::GlobalTropModel::correction
virtual double correction(double elevation) const
Definition: GlobalTropModel.cpp:301
std::sin
double sin(gnsstk::Angle x)
Definition: Angle.hpp:144
gnsstk::TropModel::CELSIUS_TO_KELVIN
static const GNSSTK_EXPORT double CELSIUS_TO_KELVIN
Definition: TropModel.hpp:108
gnsstk::GlobalTropModel::ADryAmp
static const double ADryAmp[55]
Definition: GlobalTropModel.hpp:287
gnsstk::TropModel::humid
double humid
latest value of relative humidity (percent)
Definition: TropModel.hpp:282
gnsstk::GlobalTropModel::aP
double aP[55]
Definition: GlobalTropModel.hpp:312
gnsstk::Position::getHeight
double getHeight() const noexcept
return height above ellipsoid (meters)
Definition: Position.hpp:474
gnsstk::GlobalTropModel::validLon
bool validLon
Definition: GlobalTropModel.hpp:313
gnsstk::GlobalTropModel::Bgeoid
static const double Bgeoid[55]
Definition: GlobalTropModel.hpp:295
gnsstk::GlobalTropModel::BWetAmp
static const double BWetAmp[55]
Definition: GlobalTropModel.hpp:292
example4.time
time
Definition: example4.py:103
gnsstk::Position::getGeodeticLatitude
double getGeodeticLatitude() const noexcept
return geodetic latitude (deg N)
Definition: Position.hpp:454
gnsstk::CommonTime
Definition: CommonTime.hpp:84
gnsstk::GlobalTropModel::BTempMean
static const double BTempMean[55]
Definition: GlobalTropModel.hpp:301
gnsstk::GlobalTropModel::height
double height
Definition: GlobalTropModel.hpp:311
gnsstk::GlobalTropModel::updateGTMCoeff
void updateGTMCoeff()
Update coefficients when latitude and/or longitude changes.
Definition: GlobalTropModel.cpp:606
gnsstk::TropModel::valid
bool valid
true only if current model parameters are valid
Definition: TropModel.hpp:279
gnsstk::GlobalTropModel::bP
double bP[55]
Definition: GlobalTropModel.hpp:312
std::cos
double cos(gnsstk::Angle x)
Definition: Angle.hpp:146
gnsstk::TropModel::press
double press
latest value of pressure (millibars)
Definition: TropModel.hpp:281
GNSSTK_RETHROW
#define GNSSTK_RETHROW(exc)
Definition: Exception.hpp:369
gnsstk::GlobalTropModel::BDryMean
static const double BDryMean[55]
Definition: GlobalTropModel.hpp:286
gnsstk::GlobalTropModel::setTime
void setTime(const double &mjd)
Definition: GlobalTropModel.cpp:568
gnsstk::GlobalTropModel::ATempAmp
static const double ATempAmp[55]
Definition: GlobalTropModel.hpp:302
gnsstk::GlobalTropModel::APressAmp
static const double APressAmp[55]
Definition: GlobalTropModel.hpp:298
gnsstk::GlobalTropModel::dry_zenith_delay
virtual double dry_zenith_delay() const
Definition: GlobalTropModel.cpp:348
gnsstk::GlobalTropModel::longitude
double longitude
Definition: GlobalTropModel.hpp:311
gnsstk::GlobalTropModel::BPressMean
static const double BPressMean[55]
Definition: GlobalTropModel.hpp:297
gnsstk::GlobalTropModel::Ageoid
static const double Ageoid[55]
Definition: GlobalTropModel.hpp:294
gnsstk::GlobalTropModel::Factorial
static const double Factorial[19]
Definition: GlobalTropModel.hpp:305
gnsstk::Position::getLongitude
double getLongitude() const noexcept
return longitude (deg E) (either geocentric or geodetic)
Definition: Position.hpp:464
gnsstk::GlobalTropModel::GlobalTropModel
GlobalTropModel()
Default constructor.
Definition: GlobalTropModel.cpp:287
gnsstk::GlobalTropModel::validLat
bool validLat
Definition: GlobalTropModel.hpp:313
gnsstk::GlobalTropModel::wet_mapping_function
virtual double wet_mapping_function(double elevation) const
Definition: GlobalTropModel.cpp:450
gnsstk::Position
Definition: Position.hpp:136
gnsstk::MJD
Definition: MJD.hpp:54
GNSSTK_THROW
#define GNSSTK_THROW(exc)
Definition: Exception.hpp:366
gnsstk::GlobalTropModel::AWetMean
static const double AWetMean[55]
Definition: GlobalTropModel.hpp:289
GlobalTropModel.hpp
gnsstk::TropModel::temp
double temp
latest value of temperature (kelvin or celsius)
Definition: TropModel.hpp:280
gnsstk::GlobalTropModel::getGPT
void getGPT(double &P, double &T, double &U)
Definition: GlobalTropModel.cpp:493
gnsstk::Position::elevationGeodetic
double elevationGeodetic(const Position &Target) const
Definition: Position.cpp:1331
gnsstk::GlobalTropModel::setReceiverLongitude
virtual void setReceiverLongitude(const double &lon)
Definition: GlobalTropModel.cpp:557
gnsstk::GlobalTropModel::BDryAmp
static const double BDryAmp[55]
Definition: GlobalTropModel.hpp:288
gnsstk::GlobalTropModel::ATempMean
static const double ATempMean[55]
Definition: GlobalTropModel.hpp:300
gnsstk::GlobalTropModel::dry_mapping_function
virtual double dry_mapping_function(double elevation) const
Definition: GlobalTropModel.cpp:364
gnsstk::GlobalTropModel::wet_zenith_delay
virtual double wet_zenith_delay() const
Definition: GlobalTropModel.cpp:356
gnsstk::GlobalTropModel::ADryMean
static const double ADryMean[55]
Definition: GlobalTropModel.hpp:285
gnsstk::DEG_TO_RAD
static const double DEG_TO_RAD
Conversion Factor from degrees to radians (units: degrees^-1)
Definition: GNSSconstants.hpp:76

gnsstk
Author(s):
autogenerated on Wed Oct 25 2023 02:40:39