+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//
+// Calculate the WGS84 gravitation value in ECEF frame. Pass in the ECEF position
+// via the position parameter. The J2Gravity value returned is in ECEF frame,
+// and therefore may need to be expressed (transformed) in another frame,
+// depending on how it is used. See Stevens and Lewis eqn. 1.4-16.
+
+FGColumnVector3 FGInertial::GetGravityJ2(FGColumnVector3 position) const
+{
+ FGColumnVector3 J2Gravity;
+
+ // Gravitation accel
+ double r = position.Magnitude();
+ double lat = FDMExec->GetPropagate()->GetLatitude();
+ double sinLat = sin(lat);
+
+ double adivr = a/r;
+ double preCommon = 1.5*J2*adivr*adivr;
+ double xy = 1.0 - 5.0*(sinLat*sinLat);
+ double z = 3.0 - 5.0*(sinLat*sinLat);
+ double GMOverr2 = GM/(r*r);
+
+ J2Gravity(1) = -GMOverr2 * ((1.0 + (preCommon * xy)) * position(eX)/r);
+ J2Gravity(2) = -GMOverr2 * ((1.0 + (preCommon * xy)) * position(eY)/r);
+ J2Gravity(3) = -GMOverr2 * ((1.0 + (preCommon * z)) * position(eZ)/r);
+
+ return J2Gravity;
+}
+