# Doppler shift using Orekit

Hello,

I am a very new Orekit-user and currently working on the Doppler shift function in my project. I am attempting to use the basic Doppler shift formula for calculations (get inspiration from Calculating Doppler Shift for satellites - Amateur Radio Stack Exchange). However, I am unsure if I am doing it correctly, as I currently can not test it with physical radio equipment. I would appreciate it if someone could help me verify if I have implemented it correctly.

``````    public static double calcDopplerFactor(TLE tle, Position gsPos, AbsoluteDate date) {
TopocentricFrame groundStationFrame = new TopocentricFrame(
Constants.WGS84_EARTH_FLATTENING,
FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
"Ground Station");

// Initialize TLE propagator
TLEPropagator propagator = TLEPropagator.selectExtrapolator(tle);

// Propagate satellite position to current date
PVCoordinates pv = propagator.propagate(date).getPVCoordinates();

// Calculate satellite velocity
Vector3D satelliteVelocity = pv.getVelocity();

// Calculate ground station velocity
Vector3D groundStationVelocity = groundStationFrame.getPVCoordinates(date, FramesFactory.getITRF(IERSConventions.IERS_2010, true)).getVelocity();

// Calculate relative velocity
Vector3D relativeVelocity = satelliteVelocity.subtract(groundStationVelocity);

double relativeVelocityNorm = relativeVelocity.getNorm();

// Calculate Doppler shift
return relativeVelocityNorm / Constants.SPEED_OF_LIGHT;
}
``````

Hi @SUIFENGSK welcome

What you are computing here is not the Doppler shift but the norm of relative velocity. The Doppler shift is projected along the line of sight, i.e. you have to compute the dot product of relative velocity with the unit vector given by normalizing the relative position.

You should also take into account the light time delay, i.e. the relative position and velocity should not be computed using the same date, but rather the position/velocity of receiver at time of signal arrival and
the position/velocity of transmitter at time of signal departure.

Also note that there is a `RangeRate` class that does all of this (and more). It was designed for orbit determination, so it takes a range rate as one of its construction parameter, which looks like a chicken and egg problem. In fact, you can use `Double.NaN `, or 0, or any number when building a first dummy `RangeRate` instance, and then call its `estimateWithoutDerivatives` method to evaluate properly the Doppler effect. This is waht is done when using the measurements generation feature.

1 Like

I am currently attempting to use the RangeRate class as you mentioned. However, Iâ€™m unsure about some parameters in the constructor.

In the documentation, it mentions:

public RangeRateâ€‹(GroundStation station, AbsoluteDate date, double rangeRate, double sigma, double baseWeight, boolean twoway, ObservableSatellite satellite)

Parameters:
`station` - ground station from which measurement is performed
`date` - date of the measurement
`rangeRate` - observed value, m/s
`sigma` - theoretical standard deviation
`baseWeight` - base weight
`twoway` - if true, this is a two-way measurement
`satellite` - satellite related to this measurement

For the `rangeRate` parameter, you mentioned I can use `Double.NaN`, 0, or any number. So, what are the `sigma` and `baseWeight` values? Which values do I need to set?

Also, how can I build an `ObservableSatellite`? I see that the constructor of this class only takes a propagator index. Currently, I only have TLE data of the satellite.

As Romain would say, our measurements class hierarchy is deeply related to the orbit determination use, which explains all these parameters. They were intended to be used by the orbit determination engine (either batch least squares or Kalman filter) and are not really useful when you are building the measurements. So here again, you can just use any values, as these parameters wonâ€™t be used by the `estimate` methods.
As for the `ObservableSatellite` parameter, just build one with an index set to 0, and when calling the estimate method juste pass a states array with one element at index 0 to represent your satellite state.

1 Like

Iâ€™ve attempted to implement it as follows:

``````    public static double calcDopplerFactor(TLE tle, Position gsPos, AbsoluteDate date) {
GroundStation groundStation = new GroundStation(new TopocentricFrame(new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
TLEPropagator tlePropagator = TLEPropagator.selectExtrapolator(tle);
SpacecraftState spacecraftState = tlePropagator.propagate(date);
SpacecraftState[] states = new SpacecraftState[1];
states[0] = spacecraftState;
RangeRate rangeRate = new RangeRate(groundStation, date, 0, 0, 0, true, new ObservableSatellite(0));
double estimatedValue = rangeRate.estimateWithoutDerivatives(1, 1, states).getEstimatedValue()[0];
System.out.println("DopplerCorrectionHelper: calcDopplerFactor: estimatedValue: " + estimatedValue);
return estimatedValue;
}
``````

Firstly, the program gets stuck after calling `estimateWithoutDerivatives` . Iâ€™m not sure if Iâ€™m calling it correctly here. Secondly, `estimateWithoutDerivatives` has three parameters; the first two are for `iteration` and `evaluation`. Which values should I set for them? Currently, Iâ€™ve set them both to 1.

The geodetic point should refer to that latitude, longitude and altitude of the ground station.

I have now first converted my position to latitude and longitude and then set it to a `GeodeticPoint`. However, it seems that the program still gets stuck at the same line:

`double estimatedValue = rangeRate.estimateWithoutDerivatives(1, 1, states).getEstimatedValue()[0];`

I have placed a breakpoint before this, so could you please check if there is still something wrong?

I donâ€™t see anything obvious (apart from the using â€śelevationâ€ť and â€śazimuthâ€ť to name latitude and longitude which is confusing to me).
Are you sure that at the date at which you compute the range rate the satellite is visible from the ground station (i.e. it is at least above horizon)?

1 Like

Yes, I did not take the visibility of the satellite into account. However, now Iâ€™m trying to specify a point in time when the satellite is over the horizon. Despite this, the program still gets stuck at that line.

I live in Denmark and Iâ€™m using the following TLE data and time for the ISS (ZARYA):

``````1 25544U 98067A   24103.85637196  .00013837  00000-0  24880-3 0  9992
2 25544  51.6381 279.3375 0004725  59.1380 301.0073 15.50131378448415
``````

Time:

``````date = new AbsoluteDate(2024, 4, 12, 8, 28, 0, TimeScalesFactory.getUTC());
``````

Ground station position:

`Latitude: 55.7704, Longitude: 12.5037, Altitude: 0`

I also found a website(https://www.satflare.com/track.asp?q=25544#TOP) that shows where the ISS is for a given time:

I have overwritten all parameters in the function and hardcoded everything.

``````    public static double calcDopplerFactor(TLE tle, Position gsPos, AbsoluteDate date) {
tle = new TLE("1 25544U 98067A   24103.85637196  .00013837  00000-0  24880-3 0  9992",
"2 25544  51.6381 279.3375 0004725  59.1380 301.0073 15.50131378448415");
date = new AbsoluteDate(2024, 4, 12, 8, 28, 0, TimeScalesFactory.getUTC());
// Location of the Ground Station (GS)
double latitude = 55.7704; // GS latitude in degrees
double longitude = 12.5037; // GS longitude in degrees
double altitude = 0; // Altitude in meters (here assumed to be sea level)
GroundStation groundStation = new GroundStation(new TopocentricFrame(new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
TLEPropagator tlePropagator = TLEPropagator.selectExtrapolator(tle);
SpacecraftState spacecraftState = tlePropagator.propagate(date);
SpacecraftState[] states = new SpacecraftState[1];
states[0] = spacecraftState;
RangeRate rangeRate = new RangeRate(groundStation, date, 0, 0, 0, true, new ObservableSatellite(0));
double estimatedValue = rangeRate.estimateWithoutDerivatives(1, 1, states).getEstimatedValue()[0];
System.out.println("DopplerCorrectionHelper: calcDopplerFactor: estimatedValue: " + estimatedValue);
return estimatedValue;
}
``````

Could you add this after you built the ground station:

``````        for (ParameterDriver driver : Arrays.asList(groundStation.getClockOffsetDriver(),
groundStation.getEastOffsetDriver(),
groundStation.getNorthOffsetDriver(),
groundStation.getZenithOffsetDriver(),
groundStation.getPrimeMeridianOffsetDriver(),
groundStation.getPrimeMeridianDriftDriver(),
groundStation.getPolarOffsetXDriver(),
groundStation.getPolarDriftXDriver(),
groundStation.getPolarOffsetYDriver(),
groundStation.getPolarDriftYDriver())) {
if (driver.getReferenceDate() == null) {
driver.setReferenceDate(date);
}
}
``````

If I do so, I get this result:

``````DopplerCorrectionHelper: calcDopplerFactor: estimatedValue: -3577.1542958506725
``````
1 Like

It works now! Thank you!

I have some additional questions to help me understand this more deeply:

1. The value we are trying to estimate here is Doppler shift, not Doppler factor, right? So, can I simply add this value to my current radio frequency?
2. Could you explain the first two parameters in the `estimateWithoutDerivatives` function - `iteration` and `evaluation`? Which values should I set for them, and what factors do I need to take into account?
3. What is the difference between using `rangeRate.estimateWithoutDerivatives` and `rangeRate.estimate` to estimate Doppler shift?"
1 Like

The computed value is a range rate in meters per seconds.
Iteration and evaluation are not used at all by `estimateWithoutDerivatives`, they are just past along to the `EstimatedMeasurement` object so they can be retrieved by further processing. As you only extract the measurement value, you can ignore them. They are used by the orbit determination process to count iterations, as some modifiers that can be applied later on to the estimated measurement may depend on these counts (for example the outlier filters only act after a warmup and needs to look at the count for that).
The difference is that the second method also computes partial derivatives, for example with respect to the state or with respect to the parameters drivers I mentioned earlier. This implies additional information is computed using algorithmic differentiation and that the parameter drivers have been set up properly. You donâ€™t need this I think.

To calculate the Doppler shift, I need to multiply it by f_0/c?

I use the following equation:

Yes, because what is computed is vrange.

2 Likes

Hi @luc

I recently applied the Doppler shift to my physical radio system in the real world, but I noticed that the sign of the Doppler factor is reversed. When the satellite is moving toward the ground station, the Doppler factor should have a positive sign and result in a positive Doppler effect for the RX radio. However, my code above produces a negative sign. The value is correct, but the sign is reversed. Do you have any idea why this might be happening?

Currently, the code is as follows:

``````    public double calcDopplerFactor(TLE tle, Position gsPos, Instant obsTime) {
// Location of the Ground Station (GS)
double latitude = gsPos.elevation(); // GS latitude in degrees
double longitude = gsPos.azimuth(); // GS longitude in degrees
double altitude = gsPos.distance(); // Altitude in meters
GroundStation groundStation = new GroundStation(new TopocentricFrame(new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true)),

AbsoluteDate observerTime = new AbsoluteDate(obsTime, TimeScalesFactory.getUTC());

// Build the ground station parameters
for (ParameterDriver driver : Arrays.asList(groundStation.getClockOffsetDriver(),
groundStation.getEastOffsetDriver(),
groundStation.getNorthOffsetDriver(),
groundStation.getZenithOffsetDriver(),
groundStation.getPrimeMeridianOffsetDriver(),
groundStation.getPrimeMeridianDriftDriver(),
groundStation.getPolarOffsetXDriver(),
groundStation.getPolarDriftXDriver(),
groundStation.getPolarOffsetYDriver(),
groundStation.getPolarDriftYDriver())) {
if (driver.getReferenceDate() == null) {
driver.setReferenceDate(observerTime);
}
}
TLEPropagator tlePropagator = TLEPropagator.selectExtrapolator(tle);
SpacecraftState spacecraftState = tlePropagator.propagate(observerTime);

// Build a range rate measurement. The parameters for range rate, sigma, and base weight are set to NaN due to the fact that they are not used in the calculation.
RangeRate rangeRate = new RangeRate(groundStation, observerTime, Double.NaN, Double.NaN, Double.NaN, true, new ObservableSatellite(0));

// Calculate range rate in meters per second. Iteration and evaluation are not used at all by estimateWithoutDerivatives,
// they are just past along to the EstimatedMeasurement object, so they can be retrieved by further processing.
// Formulas for the Doppler factor: v_range / c, where v_range is the range rate and c is the speed of light.
return rangeRate.estimateWithoutDerivatives(0, 0, new SpacecraftState[]{spacecraftState}).getEstimatedValue()[0] / Constants.SPEED_OF_LIGHT;
}
``````

I think I have figured it out myself, but please correct me if I am wrong.

When the satellite is moving toward the ground station, the range rate is negative, resulting in a negative value. Therefore, the solution is to manually add a negative sign in front of the final result to obtain the correct Doppler factor.

Yes, this is exactly what happens. The sign convention was chosen to be the velocity itself, not according to the frequency shift.

1 Like