We want to compute roughly the Doppler effect of a satellite with respect to a ground station. Note that for an accurate computation, including light time delay, tropospheric effects and others, the measurements generation feature would be better suited. This example is only meant as a basic tutorial.
On the one hand, let's first define some initial state for the satellite with:
The initial orbit is just set as a CartesianOrbit.
More details on the orbit representation can be found
in the orbits section
of the library architecture documentation.
Frame inertialFrame = FramesFactory.getEME2000();
TimeScale utc = TimeScalesFactory.getUTC();
AbsoluteDate initialDate = new AbsoluteDate(2008, 10, 01, 0, 0, 00.000, utc);
double mu = 3.986004415e+14;
Vector3D position = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
Vector3D velocity = new Vector3D(505.8479685, 942.7809215, 7435.922231);
PVCoordinates pvsat = new PVCoordinates(position, velocity);
Orbit initialOrbit = new CartesianOrbit(pvsat, inertialFrame, initialDate, mu);
As a propagator, we consider a simple KeplerianPropagator that will propagate orbit in the inertial
frame in which the initial orbit was defined (here EME2000).
Propagator kepler = new KeplerianPropagator(initialOrbit);
On the other hand, let's define the ground station by its coordinates as a GeodeticPoint
in its own TopocentricFrame related to a BodyShape in some terrestrial frame.
Frame earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING,
earthFrame);
double longitude = FastMath.toRadians(45.);
double latitude = FastMath.toRadians(25.);
double altitude = 0.;
GeodeticPoint station = new GeodeticPoint(latitude, longitude, altitude);
TopocentricFrame staF = new TopocentricFrame(earth, station, "station");
More details on BodyShape and GeodeticPoint can be found
in the bodies section
of the library architecture documentation.
More details on TopocentricFrame can be found
in the frames section
of the library architecture documentation.
Finally, we can get the Doppler measurement in a simple propagation loop
AbsoluteDate finalDate = initialDate.shiftedBy(6000);
AbsoluteDate extrapDate = initialDate;
while (extrapDate.compareTo(finalDate) <= 0) {
// We can simply get the position and velocity of spacecraft in station frame at any time
PVCoordinates pvInert = kepler.propagate(extrapDate).getPVCoordinates();
PVCoordinates pvStation = inertialFrame.getTransformTo(staF, extrapDate).transformPVCoordinates(pvInert);
// And then calculate the Doppler signal
double doppler = Vector3D.dotProduct(pvStation.getPosition(), pvStation.getVelocity()) / pvStation.getPosition().getNorm();
System.out.format(Locale.US, "%s %9.3f%n", extrapDate, doppler);
extrapDate = extrapDate.shiftedBy(600);
}
Here are the results displayed by this program:
time doppler (m/s)
2008-10-01T00:00:00.000 -2925.114
2008-10-01T00:10:00.000 -3069.084
2008-10-01T00:20:00.000 -1331.910
2008-10-01T00:30:00.000 1664.611
2008-10-01T00:40:00.000 3143.549
2008-10-01T00:50:00.000 2832.906
2008-10-01T01:00:00.000 1556.662
2008-10-01T01:10:00.000 -140.889
2008-10-01T01:20:00.000 -1860.637
2008-10-01T01:30:00.000 -3195.728
2008-10-01T01:40:00.000 -3538.053
The complete code for this example can be found in the source tree of the tutorials,
in file src/main/java/org/orekit/tutorials/frames/Frames1.java.
The problem is related to the one exposed in the User defined frames subsection of the frames section of the library architecture documentation.
It can be summarized by the following schema.
For a given satellite, GPS measurements for position and velocity, expressed in the ITRF frame, are available at any moment. The GPS antenna is fixed with some offset with respect to the satellite reference frame. The attitude of the satellite reference frame with respect to some moving frame related to the satellite center of gravity (CoG) is known at any moment. We want to compute for some instant the position and velocity of the CoG in the EME2000 inertial frame.
A possible solution provided by OREKIT is detailed above.
The CoG frame has its origin at the satellite center of gravity and is aligned with EME2000. It is linked to its parent EME2000 frame by an a priori unknown transform which depends on the current position and velocity. This transform is what we want to compute.
We first build the frame. We use the identity transform as a simple dummy value, the real value which is time-dependent will be recomputed when time is reset:
Frame cogFrame = new UpdatableFrame(FramesFactory.getEME2000(), Transform.IDENTITY, "LOF", false);
The satellite frame, with origin also at the CoG, depends on attitude. For the sake of this tutorial, we consider a simple inertial attitude here:
Transform cogToSat = new Transform(new Rotation(0.6, 0.48, 0, 0.64, false));
Frame satFrame = new Frame(cogFrame, cogToSat, "sat", false);
Finally, the GPS antenna frame is always defined from the satellite frame by 2 transforms: a translation and a rotation. Let's set some values:
Transform translateGPS = new Transform(new Vector3D(0, 0, 1));
Transform rotateGPS = new Transform(date,
new Rotation(new Vector3D(0, 1, 3),
FastMath.toRadians(10),
RotationConvention.VECTOR_OPERATOR));
Frame gpsFrame = new Frame(satFrame,
new Transform(date, translateGPS, rotateGPS),
"GPS", false);
Let's consider some measurement date in UTC time scale:
TimeScale utc = TimeScalesFactory.getUTC();
AbsoluteDate date = new AbsoluteDate(2008, 10, 01, 12, 00, 00.000, utc);
Then let's get some satellite position and velocity in ITRF as measured by GPS antenna at this moment:
Vector3D position = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
Vector3D velocity = new Vector3D(505.8479685, 942.7809215, 7435.922231);
The transform from GPS frame to ITRF frame at this moment is defined by a translation and a rotation. The translation is directly provided by the GPS measurement above. The rotation is extracted from the existing tree, where we know all rotations are already up to date, even if one translation is still unknown. We combine the extracted rotation and the measured translation by applying the rotation first because the position/velocity vector are given in ITRF frame not in GPS antenna frame:
Transform measuredTranslation = new Transform(date, position, velocity);
Transform formerTransform = gpsFrame.getTransformTo(FramesFactory.getITRF(IERSConventions.IERS_2010, true), date);
Transform preservedRotation = new Transform(date,
formerTransform.getRotation(),
formerTransform.getRotationRate());
Transform gpsToItrf = new Transform(date, preservedRotation, measuredTranslation);
So we can now update the transform from EME2000 to CoG frame.
The updateTransform method automatically locates the frames
in the global tree, combines all transforms and updates the single
transform between CoG itself and its parent EME2000 frame.
cogFrame.updateTransform(gpsFrame, FramesFactory.getITRF(IERSConventions.IERS_2010, true), gpsToItrf, date);
The frame tree is now in sync. We can compute all position and velocity pairs we want:
PVCoordinates origin = PVCoordinates.ZERO;
Transform cogToItrf = cogFrame.getTransformTo(FramesFactory.getITRF(IERSConventions.IERS_2010, true), date);
PVCoordinates satItrf = cogToItrf.transformPVCoordinates(origin);
Transform cogToEme2000 = cogFrame.getTransformTo(FramesFactory.getEME2000(), date);
PVCoordinates satEME2000 = cogToEme2000.transformPVCoordinates(origin);
As a result, we can compare the GPS measurements to the computed values:
GPS antenna position in ITRF: -6142438.668 3492467.560 -25767.257
GPS antenna velocity in ITRF: 505.8479685 942.7809215 7435.9222310
Satellite position in ITRF: -6142439.167 3492468.238 -25766.717
Satellite velocity in ITRF: 505.8480179 942.7809579 7435.9222310
Satellite position in EME2000: 6675017.354 -2317478.630 -31517.561
Satellite velocity in EME2000: -150.5212935 -532.0401727 7436.0739031
The complete code for this example can be found in the source tree of the tutorials,
in file src/main/java/org/orekit/tutorials/frames/Frames2.java.
Let's consider a spacecraft which attitude is guided by a yaw steering law, we want to plot some of the effects of such a law on the spacecraft motion.
First, an initial orbit for the spacecraft is defined as follows:
Frame eme2000 = FramesFactory.getEME2000();
AbsoluteDate initialDate = new AbsoluteDate(2003, 4, 7, 10, 55, 21.575,
TimeScalesFactory.getUTC());
double mu = 3.986004415e+14;
Orbit orbit = new CircularOrbit(7178000.0, 0.5e-4, -0.5e-4,
FastMath.toRadians(50.),
FastMath.toRadians(220.),
FastMath.toRadians(5.300),
PositionAngleType.MEAN,
eme2000, initialDate, mu);
More details on the orbit representation can be found in the orbits section of the library architecture documentation.
Then the attitude law for the spacecraft is constructed:
The yaw steering law wraps the nadir point law in order to give maximal lighting to the solar arrays, which rotation axis is Y in spacecraft frame:
Frame earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING,
earthFrame);
NadirPointing nadirLaw = new NadirPointing(eme2000, earth);
final PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
YawSteering yawSteeringLaw = new YawSteering(eme2000, nadirLaw, sun, Vector3D.MINUS_I);
More details on the attitude representation can be found in the attitudes section of the library architecture documentation.
As a propagator, we consider an analytic EcksteinHechlerPropagator.
Propagator propagator =
new EcksteinHechlerPropagator(orbit, yawSteeringLaw,
Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
Constants.EIGEN5C_EARTH_MU,
Constants.EIGEN5C_EARTH_C20,
Constants.EIGEN5C_EARTH_C30,
Constants.EIGEN5C_EARTH_C40,
Constants.EIGEN5C_EARTH_C50,
Constants.EIGEN5C_EARTH_C60);
More details on Propagator can be found in the
propagation section of the library
architecture documentation.
Then, we associate a step handler with the propagator which directly computes the desired data in spacecraft frame using the current spacecraft state
propagator.getMultiplexer().add(10, new OrekitFixedStepHandler() {
public void init(SpacecraftState s0, AbsoluteDate t, final double step)
// write file header
}
public void handleStep(SpacecraftState currentState) {
Transform inertToSpacecraft = currentState.toTransform();
Vector3D sunInert = sun.getPVCoordinates(currentState.getDate(), currentState.getFrame()).getPosition();
Vector3D sunSat = inertToSpacecraft.transformPosition(sunInert);
Vector3D spin = inertToSpacecraft.getRotationRate();
// write data to file
}
public void finish(final SpacecraftState finalState) {
// close the file
}
}
So, computing over an orbit, we can plot the following results showing clearly a typical yaw steering behavior:
The complete code for this example can be found in the source tree of the tutorials,
in file src/main/java/org/orekit/tutorials/frames/Frames3.java.