Commit 0c944b19 authored by Bryan Cazabonne's avatar Bryan Cazabonne

Fixed calculation of CR3BP constants.

Fixes #744
parent d403cd18
Pipeline #855 passed with stage
in 28 minutes and 33 seconds
......@@ -21,6 +21,9 @@
</properties>
<body>
<release version="11.0" date="TBD" description="TBD">
<action dev="bryan" type="fix" issue="744">
Fixed calculation of CR3BP constants.
</action>
<action dev="bryan" type="update" issue="743">
Updated JUnit version to 4.13.1.
</action>
......
......@@ -18,6 +18,8 @@ package org.orekit.bodies;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.propagation.numerical.cr3bp.CR3BPConstants;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
/**
* Factory class creating predefined CR3BP system using CR3BPSystem class. For example, Earth-Moon CR3BP
......@@ -37,19 +39,23 @@ public class CR3BPFactory {
}
/** Get the Sun-Jupiter CR3BP singleton bodies pair.
* @param date date
* @param timeScale time scale
* @return Sun-Jupiter CR3BP system
*/
@DefaultDataContext
public static CR3BPSystem getSunJupiterCR3BP() {
return getSystem(CelestialBodyFactory.getSun(), CelestialBodyFactory.getJupiter(), CR3BPConstants.JUPITER_SEMI_MAJOR_AXIS);
public static CR3BPSystem getSunJupiterCR3BP(final AbsoluteDate date, final TimeScale timeScale) {
return getSystem(CelestialBodyFactory.getSun(), CelestialBodyFactory.getJupiter(), CR3BPConstants.getJupiterSemiMajorAxis(date, timeScale));
}
/** Get the Sun-Earth CR3BP singleton bodies pair.
* @param date date
* @param timeScale time scale
* @return Sun-Earth CR3BP system
*/
@DefaultDataContext
public static CR3BPSystem getSunEarthCR3BP() {
return getSystem(CelestialBodyFactory.getSun(), CelestialBodyFactory.getEarthMoonBarycenter(), CR3BPConstants.EARTH_MOON_BARYCENTER_SEMI_MAJOR_AXIS);
public static CR3BPSystem getSunEarthCR3BP(final AbsoluteDate date, final TimeScale timeScale) {
return getSystem(CelestialBodyFactory.getSun(), CelestialBodyFactory.getEarthMoonBarycenter(), CR3BPConstants.getEarthMoonBarycenterSemiMajorAxis(date, timeScale));
}
/** Get the Earth-Moon CR3BP singleton bodies pair.
......@@ -57,7 +63,7 @@ public class CR3BPFactory {
*/
@DefaultDataContext
public static CR3BPSystem getEarthMoonCR3BP() {
return getSystem(CelestialBodyFactory.getEarth(), CelestialBodyFactory.getMoon(), CR3BPConstants.MOON_SEMI_MAJOR_AXIS);
return getSystem(CelestialBodyFactory.getEarth(), CelestialBodyFactory.getMoon(), CR3BPConstants.getMoonSemiMajorAxis());
}
/** Get the corresponding CR3BP System.
......
......@@ -16,26 +16,74 @@
*/
package org.orekit.propagation.numerical.cr3bp;
import java.util.Calendar;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.DateComponents;
import org.orekit.time.TimeScale;
import org.orekit.utils.Constants;
/**
* Set of useful physical CR3BP constants using JPL data.
* @author Vincent Mouraux
* @since 11.0
*/
public class CR3BPConstants {
/** Private constructor.
* <p>This class is a utility class, it should neither have a public
* nor a default constructor. This private constructor prevents
* the compiler from generating one automatically.</p>
*/
private CR3BPConstants() {
}
public interface CR3BPConstants {
/**
* Get the Moon semi-major axis.
* @return the Moon semi-major axis in meters
*/
public static double getMoonSemiMajorAxis() {
return 384400000.0;
}
/** Moon semi-major axis = 384 400 000 m. */
double CENTURY = (Calendar.getInstance().get(Calendar.YEAR) - 2000.0) / 100.0;
/**
* Get the Earth-Moon barycenter semi-major axis.
* @param date date
* @param timeScale time scale
* @return the Earth-Moon barycenter semi-major axis in meters
*/
public static double getEarthMoonBarycenterSemiMajorAxis(final AbsoluteDate date,
final TimeScale timeScale) {
// Century
final double century = getCentury(date, timeScale);
// Return the Earth - Moon barycenter semi-major axis
return (1.00000261 + 0.00000562 * century) * Constants.IAU_2012_ASTRONOMICAL_UNIT;
}
/** Moon semi-major axis in meters. */
double MOON_SEMI_MAJOR_AXIS = 384400000.0;
/**
* Get the Jupiter semi-major axis.
* @param date date
* @param timeScale time scale
* @return the Jupiter semi-major axis in meters
*/
public static double getJupiterSemiMajorAxis(final AbsoluteDate date,
final TimeScale timeScale) {
// Century
final double century = getCentury(date, timeScale);
// Return the Jupiter semi-major axis
return (5.20288700 - 0.00011607 * century) * Constants.IAU_2012_ASTRONOMICAL_UNIT;
}
/** Earth-Moon barycenter semi-major axis in meters. */
double EARTH_MOON_BARYCENTER_SEMI_MAJOR_AXIS = (1.00000261 + 0.00000562 * CENTURY) * Constants.IAU_2012_ASTRONOMICAL_UNIT;
/**
* Get the current century as a floating value.
* @param date date
* @param timeScale time scale
* @return the current century as a floating value
*/
private static double getCentury(final AbsoluteDate date,
final TimeScale timeScale) {
// Get the date component
final DateComponents dc = date.getComponents(timeScale).getDate();
// Return the current century as a floating value
return 0.01 * (dc.getYear() - 2000.0);
}
/** Jupiter semi-major axis in meters. */
double JUPITER_SEMI_MAJOR_AXIS = (5.20288700 - 0.00011607 * CENTURY) * Constants.IAU_2012_ASTRONOMICAL_UNIT;
}
......@@ -20,12 +20,17 @@ import org.junit.Assert;
import org.junit.Before;
import org.junit.Test;
import org.orekit.Utils;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
public class CR3BPFactoryTest {
@Test
public void getSunJupiterCR3BP() {
CR3BPSystem sunJupiterCR3BP = CR3BPFactory.getSunJupiterCR3BP();
AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
TimeScale timeScale = TimeScalesFactory.getUTC();
CR3BPSystem sunJupiterCR3BP = CR3BPFactory.getSunJupiterCR3BP(date, timeScale);
Assert.assertNotNull(sunJupiterCR3BP);
}
......@@ -37,7 +42,9 @@ public class CR3BPFactoryTest {
@Test
public void getSunEarthCR3BP() {
CR3BPSystem sunEarthCR3BP = CR3BPFactory.getSunEarthCR3BP();
AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
TimeScale timeScale = TimeScalesFactory.getUTC();
CR3BPSystem sunEarthCR3BP = CR3BPFactory.getSunEarthCR3BP(date, timeScale);
Assert.assertNotNull(sunEarthCR3BP);
}
......
......@@ -26,6 +26,7 @@ import org.orekit.Utils;
import org.orekit.frames.Frame;
import org.orekit.frames.Transform;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.AbsolutePVCoordinates;
import org.orekit.utils.LagrangianPoints;
......@@ -40,7 +41,9 @@ public class CR3BPSystemTest {
@Test
public void testCR3BPSystem() {
final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP();
AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
TimeScale timeScale = TimeScalesFactory.getUTC();
final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP(date, timeScale);
final double lDim = syst.getDdim();
Assert.assertNotNull(lDim);
......@@ -53,31 +56,41 @@ public class CR3BPSystemTest {
@Test
public void testGetRotatingFrame() {
final Frame baryFrame = CR3BPFactory.getSunEarthCR3BP().getRotatingFrame();
AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
TimeScale timeScale = TimeScalesFactory.getUTC();
final Frame baryFrame = CR3BPFactory.getSunEarthCR3BP(date, timeScale).getRotatingFrame();
Assert.assertNotNull(baryFrame);
}
@Test
public void testGetPrimary() {
final CelestialBody primaryBody = CR3BPFactory.getSunEarthCR3BP().getPrimary();
AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
TimeScale timeScale = TimeScalesFactory.getUTC();
final CelestialBody primaryBody = CR3BPFactory.getSunEarthCR3BP(date, timeScale).getPrimary();
Assert.assertNotNull(primaryBody);
}
@Test
public void testGetSecondary() {
final CelestialBody secondaryBody = CR3BPFactory.getSunEarthCR3BP().getSecondary();
AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
TimeScale timeScale = TimeScalesFactory.getUTC();
final CelestialBody secondaryBody = CR3BPFactory.getSunEarthCR3BP(date, timeScale).getSecondary();
Assert.assertNotNull(secondaryBody);
}
@Test
public void testGetMu() {
final double mu = CR3BPFactory.getSunJupiterCR3BP().getMassRatio();
AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
TimeScale timeScale = TimeScalesFactory.getUTC();
final double mu = CR3BPFactory.getSunJupiterCR3BP(date, timeScale).getMassRatio();
Assert.assertNotNull(mu);
}
@Test
public void testGetName() {
final String name = CR3BPFactory.getSunEarthCR3BP().getName();
AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
TimeScale timeScale = TimeScalesFactory.getUTC();
final String name = CR3BPFactory.getSunEarthCR3BP(date, timeScale).getName();
Assert.assertNotNull(name);
}
......@@ -113,7 +126,9 @@ public class CR3BPSystemTest {
@Test
public void testGetGamma() {
final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP();
TimeScale timeScale = TimeScalesFactory.getUTC();
AbsoluteDate date = new AbsoluteDate(2020, 7, 5, 12, 0, 0.0, timeScale);
final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP(date, timeScale);
final double l1Gamma = syst.getGamma(LagrangianPoints.L1);
Assert.assertEquals(1.497655E9, l1Gamma * syst.getDdim(),1E3);
......@@ -131,11 +146,12 @@ public class CR3BPSystemTest {
public void testGetRealAPV() {
// Time settings
final TimeScale timeScale = TimeScalesFactory.getUTC();
final AbsoluteDate initialDate =
new AbsoluteDate(1996, 06, 26, 0, 0, 00.000,
TimeScalesFactory.getUTC());
final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP();
final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP(initialDate, timeScale);
final CelestialBody primaryBody = syst.getPrimary();
final CelestialBody secondaryBody = syst.getSecondary();
......
......@@ -32,6 +32,7 @@ import org.orekit.bodies.CelestialBody;
import org.orekit.bodies.CelestialBodyFactory;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.Constants;
import org.orekit.utils.FieldPVCoordinates;
......@@ -112,15 +113,16 @@ public class CR3BPRotatingTransformProviderTest {
final CelestialBody sun = CelestialBodyFactory.getSun();
final CelestialBody earth = CelestialBodyFactory.getEarth();
// Time settings
final TimeScale timeScale = TimeScalesFactory.getUTC();
final AbsoluteDate date = new AbsoluteDate(2000, 01, 01, 0, 0, 00.000,
timeScale);
// Set frames
final Frame sunFrame = sun.getInertiallyOrientedFrame();
final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP();
final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP(date, timeScale);
final Frame baryFrame = syst.getRotatingFrame();
// Time settings
final AbsoluteDate date = new AbsoluteDate(2000, 01, 01, 0, 0, 00.000,
TimeScalesFactory.getUTC());
// Compute Earth position in Sun centered frame
PVCoordinates pvEarth = earth.getPVCoordinates(date, sunFrame);
Vector3D posEarth = pvEarth.getPosition();
......@@ -146,15 +148,16 @@ public class CR3BPRotatingTransformProviderTest {
final CelestialBody sun = CelestialBodyFactory.getSun();
final CelestialBody earth = CelestialBodyFactory.getEarth();
// Time settings
final TimeScale timeScale = TimeScalesFactory.getUTC();
final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2000, 01, 01, 0, 0, 00.000,
timeScale);
// Set frames
final Frame sunFrame = sun.getInertiallyOrientedFrame();
final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP();
final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP(date.toAbsoluteDate(), timeScale);
final Frame baryFrame = syst.getRotatingFrame();
// Time settings
final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2000, 01, 01, 0, 0, 00.000,
TimeScalesFactory.getUTC());
// Compute Earth position in Sun centered frame
FieldPVCoordinates<T> pvEarth = earth.getPVCoordinates(date, sunFrame);
FieldVector3D<T> posEarth = pvEarth.getPosition();
......@@ -176,15 +179,16 @@ public class CR3BPRotatingTransformProviderTest {
final CelestialBody sun = CelestialBodyFactory.getSun();
final CelestialBody jupiter = CelestialBodyFactory.getJupiter();
// Time settings
final TimeScale timeScale = TimeScalesFactory.getUTC();
final AbsoluteDate date = new AbsoluteDate(2000, 01, 01, 0, 0, 00.000,
timeScale);
// Set frames
final Frame sunFrame = sun.getInertiallyOrientedFrame();
final CR3BPSystem syst = CR3BPFactory.getSunJupiterCR3BP();
final CR3BPSystem syst = CR3BPFactory.getSunJupiterCR3BP(date, timeScale);
final Frame baryFrame = syst.getRotatingFrame();
// Time settings
final AbsoluteDate date = new AbsoluteDate(2000, 01, 01, 0, 0, 00.000,
TimeScalesFactory.getUTC());
// Compute Jupiter position in Sun centered frame
PVCoordinates pvJupiter = jupiter.getPVCoordinates(date, sunFrame);
Vector3D posJupiter = pvJupiter.getPosition();
......@@ -210,15 +214,16 @@ public class CR3BPRotatingTransformProviderTest {
final CelestialBody sun = CelestialBodyFactory.getSun();
final CelestialBody jupiter = CelestialBodyFactory.getJupiter();
// Time settings
final TimeScale timeScale = TimeScalesFactory.getUTC();
final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2000, 01, 01, 0, 0, 00.000,
timeScale);
// Set frames
final Frame sunFrame = sun.getInertiallyOrientedFrame();
final CR3BPSystem syst = CR3BPFactory.getSunJupiterCR3BP();
final CR3BPSystem syst = CR3BPFactory.getSunJupiterCR3BP(date.toAbsoluteDate(), timeScale);
final Frame baryFrame = syst.getRotatingFrame();
// Time settings
final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2000, 01, 01, 0, 0, 00.000,
TimeScalesFactory.getUTC());
// Compute Jupiter position in Sun centered frame
FieldPVCoordinates<T> pvJupiter = jupiter.getPVCoordinates(date, sunFrame);
FieldVector3D<T> posJupiter = pvJupiter.getPosition();
......@@ -236,11 +241,12 @@ public class CR3BPRotatingTransformProviderTest {
@Test
public void testBaryOrientation() {
final TimeScale timeScale = TimeScalesFactory.getUTC();
final AbsoluteDate date0 = new AbsoluteDate(2000, 01, 1, 11, 58, 20.000,
TimeScalesFactory.getUTC());
timeScale);
final CelestialBody sun = CelestialBodyFactory.getSun();
final CelestialBody earth = CelestialBodyFactory.getEarth();
final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP();
final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP(date0, timeScale);
final Frame baryFrame = syst.getRotatingFrame();
for (double dt = -Constants.JULIAN_DAY; dt <= Constants.JULIAN_DAY; dt += 3600.0) {
final AbsoluteDate date = date0.shiftedBy(dt);
......@@ -258,11 +264,12 @@ public class CR3BPRotatingTransformProviderTest {
private <T extends RealFieldElement<T>> void doTestFieldBaryOrientation(final Field<T> field) {
final TimeScale timeScale = TimeScalesFactory.getUTC();
final FieldAbsoluteDate<T> date0 = new FieldAbsoluteDate<>(field, 2000, 01, 1, 11, 58, 20.000,
TimeScalesFactory.getUTC());
timeScale);
final CelestialBody sun = CelestialBodyFactory.getSun();
final CelestialBody earth = CelestialBodyFactory.getEarth();
final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP();
final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP(date0.toAbsoluteDate(), timeScale);
final Frame baryFrame = syst.getRotatingFrame();
for (double dt = -Constants.JULIAN_DAY; dt <= Constants.JULIAN_DAY; dt += 3600.0) {
final FieldAbsoluteDate<T> date = date0.shiftedBy(dt);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment