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.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.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.Frame; import org.orekit.frames.FramesFactory; import org.orekit.frames.Transform; import org.orekit.orbits.CircularOrbit; import org.orekit.orbits.KeplerianOrbit; import org.orekit.orbits.Orbit; import org.orekit.orbits.PositionAngle; import org.orekit.propagation.Propagator; import org.orekit.propagation.SpacecraftState; import org.orekit.propagation.analytical.KeplerianPropagator; 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.CartesianDerivativesFilter; import org.orekit.utils.Constants; import org.orekit.utils.IERSConventions; import org.orekit.utils.PVCoordinates; import org.orekit.utils.TimeStampedPVCoordinates; public class YawCompensationTest { // Computation date private AbsoluteDate date; // Reference frame = ITRF 2005C private Frame itrf; // Satellite position CircularOrbit circOrbit; // Earth shape OneAxisEllipsoid earthShape; /** Test that pointed target remains the same with or without yaw compensation */ @Test public void testTarget() throws OrekitException { // Attitude laws // ************** // Target pointing attitude provider without yaw compensation NadirPointing nadirLaw = new NadirPointing(circOrbit.getFrame(), earthShape); // Target pointing attitude provider with yaw compensation YawCompensation yawCompensLaw = new YawCompensation(circOrbit.getFrame(), nadirLaw); // Check target // ************* // without yaw compensation TimeStampedPVCoordinates noYawObserved = nadirLaw.getTargetPV(circOrbit, date, itrf); // with yaw compensation TimeStampedPVCoordinates yawObserved = yawCompensLaw.getTargetPV(circOrbit, date, itrf); // Check difference PVCoordinates observedDiff = new PVCoordinates(yawObserved, noYawObserved); Assert.assertEquals(0.0, observedDiff.getPosition().getNorm(), Utils.epsilonTest); Assert.assertEquals(0.0, observedDiff.getVelocity().getNorm(), Utils.epsilonTest); Assert.assertEquals(0.0, observedDiff.getAcceleration().getNorm(), Utils.epsilonTest); } /** Test the derivatives of the sliding target */ @Test public void testSlidingDerivatives() throws OrekitException { GroundPointing law = new YawCompensation(circOrbit.getFrame(), new NadirPointing(circOrbit.getFrame(), earthShape)); List<TimeStampedPVCoordinates> sample = new ArrayList<TimeStampedPVCoordinates>(); for (double dt = -0.1; dt < 0.1; dt += 0.01) { Orbit o = circOrbit.shiftedBy(dt); sample.add(law.getTargetPV(o, o.getDate(), o.getFrame())); } TimeStampedPVCoordinates reference = TimeStampedPVCoordinates.interpolate(circOrbit.getDate(), CartesianDerivativesFilter.USE_P, sample); TimeStampedPVCoordinates target = law.getTargetPV(circOrbit, circOrbit.getDate(), circOrbit.getFrame()); Assert.assertEquals(0.0, Vector3D.distance(reference.getPosition(), target.getPosition()), 1.0e-15 * reference.getPosition().getNorm()); Assert.assertEquals(0.0, Vector3D.distance(reference.getVelocity(), target.getVelocity()), 3.0e-11 * reference.getVelocity().getNorm()); Assert.assertEquals(0.0, Vector3D.distance(reference.getAcceleration(), target.getAcceleration()), 7.0e-6 * reference.getAcceleration().getNorm()); } /** Test that pointed target motion is along -X sat axis */ @Test public void testAlignment() throws OrekitException { GroundPointing notCompensated = new NadirPointing(circOrbit.getFrame(), earthShape); YawCompensation compensated = new YawCompensation(circOrbit.getFrame(), notCompensated); Attitude att0 = compensated.getAttitude(circOrbit, circOrbit.getDate(), circOrbit.getFrame()); // ground point in satellite Z direction Vector3D satInert = circOrbit.getPVCoordinates().getPosition(); Vector3D zInert = att0.getRotation().applyInverseTo(Vector3D.PLUS_K); GeodeticPoint gp = earthShape.getIntersectionPoint( new Line(satInert, satInert.add(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, zInert), 1.0e-10), satInert, circOrbit.getFrame(), circOrbit.getDate()); PVCoordinates pEarth = new PVCoordinates(earthShape.transform(gp), Vector3D.ZERO, Vector3D.ZERO); double minYWithoutCompensation = Double.POSITIVE_INFINITY; double maxYWithoutCompensation = Double.NEGATIVE_INFINITY; double minYDotWithoutCompensation = Double.POSITIVE_INFINITY; double maxYDotWithoutCompensation = Double.NEGATIVE_INFINITY; double minYWithCompensation = Double.POSITIVE_INFINITY; double maxYWithCompensation = Double.NEGATIVE_INFINITY; double minYDotWithCompensation = Double.POSITIVE_INFINITY; double maxYDotWithCompensation = Double.NEGATIVE_INFINITY; for (double dt = -0.2; dt < 0.2; dt += 0.002) { PVCoordinates withoutCompensation = toSpacecraft(pEarth, circOrbit.shiftedBy(dt), notCompensated); if (FastMath.abs(withoutCompensation.getPosition().getX()) <= 1000.0) { minYWithoutCompensation = FastMath.min(minYWithoutCompensation, withoutCompensation.getPosition().getY()); maxYWithoutCompensation = FastMath.max(maxYWithoutCompensation, withoutCompensation.getPosition().getY()); minYDotWithoutCompensation = FastMath.min(minYDotWithoutCompensation, withoutCompensation.getVelocity().getY()); maxYDotWithoutCompensation = FastMath.max(maxYDotWithoutCompensation, withoutCompensation.getVelocity().getY()); } PVCoordinates withCompensation = toSpacecraft(pEarth, circOrbit.shiftedBy(dt), compensated); if (FastMath.abs(withCompensation.getPosition().getX()) <= 1000.0) { minYWithCompensation = FastMath.min(minYWithCompensation, withCompensation.getPosition().getY()); maxYWithCompensation = FastMath.max(maxYWithCompensation, withCompensation.getPosition().getY()); minYDotWithCompensation = FastMath.min(minYDotWithCompensation, withCompensation.getVelocity().getY()); maxYDotWithCompensation = FastMath.max(maxYDotWithCompensation, withCompensation.getVelocity().getY()); } } // when the ground point is close to cross the push-broom line (i.e. when x decreases from +1000m to -1000m) // it will drift along the Y axis if we don't apply compensation // but will remain nearly at y=0 if we do apply compensation // in fact, as the yaw compensation mode removes the linear drift, // what remains is a parabola y = a ux Assert.assertEquals(-55.7056, minYWithoutCompensation, 0.0001); Assert.assertEquals(+55.7056, maxYWithoutCompensation, 0.0001); Assert.assertEquals(352.5667, minYDotWithoutCompensation, 0.0001); Assert.assertEquals(352.5677, maxYDotWithoutCompensation, 0.0001); Assert.assertEquals(0.0000, minYWithCompensation, 0.0001); Assert.assertEquals(0.0008, maxYWithCompensation, 0.0001); Assert.assertEquals(-0.0101, minYDotWithCompensation, 0.0001); Assert.assertEquals(0.0102, maxYDotWithCompensation, 0.0001); } PVCoordinates toSpacecraft(PVCoordinates groundPoint, Orbit orbit, AttitudeProvider attitudeProvider) throws OrekitException { SpacecraftState state = new SpacecraftState(orbit, attitudeProvider.getAttitude(orbit, orbit.getDate(), orbit.getFrame())); Transform earthToSc = new Transform(orbit.getDate(), earthShape.getBodyFrame().getTransformTo(orbit.getFrame(), orbit.getDate()), state.toTransform()); return earthToSc.transformPVCoordinates(groundPoint); } /** Test that maximum yaw compensation is at ascending/descending node, * and minimum yaw compensation is at maximum latitude. */ @Test public void testCompensMinMax() throws OrekitException { // Attitude laws // ************** // Target pointing attitude provider over satellite nadir at date, without yaw compensation NadirPointing nadirLaw = new NadirPointing(circOrbit.getFrame(), earthShape); // Target pointing attitude provider with yaw compensation YawCompensation yawCompensLaw = new YawCompensation(circOrbit.getFrame(), nadirLaw); // Extrapolation over one orbital period (sec) double duration = circOrbit.getKeplerianPeriod(); KeplerianPropagator extrapolator = new KeplerianPropagator(circOrbit); // Extrapolation initializations double delta_t = 15.0; // extrapolation duration in seconds AbsoluteDate extrapDate = date; // extrapolation start date // Min initialization double yawMin = 1.e+12; double latMin = 0.; while (extrapDate.durationFrom(date) < duration) { extrapDate = extrapDate.shiftedBy(delta_t); // Extrapolated orbit state at date Orbit extrapOrbit = extrapolator.propagate(extrapDate).getOrbit(); PVCoordinates extrapPvSatEME2000 = extrapOrbit.getPVCoordinates(); // Satellite latitude at date double extrapLat = earthShape .transform(extrapPvSatEME2000.getPosition(), FramesFactory.getEME2000(), extrapDate) .getLatitude(); // Compute yaw compensation angle -- rotations composition double yawAngle = yawCompensLaw.getYawAngle(extrapOrbit, extrapDate, extrapOrbit.getFrame()); // Update minimum yaw compensation angle if (FastMath.abs(yawAngle) <= yawMin) { yawMin = FastMath.abs(yawAngle); latMin = extrapLat; } // Checks // ------------------ // 1/ Check yaw values around ascending node (max yaw) if ((FastMath.abs(extrapLat) < FastMath.toRadians(2.)) && (extrapPvSatEME2000.getVelocity().getZ() >= 0.)) { Assert.assertEquals(-3.206, FastMath.toDegrees(yawAngle), 0.003); } // 2/ Check yaw values around maximum positive latitude (min yaw) if (extrapLat > FastMath.toRadians(50.15)) { Assert.assertEquals(0, FastMath.toDegrees(yawAngle), 0.15); } // 3/ Check yaw values around descending node (max yaw) if ((FastMath.abs(extrapLat) < FastMath.toRadians(2.)) && (extrapPvSatEME2000.getVelocity().getZ() <= 0.)) { Assert.assertEquals(3.206, FastMath.toDegrees(yawAngle), 0.003); } // 4/ Check yaw values around maximum negative latitude (min yaw) if (extrapLat < FastMath.toRadians(-50.15)) { Assert.assertEquals(0, FastMath.toDegrees(yawAngle), 0.15); } } // 5/ Check that minimum yaw compensation value is around maximum latitude Assert.assertEquals(0.0, FastMath.toDegrees(yawMin), 0.004); Assert.assertEquals(50.0, FastMath.toDegrees(latMin), 0.22); } /** Test that compensation rotation axis is Zsat, yaw axis */ @Test public void testCompensAxis() throws OrekitException { // Attitude laws // ************** // Target pointing attitude provider over satellite nadir at date, without yaw compensation NadirPointing nadirLaw = new NadirPointing(circOrbit.getFrame(), earthShape); // Target pointing attitude provider with yaw compensation YawCompensation yawCompensLaw = new YawCompensation(circOrbit.getFrame(), nadirLaw); // Get attitude rotations from non yaw compensated / yaw compensated laws Rotation rotNoYaw = nadirLaw.getAttitude(circOrbit, date, circOrbit.getFrame()).getRotation(); Rotation rotYaw = yawCompensLaw.getAttitude(circOrbit, date, circOrbit.getFrame()).getRotation(); // Compose rotations composition Rotation compoRot = rotYaw.applyTo(rotNoYaw.revert()); Vector3D yawAxis = compoRot.getAxis(); // Check axis Assert.assertEquals(0., yawAxis.subtract(Vector3D.PLUS_K).getNorm(), Utils.epsilonTest); } @Test public void testSpin() throws OrekitException { NadirPointing nadirLaw = new NadirPointing(circOrbit.getFrame(), earthShape); // Target pointing attitude provider with yaw compensation YawCompensation law = new YawCompensation(circOrbit.getFrame(), nadirLaw); KeplerianOrbit orbit = new KeplerianOrbit(7178000.0, 1.e-4, FastMath.toRadians(50.), FastMath.toRadians(10.), FastMath.toRadians(20.), FastMath.toRadians(30.), PositionAngle.MEAN, FramesFactory.getEME2000(), date.shiftedBy(-300.0), 3.986004415e14); Propagator propagator = new KeplerianPropagator(orbit, law); double h = 0.01; SpacecraftState sMinus = propagator.propagate(date.shiftedBy(-h)); SpacecraftState s0 = propagator.propagate(date); 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, 7.5e-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, 2.0e-5 * 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(), 2.0e-8); } @Before public void setUp() { try { Utils.setDataRoot("regular-data"); // Computation date date = new AbsoluteDate(new DateComponents(2008, 04, 07), TimeComponents.H00, TimeScalesFactory.getUTC()); // Body mu final double mu = 3.9860047e14; // Reference frame = ITRF 2005 itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true); // Satellite position circOrbit = new CircularOrbit(7178000.0, 0.5e-4, -0.5e-4, FastMath.toRadians(50.), FastMath.toRadians(270.), FastMath.toRadians(5.300), PositionAngle.MEAN, FramesFactory.getEME2000(), date, mu); // Elliptic earth shape */ earthShape = new OneAxisEllipsoid(6378136.460, 1 / 298.257222101, itrf); } catch (OrekitException oe) { Assert.fail(oe.getMessage()); } } @After public void tearDown() { date = null; itrf = null; circOrbit = null; earthShape = null; } }