Java tutorial
/* Copyright 2002-2015 CS Systmes d'Information * Licensed to CS Systmes d'Information (CS) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. * CS licenses this file to You under the Apache License, Version 2.0 * (the "License"); you may not use this file except in compliance with * the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ package org.orekit.attitudes; import java.util.ArrayList; import java.util.List; import org.apache.commons.math3.fitting.PolynomialCurveFitter; import org.apache.commons.math3.fitting.WeightedObservedPoint; import org.apache.commons.math3.geometry.euclidean.threed.Line; import org.apache.commons.math3.geometry.euclidean.threed.Rotation; import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; import org.apache.commons.math3.util.FastMath; import org.apache.commons.math3.util.MathArrays; import org.junit.After; import org.junit.Assert; import org.junit.Before; import org.junit.Test; import org.orekit.Utils; import org.orekit.bodies.GeodeticPoint; import org.orekit.bodies.OneAxisEllipsoid; import org.orekit.errors.OrekitException; import org.orekit.frames.FramesFactory; import org.orekit.frames.Transform; import org.orekit.orbits.CircularOrbit; import org.orekit.orbits.PositionAngle; import org.orekit.propagation.SpacecraftState; import org.orekit.propagation.analytical.EcksteinHechlerPropagator; import org.orekit.time.AbsoluteDate; import org.orekit.time.DateComponents; import org.orekit.time.TimeComponents; import org.orekit.time.TimeScalesFactory; import org.orekit.utils.AngularCoordinates; import org.orekit.utils.Constants; import org.orekit.utils.IERSConventions; import org.orekit.utils.PVCoordinates; import org.orekit.utils.TimeStampedPVCoordinates; public class BodyCenterPointingTest { // Computation date private AbsoluteDate date; // Orbit private CircularOrbit circ; // WGS84 Earth model private OneAxisEllipsoid earth; // Transform from EME2000 to ITRF2008 private Transform eme2000ToItrf; // Earth center pointing attitude provider private BodyCenterPointing earthCenterAttitudeLaw; /** Test if target is on Earth surface */ @Test public void testTarget() throws OrekitException { // Call get target method TimeStampedPVCoordinates target = earthCenterAttitudeLaw.getTargetPV(circ, date, circ.getFrame()); // Check that target is on Earth surface GeodeticPoint gp = earth.transform(target.getPosition(), circ.getFrame(), date); Assert.assertEquals(0.0, gp.getAltitude(), 1.0e-10); Assert.assertEquals(date, target.getDate()); } /** Test if body center belongs to the direction pointed by the satellite */ @Test public void testBodyCenterInPointingDirection() throws OrekitException { // Transform satellite position to position/velocity parameters in EME2000 frame PVCoordinates pvSatEME2000 = circ.getPVCoordinates(); // Pointing direction // ******************** // Get satellite attitude rotation, i.e rotation from EME2000 frame to satellite frame Rotation rotSatEME2000 = earthCenterAttitudeLaw.getAttitude(circ, date, circ.getFrame()).getRotation(); // Transform Z axis from satellite frame to EME2000 Vector3D zSatEME2000 = rotSatEME2000.applyInverseTo(Vector3D.PLUS_K); // Transform Z axis from EME2000 to ITRF2008 Vector3D zSatITRF2008C = eme2000ToItrf.transformVector(zSatEME2000); // Transform satellite position/velocity from EME2000 to ITRF2008 PVCoordinates pvSatITRF2008C = eme2000ToItrf.transformPVCoordinates(pvSatEME2000); // Line containing satellite point and following pointing direction Line pointingLine = new Line(pvSatITRF2008C.getPosition(), pvSatITRF2008C.getPosition().add(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, zSatITRF2008C), 2.0e-8); // Check that the line contains Earth center Assert.assertTrue(pointingLine.contains(Vector3D.ZERO)); } @Test public void testQDot() throws OrekitException { Utils.setDataRoot("regular-data"); final double ehMu = 3.9860047e14; final double ae = 6.378137e6; final double c20 = -1.08263e-3; final double c30 = 2.54e-6; final double c40 = 1.62e-6; final double c50 = 2.3e-7; final double c60 = -5.5e-7; final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.); final Vector3D position = new Vector3D(3220103., 69623., 6449822.); final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.); final CircularOrbit initialOrbit = new CircularOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), date, ehMu); EcksteinHechlerPropagator propagator = new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60); propagator.setAttitudeProvider(earthCenterAttitudeLaw); List<WeightedObservedPoint> w0 = new ArrayList<WeightedObservedPoint>(); List<WeightedObservedPoint> w1 = new ArrayList<WeightedObservedPoint>(); List<WeightedObservedPoint> w2 = new ArrayList<WeightedObservedPoint>(); List<WeightedObservedPoint> w3 = new ArrayList<WeightedObservedPoint>(); for (double dt = -1; dt < 1; dt += 0.01) { Rotation rP = propagator.propagate(date.shiftedBy(dt)).getAttitude().getRotation(); w0.add(new WeightedObservedPoint(1, dt, rP.getQ0())); w1.add(new WeightedObservedPoint(1, dt, rP.getQ1())); w2.add(new WeightedObservedPoint(1, dt, rP.getQ2())); w3.add(new WeightedObservedPoint(1, dt, rP.getQ3())); } double q0DotRef = PolynomialCurveFitter.create(2).fit(w0)[1]; double q1DotRef = PolynomialCurveFitter.create(2).fit(w1)[1]; double q2DotRef = PolynomialCurveFitter.create(2).fit(w2)[1]; double q3DotRef = PolynomialCurveFitter.create(2).fit(w3)[1]; Attitude a0 = propagator.propagate(date).getAttitude(); double q0 = a0.getRotation().getQ0(); double q1 = a0.getRotation().getQ1(); double q2 = a0.getRotation().getQ2(); double q3 = a0.getRotation().getQ3(); double oX = a0.getSpin().getX(); double oY = a0.getSpin().getY(); double oZ = a0.getSpin().getZ(); // first time-derivatives of the quaternion double q0Dot = 0.5 * MathArrays.linearCombination(-q1, oX, -q2, oY, -q3, oZ); double q1Dot = 0.5 * MathArrays.linearCombination(q0, oX, -q3, oY, q2, oZ); double q2Dot = 0.5 * MathArrays.linearCombination(q3, oX, q0, oY, -q1, oZ); double q3Dot = 0.5 * MathArrays.linearCombination(-q2, oX, q1, oY, q0, oZ); Assert.assertEquals(q0DotRef, q0Dot, 5.0e-9); Assert.assertEquals(q1DotRef, q1Dot, 5.0e-9); Assert.assertEquals(q2DotRef, q2Dot, 5.0e-9); Assert.assertEquals(q3DotRef, q3Dot, 5.0e-9); } @Test public void testSpin() throws OrekitException { Utils.setDataRoot("regular-data"); final double ehMu = 3.9860047e14; final double ae = 6.378137e6; final double c20 = -1.08263e-3; final double c30 = 2.54e-6; final double c40 = 1.62e-6; final double c50 = 2.3e-7; final double c60 = -5.5e-7; final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.); final Vector3D position = new Vector3D(3220103., 69623., 6449822.); final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.); final CircularOrbit initialOrbit = new CircularOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), date, ehMu); EcksteinHechlerPropagator propagator = new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60); propagator.setAttitudeProvider(earthCenterAttitudeLaw); double h = 0.01; SpacecraftState s0 = propagator.propagate(date); SpacecraftState sMinus = propagator.propagate(date.shiftedBy(-h)); SpacecraftState sPlus = propagator.propagate(date.shiftedBy(h)); // check spin is consistent with attitude evolution double errorAngleMinus = Rotation.distance(sMinus.shiftedBy(h).getAttitude().getRotation(), s0.getAttitude().getRotation()); double evolutionAngleMinus = Rotation.distance(sMinus.getAttitude().getRotation(), s0.getAttitude().getRotation()); Assert.assertEquals(0.0, errorAngleMinus, 1.0e-6 * evolutionAngleMinus); double errorAnglePlus = Rotation.distance(s0.getAttitude().getRotation(), sPlus.shiftedBy(-h).getAttitude().getRotation()); double evolutionAnglePlus = Rotation.distance(s0.getAttitude().getRotation(), sPlus.getAttitude().getRotation()); Assert.assertEquals(0.0, errorAnglePlus, 1.0e-6 * evolutionAnglePlus); Vector3D spin0 = s0.getAttitude().getSpin(); Vector3D reference = AngularCoordinates.estimateRate(sMinus.getAttitude().getRotation(), sPlus.getAttitude().getRotation(), 2 * h); Assert.assertTrue(spin0.getNorm() > 1.0e-3); Assert.assertEquals(0.0, spin0.subtract(reference).getNorm(), 1.0e-13); } @Before public void setUp() { try { Utils.setDataRoot("regular-data"); // Computation date date = new AbsoluteDate(new DateComponents(2008, 04, 07), TimeComponents.H00, TimeScalesFactory.getUTC()); // Satellite position as circular parameters final double mu = 3.9860047e14; final double raan = 270.; circ = new CircularOrbit(7178000.0, 0.5e-4, -0.5e-4, FastMath.toRadians(50.), FastMath.toRadians(raan), FastMath.toRadians(5.300 - raan), PositionAngle.MEAN, FramesFactory.getEME2000(), date, mu); // WGS84 Earth model earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true)); // Transform from EME2000 to ITRF2008 eme2000ToItrf = FramesFactory.getEME2000().getTransformTo(earth.getBodyFrame(), date); // Create earth center pointing attitude provider */ earthCenterAttitudeLaw = new BodyCenterPointing(circ.getFrame(), earth); } catch (OrekitException oe) { Assert.fail(oe.getMessage()); } } @After public void tearDown() { date = null; earth = null; eme2000ToItrf = null; earthCenterAttitudeLaw = null; circ = null; } }