1 //==============================================================================
2 //
3 // This file is part of GPSTk, the GPS Toolkit.
4 //
5 // The GPSTk 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 GPSTk 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 GPSTk; 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-2020, 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 <iostream>
40 #include "TestUtil.hpp"
41 #include "OrbitEph.hpp"
42 #include "CivilTime.hpp"
43 #include "GPSWeekSecond.hpp"
44
45 using namespace std;
46
47 /** Threshold for how much different our velocities can be between
48 * being computed directly via svXvt and computed via differencing
49 * svXvt positions over time. */
50 double velDiffThresh = 0.0008;
51
52
53 class OrbitEph_T
54 {
55 public:
56 unsigned testSvXvt();
57 };
58
59
60 unsigned OrbitEph_T ::
testSvXvt()61 testSvXvt()
62 {
63 TUDEF("OrbitEph", "svXvt");
64 // Hard code orbital parameters mostly so we can copy and paste
65 // the data into other similar tests with minimal changes.
66 gpstk::OrbitEph oe;
67 oe.Cuc = -.324845314026e-05;
68 oe.Cus = .101532787085e-04;
69 oe.Crc = .168968750000e+03;
70 oe.Crs = -.646250000000e+02;
71 oe.Cic = .320374965668e-06;
72 oe.Cis = .117346644402e-06;
73 oe.M0 = -.136404614938e+01;
74 oe.dn = .489591822036e-08;
75 oe.dndot = 0;
76 oe.ecc = .146582192974e-01;
77 oe.A = .515359719276e+04 * .515359719276e+04;
78 oe.Adot = 0;
79 oe.OMEGA0 = -.296605403382e+01;
80 oe.i0 = .941587707856e+00;
81 oe.w = -.224753761329e+01;
82 oe.OMEGAdot = -.804390648956e-08;
83 oe.idot = .789318592573e-10;
84 oe.ctToc = gpstk::CivilTime(2015,7,19,1,59,28.0,gpstk::TimeSystem::GPS);
85 oe.af0 = .579084269702e-03;
86 oe.af1 = .227373675443e-11;
87 oe.af2 = .000000000000e+00;
88 oe.dataLoadedFlag = true;
89 oe.satID = gpstk::SatID(2, gpstk::SatelliteSystem::GPS);
90 oe.ctToe = gpstk::GPSWeekSecond(1854,.716800000000e+04);
91 // iode .700000000000e+01
92 // codes on L2 .100000000000e+01
93 // L2 P data .000000000000e+00
94 // sv accuracy .240000000000e+01
95 // sv health .000000000000e+00
96 // tgd -.204890966415e-07
97 // iodc .700000000000e+01
98 // xmit time .360000000000e+04
99 // fit int .400000000000e+01
100 bool testFailed = false;
101 try
102 {
103 // first compute Xvt
104 static const unsigned SECONDS = 7200;
105 gpstk::Xvt zeroth_array[SECONDS];
106 for (unsigned ii = 0; ii < SECONDS; ii++)
107 {
108 zeroth_array[ii] = oe.svXvt(oe.ctToc + ii);
109 }
110 // then compute first derivative of position, i.e. velocity
111 gpstk::Triple deriv[SECONDS];
112 double h = 1; // time step size in seconds
113 for (unsigned ii = 0; ii < SECONDS; ii++)
114 {
115 if (ii == 0)
116 {
117 deriv[ii] = (1/h)*(-1.5*zeroth_array[ii].getPos() +
118 2.0*zeroth_array[ii+1].getPos() -
119 0.5*zeroth_array[ii+2].getPos());
120 }
121 else if ((ii == 1) || (ii == (SECONDS-2)))
122 {
123 deriv[ii] = (1/h)*(-0.5*zeroth_array[ii-1].getPos() +
124 0.5*zeroth_array[ii+1].getPos());
125 }
126 else if (ii == (SECONDS-1))
127 {
128 deriv[ii] = (1/h)*(0.5*zeroth_array[ii-2].getPos() -
129 2.0*zeroth_array[ii-1].getPos() +
130 1.5*zeroth_array[ii].getPos());
131 }
132 else
133 {
134 deriv[ii] = (1/h)*((1.0/12.0)*zeroth_array[ii-2].getPos() -
135 (2.0/3.0)*zeroth_array[ii-1].getPos() +
136 (2.0/3.0)*zeroth_array[ii+1].getPos() -
137 (1.0/12.0)*zeroth_array[ii+2].getPos());
138 }
139 }
140 // then check the difference between derived and computed velocity
141 for (unsigned ii = 0; ii < SECONDS; ii++)
142 {
143 double derivedMag = deriv[ii].mag();
144 double computedMag = zeroth_array[ii].getVel().mag();
145 // If you want to print the data e.g. to plot, uncomment
146 // this stream output statement and comment out tbe break
147 // statement below.
148 // Just don't check it in that way, please.
149 // cerr << ii << " " << (computedMag - derivedMag) << endl;
150 if (fabs(computedMag - derivedMag) > velDiffThresh)
151 {
152 // no sense in printing 7200 success/fail messages.
153 testFailed = true;
154 break;
155 }
156 }
157 if (testFailed)
158 {
159 TUFAIL("computed velocity is significantly different from derived"
160 " velocity");
161 }
162 else
163 {
164 TUPASS("velocity check");
165 }
166 }
167 catch (gpstk::Exception& exc)
168 {
169 cerr << exc;
170 TUFAIL("Exception");
171 }
172 catch (...)
173 {
174 TUFAIL("Exception");
175 }
176 TURETURN();
177 }
178
179
main()180 int main()
181 {
182 unsigned total = 0;
183 OrbitEph_T testClass;
184 total += testClass.testSvXvt();
185
186 cout << "Total Failures for " << __FILE__ << ": " << total << endl;
187 return total;
188 }
189