Skip to content

Commit d009f0d

Browse files
xanthospapxanthospap
committed
gravity field error
1 parent a1e64a2 commit d009f0d

File tree

6 files changed

+120
-84
lines changed

6 files changed

+120
-84
lines changed

‎bin/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
# bin/CMakeLists.txt
22

33
set(BIN_SOURCES
4-
equations_of_motion.cpp
4+
#equations_of_motion.cpp
55
test_deriv_costg.cpp
66
)
77

‎bin/config.yaml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ planetary-ephemeris:
77

88
gravity:
99
# Filename of gravity field model in ICGEM/gfc format
10-
model: data/EIGEN-GRGS.RL04.MEAN-FIELD.gfc
10+
model: data/EIGEN6-C4.gfc
1111
degree: 180
1212
order: 180
1313

‎bin/test_deriv_costg.cpp‎

Lines changed: 61 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -66,12 +66,11 @@ int gcrf2ecef(const dso::MjdEpoch &tai, dso::EopSeries &eops,
6666
}
6767

6868
int deriv(double tsec, Eigen::Ref<const Eigen::VectorXd> y0,
69-
const dso::IntegrationParameters *params,
69+
dso::IntegrationParameters *params,
7070
Eigen::Ref<Eigen::VectorXd> y) noexcept {
7171

7272
/* current time in TAI */
73-
dso::MjdEpoch t = params->t0();
74-
t.add_seconds(dso::FractionalSeconds(tsec));
73+
dso::MjdEpoch t = params->t0().add_seconds(dso::FractionalSeconds(tsec));
7574

7675
/* current time in GPST (debugging) */
7776
const auto tgps = t.tai2gps();
@@ -82,18 +81,21 @@ int deriv(double tsec, Eigen::Ref<const Eigen::VectorXd> y0,
8281
/* Celestial to Terrestrial Matrix */
8382
dso::EopRecord eopr;
8483
Eigen::Matrix<double, 3, 3> R, dRdt;
85-
if (gcrf2ecef(dso::MjdEpoch(t), *(params->eops()), R, dRdt, fargs, eopr)) {
84+
if (gcrf2ecef(dso::MjdEpoch(t), params->eops(), R, dRdt, fargs, eopr)) {
8685
return 8;
8786
}
8887

88+
/* satellite position in ECEF */
89+
const Eigen::VectorXd r_ecef = R * y0.segment<3>(0);
90+
8991
/* Gravity field stokes coefficients */
90-
auto stokes = params->grav();
91-
stokes.C(0, 0) = stokes.C(1, 0) = stokes.C(1, 1) = 0e0;
92+
dso::StokesCoeffs stokes(params->earth_gravity());
93+
stokes.C(0, 0) = 0e0;
94+
stokes.C(1, 0) = stokes.C(1, 1) = 0e0;
9295
stokes.S(1, 1) = 0e0;
9396

9497
/* compute acceleration for given epoch/position (ECEF) */
9598
[[maybe_unused]] Eigen::Matrix<double, 3, 3> g;
96-
const Eigen::VectorXd r_ecef = R * y0.segment<3>(0);
9799
Eigen::Matrix<double, 3, 1> acc_grav;
98100
if (dso::sh2gradient_cunningham(stokes, r_ecef, acc_grav, g,
99101
stokes.max_degree(), stokes.max_order(), -1,
@@ -142,12 +144,12 @@ int deriv(double tsec, Eigen::Ref<const Eigen::VectorXd> y0,
142144

143145
/* Solid Earth Tide (ITRF) */
144146
Eigen::Matrix<double, 3, 1> acc_set;
145-
params->mse_tide->stokes_coeffs(t.tai2tt(), t.tai2ut1(eopr.dut()),
146-
R * rtb_moon, R * rtb_sun.segment<3>(0),
147-
fargs);
148-
if (dso::sh2gradient_cunningham(params->mse_tide->stokes_coeffs(), r_ecef,
149-
acc_set, g, -1, -1, -1, -1, &(params->tw()),
150-
&(params->tm()))) {
147+
params->solid_earth_tide()->stokes_coeffs(t.tai2tt(), t.tai2ut1(eopr.dut()),
148+
R * rtb_moon,
149+
R * rtb_sun.segment<3>(0), fargs);
150+
if (dso::sh2gradient_cunningham(params->solid_earth_tide()->stokes_coeffs(),
151+
r_ecef, acc_set, g, -1, -1, -1, -1,
152+
&(params->tw()), &(params->tm()))) {
151153
fprintf(stderr, "ERROR Failed computing acceleration/gradient\n");
152154
return 1;
153155
}
@@ -180,20 +182,22 @@ int deriv(double tsec, Eigen::Ref<const Eigen::VectorXd> y0,
180182

181183
/* Ocean Pole Tide */
182184
Eigen::Matrix<double, 3, 1> acc_otp;
183-
if (params->mop_tide->stokes_coeffs(t.tai2tt(), eopr.xp(), eopr.yp(),
184-
params->mop_tide->max_degree(),
185-
params->mop_tide->max_order())) {
185+
if (params->ocean_pole_tide()->stokes_coeffs(
186+
t.tai2tt(), eopr.xp(), eopr.yp(),
187+
params->ocean_pole_tide()->max_degree(),
188+
params->ocean_pole_tide()->max_order())) {
186189
fprintf(stderr, "ERROR Failed computing Stokes Coefficients\n");
187190
return 1;
188191
}
189-
params->mop_tide->stokes_coeffs().C(0, 0) =
190-
params->mop_tide->stokes_coeffs().C(1, 0) =
191-
params->mop_tide->stokes_coeffs().C(1, 1) = 0e0;
192-
params->mop_tide->stokes_coeffs().S(1, 1) = 0e0;
193-
if (dso::sh2gradient_cunningham(params->mop_tide->stokes_coeffs(), r_ecef,
194-
acc_otp, g, params->mop_tide->max_degree(),
195-
params->mop_tide->max_order(), -1, -1,
196-
&(params->tw()), &(params->tm()))) {
192+
params->ocean_pole_tide()->stokes_coeffs().C(0, 0) =
193+
params->ocean_pole_tide()->stokes_coeffs().C(1, 0) =
194+
params->ocean_pole_tide()->stokes_coeffs().C(1, 1) = 0e0;
195+
params->ocean_pole_tide()->stokes_coeffs().S(1, 1) = 0e0;
196+
if (dso::sh2gradient_cunningham(params->ocean_pole_tide()->stokes_coeffs(),
197+
r_ecef, acc_otp, g,
198+
params->ocean_pole_tide()->max_degree(),
199+
params->ocean_pole_tide()->max_order(), -1,
200+
-1, &(params->tw()), &(params->tm()))) {
197201
fprintf(stderr, "ERROR Failed computing acceleration/gradient\n");
198202
return 1;
199203
}
@@ -204,20 +208,21 @@ int deriv(double tsec, Eigen::Ref<const Eigen::VectorXd> y0,
204208
y.segment<3>(0) += acc;
205209
}
206210

207-
/* Dealiasing */
211+
/* Dealiasing
212+
* WARNING !!
213+
* Here, we are storing the Stokes coefficients in stokes (created earlier).
214+
*/
208215
Eigen::Matrix<double, 3, 1> acc_das;
209-
if (params->mdealias->coefficients_at(
216+
if (params->dealias()->coefficients_at(
210217
dso::from_mjdepoch<dso::nanoseconds>(t.tai2tt()), stokes)) {
211218
fprintf(stderr, "Failed interpolating dealiasing coefficients\n");
212219
return 1;
213220
}
214-
stokes.C(0, 0) = 0e0;
215-
stokes.C(1, 0) = stokes.C(1, 1) = 0e0;
221+
stokes.C(0, 0) = stokes.C(1, 0) = stokes.C(1, 1) = 0e0;
216222
stokes.S(1, 1) = 0e0;
217223
if (dso::sh2gradient_cunningham(stokes, r_ecef, acc_das, g,
218-
params->mdealias_maxdegree,
219-
params->mdealias_maxorder, -1, -1,
220-
&(params->tw()), &(params->tm()))) {
224+
stokes.max_degree(), stokes.max_order(), -1,
225+
-1, &(params->tw()), &(params->tm()))) {
221226
fprintf(stderr, "ERROR Failed computing acceleration/gradient\n");
222227
return 1;
223228
}
@@ -231,17 +236,19 @@ int deriv(double tsec, Eigen::Ref<const Eigen::VectorXd> y0,
231236
/* Atmospheric Tide */
232237
Eigen::Matrix<double, 3, 1> acc_atm;
233238
/* compute Stokes coeffs (for atm. tides) */
234-
params->mat_tide->stokes_coeffs(t.tai2tt(), t.tai2ut1(eopr.dut()), fargs);
239+
params->atmospheric_tide()->stokes_coeffs(t.tai2tt(), t.tai2ut1(eopr.dut()),
240+
fargs);
235241
/* for the test, degree one coefficients are not taken into account */
236-
params->mat_tide->stokes_coeffs().C(0, 0) =
237-
params->mat_tide->stokes_coeffs().C(1, 0) =
238-
params->mat_tide->stokes_coeffs().C(1, 1) = 0e0;
239-
params->mat_tide->stokes_coeffs().S(1, 1) = 0e0;
242+
params->atmospheric_tide()->stokes_coeffs().C(0, 0) =
243+
params->atmospheric_tide()->stokes_coeffs().C(1, 0) =
244+
params->atmospheric_tide()->stokes_coeffs().C(1, 1) = 0e0;
245+
params->atmospheric_tide()->stokes_coeffs().S(1, 1) = 0e0;
240246
/* compute acceleration for given epoch/position (ITRF) */
241-
if (dso::sh2gradient_cunningham(params->mat_tide->stokes_coeffs(), r_ecef,
242-
acc_atm, g, params->matm_maxdegree,
243-
params->matm_maxorder, -1, -1,
244-
&(params->tw()), &(params->tm()))) {
247+
if (dso::sh2gradient_cunningham(params->atmospheric_tide()->stokes_coeffs(),
248+
r_ecef, acc_atm, g,
249+
params->atmospheric_tide()->max_degree(),
250+
params->atmospheric_tide()->max_order(), -1,
251+
-1, &(params->tw()), &(params->tm()))) {
245252
fprintf(stderr, "ERROR Failed computing acceleration/gradient\n");
246253
return 1;
247254
}
@@ -255,16 +262,16 @@ int deriv(double tsec, Eigen::Ref<const Eigen::VectorXd> y0,
255262
/* Ocean Tide */
256263
Eigen::Matrix<double, 3, 1> acc_ot;
257264
/* compute Stokes coeffs (for ocean tides) */
258-
params->moc_tide->stokes_coeffs(t.tai2tt(), t.tai2ut1(eopr.dut()), fargs);
265+
params->ocean_tide()->stokes_coeffs(t.tai2tt(), t.tai2ut1(eopr.dut()), fargs);
259266
/* for the test, degree one coefficients are not taken into account */
260-
params->moc_tide->stokes_coeffs().C(0, 0) =
261-
params->moc_tide->stokes_coeffs().C(1, 0) =
262-
params->moc_tide->stokes_coeffs().C(1, 1) = 0e0;
263-
params->moc_tide->stokes_coeffs().S(1, 1) = 0e0;
267+
params->ocean_tide()->stokes_coeffs().C(0, 0) =
268+
params->ocean_tide()->stokes_coeffs().C(1, 0) =
269+
params->ocean_tide()->stokes_coeffs().C(1, 1) = 0e0;
270+
params->ocean_tide()->stokes_coeffs().S(1, 1) = 0e0;
264271
/* compute acceleration for given epoch/position (ITRF) */
265-
if (dso::sh2gradient_cunningham(params->moc_tide->stokes_coeffs(), r_ecef,
266-
acc_ot, g, params->moc_tide->max_degree(),
267-
params->moc_tide->max_order(), -1, -1,
272+
if (dso::sh2gradient_cunningham(params->ocean_tide()->stokes_coeffs(), r_ecef,
273+
acc_ot, g, params->ocean_tide()->max_degree(),
274+
params->ocean_tide()->max_order(), -1, -1,
268275
&(params->tw()), &(params->tm()))) {
269276
fprintf(stderr, "ERROR Failed computing acceleration/gradient\n");
270277
return 1;
@@ -315,16 +322,15 @@ int main(int argc, char *argv[]) {
315322
sv = dso::sp3::SatelliteId{argv[3]};
316323
}
317324

318-
/* get starting epoch in TAI */
325+
/* get starting epoch in TT */
319326
auto start_t = sp3.start_epoch();
320327
if (!std::strcmp(sp3.time_sys(), "GPS")) {
321328
start_t = start_t.gps2tai();
322329
}
323330

324-
auto t2 = start_t;
325-
t2.add_seconds(dso::seconds(5 * 86400));
326-
dso::IntegrationParameters params =
327-
dso::IntegrationParameters::from_config(argv[1], start_t, t2);
331+
auto t2 = start_t.add_seconds(dso::seconds(5 * 86400));
332+
dso::IntegrationParameters params = dso::IntegrationParameters::from_config(
333+
argv[1], dso::MjdEpoch(start_t), dso::MjdEpoch(t2));
328334

329335
///* EOPs */
330336
// dso::EopSeries eop;
@@ -476,7 +482,8 @@ int main(int argc, char *argv[]) {
476482
block_tai = block_tai.gps2tai();
477483
}
478484
/* GCRF to ITRF rotation matrix */
479-
if (gcrf2ecef(dso::MjdEpoch(block_tai), eop, R, dRdt, fargs, eopr)) {
485+
if (gcrf2ecef(dso::MjdEpoch(block_tai), params.eops(), R, dRdt, fargs,
486+
eopr)) {
480487
return 8;
481488
}
482489
/* get state for current epoch ITRF */

‎include/integration_parameters.hpp‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -113,6 +113,26 @@ class IntegrationParameters {
113113

114114
const MjdEpoch &t0() const noexcept { return mtai0; }
115115
MjdEpoch &t0() noexcept { return mtai0; }
116+
const EopSeries &eops() const noexcept { return meops; }
117+
EopSeries &eops() noexcept { return meops; }
118+
const StokesCoeffs &earth_gravity() const noexcept { return mgrav; }
119+
StokesCoeffs &earth_gravity() noexcept { return mgrav; }
120+
const SolidEarthTide *solid_earth_tide() const noexcept { return mse_tide; }
121+
SolidEarthTide *solid_earth_tide() noexcept { return mse_tide; }
122+
const OceanTide *ocean_tide() const noexcept { return moc_tide; }
123+
OceanTide *ocean_tide() noexcept { return moc_tide; }
124+
const PoleTide *pole_tide() const noexcept { return mep_tide; }
125+
PoleTide *pole_tide() noexcept { return mep_tide; }
126+
const OceanPoleTide *ocean_pole_tide() const noexcept { return mop_tide; }
127+
OceanPoleTide *ocean_pole_tide() noexcept { return mop_tide; }
128+
const AtmosphericTide *atmospheric_tide() const noexcept { return mat_tide; }
129+
AtmosphericTide *atmospheric_tide() noexcept { return mat_tide; }
130+
const Aod1bDataStream<AOD1BCoefficientType::GLO> *dealias() const noexcept {
131+
return mdealias;
132+
}
133+
Aod1bDataStream<AOD1BCoefficientType::GLO> *dealias() noexcept {
134+
return mdealias;
135+
}
116136

117137
/** @brief Construct instance from a YAML configuration file. */
118138
static IntegrationParameters from_config(const char *config_fn,

‎src/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,4 +2,5 @@ target_sources(integrator
22
PRIVATE
33
${CMAKE_SOURCE_DIR}/src/dp86co.cpp
44
${CMAKE_SOURCE_DIR}/src/hinit853.cpp
5+
${CMAKE_SOURCE_DIR}/src/integration_parameters_from_config.cpp
56
)

0 commit comments

Comments
 (0)