Skip to content

Commit 6f9a8f7

Browse files
committed
Simplifed code by removing R() from the code. Resolved some extra declarations.
Rotated omega (earth angular rate, not IMU angular velocities) to the body frame for simplification Cannot rename "D_t_t" to say "D_p_p" or "D_v_t" to "D_v_p" because of NavState derivatives?
1 parent 50e8c8c commit 6f9a8f7

File tree

1 file changed

+18
-17
lines changed

1 file changed

+18
-17
lines changed

gtsam/navigation/NavState.cpp

Lines changed: 18 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -436,6 +436,13 @@ NavState NavState::update(const Vector3& b_acceleration, const Vector3& b_omega,
436436
}
437437

438438
//------------------------------------------------------------------------------
439+
440+
// Because our navigation frames are placed on a spinning Earth, we experience two apparent forces on our inertials
441+
// Let Omega be the Earth's rotation rate in the navigation frame
442+
// Coriolis acceleration = -2 * (omega X n_v)
443+
// Centrifugal acceleration (secondOrder) = -omega X (omega X n_t)
444+
// We would also experience a rotation of (omega*dt) over time - so, counteract by compensating rotation by (-omega * dt)
445+
// Integrate centrifugal & coriolis accelerations to yield position, velocity perturbations
439446
Vector9 NavState::coriolis(double dt, const Vector3& omega, bool secondOrder,
440447
OptionalJacobian<9, 9> H) const {
441448
Rot3 nRb = R_;
@@ -444,48 +451,42 @@ Vector9 NavState::coriolis(double dt, const Vector3& omega, bool secondOrder,
444451
const double dt2 = dt * dt;
445452
const Vector3 omega_cross_vel = omega.cross(n_v);
446453

447-
// Because our navigation frames are placed on a spinning Earth, we experience two apparent forces on our inertials
448-
// Let Omega be the Earth's rotation rate in the navigation frame
449-
// Coriolis acceleration = -2 * (omega X n_v)
450-
// Centrifugal acceleration (secondOrder) = -omega X (omega X n_t)
451-
// We would also experience a rotation of (omega*dt) over time - so, counteract by compensating rotation by (-omega * dt)
452-
// Integrate centrifugal & coriolis accelerations to yield position, velocity perturbations
453-
454454
Vector9 n_xi;
455455
// Coriolis (first order) acceleration corrections
456456
dR(n_xi) << ((-dt) * omega);
457457
dP(n_xi) << ((-dt2) * omega_cross_vel); // NOTE(luca): we got rid of the 2 wrt INS paper
458458
dV(n_xi) << ((-2.0 * dt) * omega_cross_vel);
459459

460460
// Centrifugal (second order) acceleration corrections
461-
Matrix3 H2_wrt_t; // To store Jacobian (if needed/desired)
461+
Matrix3 D_c2_nt; // To store Jacobian (if needed/desired)
462462
if (secondOrder) {
463-
const Vector3 omega_cross2_t = doubleCross(omega, n_t, nullptr, H ? &H2_wrt_t : nullptr);
463+
const Vector3 omega_cross2_t = doubleCross(omega, n_t, nullptr, H ? &D_c2_nt : nullptr);
464464
dP(n_xi) -= (0.5 * dt2) * omega_cross2_t;
465465
dV(n_xi) -= dt * omega_cross2_t;
466466
}
467467

468468
// Transform correction from navigation frame -> body frame and get Jacobians
469469
Vector9 xi;
470-
Matrix3 D_dR_R, D_dP_R, D_dV_R, D_body_nav;
471-
dR(xi) = nRb.unrotate(dR(n_xi), H ? &D_dR_R : 0, H ? &D_body_nav : 0);
470+
Matrix3 D_dR_R, D_dP_R, D_dV_R;
471+
dR(xi) = nRb.unrotate(dR(n_xi), H ? &D_dR_R : 0);
472472
dP(xi) = nRb.unrotate(dP(n_xi), H ? &D_dP_R : 0);
473473
dV(xi) = nRb.unrotate(dV(n_xi), H ? &D_dV_R : 0);
474474

475475
// Assemble Jacobians
476476
if (H) {
477477
H->setZero();
478-
const Matrix3 Omega = skewSymmetric(omega);
479-
const Matrix3 D_cross_state = Omega * R();
478+
const Vector3 omega_b = nRb.unrotate(omega);
479+
const Matrix3 D_c1_b = skewSymmetric(omega_b);
480480
H->setZero();
481481
D_R_R(H) << D_dR_R;
482-
D_t_v(H) << D_body_nav * (-dt2) * D_cross_state;
482+
D_t_v(H) << (-dt2) * D_c1_b;
483483
D_t_R(H) << D_dP_R;
484-
D_v_v(H) << D_body_nav * (-2.0 * dt) * D_cross_state;
484+
D_v_v(H) << (-2.0 * dt) * D_c1_b;
485485
D_v_R(H) << D_dV_R;
486486
if (secondOrder) {
487-
D_t_t(H) -= (0.5 * dt2) * D_body_nav * (H2_wrt_t * R());
488-
D_v_t(H) -= dt * D_body_nav * (H2_wrt_t * R());
487+
Matrix3 D_c2_b = D_c1_b * D_c1_b;
488+
D_t_t(H) -= (0.5 * dt2) * D_c2_b;
489+
D_v_t(H) -= dt * D_c2_b;
489490
}
490491
}
491492
return xi;

0 commit comments

Comments
 (0)