Doppler shift using Orekit

SUIFENGSK · April 12, 2024, 1:37pm

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(
                new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                        Constants.WGS84_EARTH_FLATTENING,
                        FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
                new GeodeticPoint(Math.toRadians(gsPos.azimuth()), Math.toRadians(gsPos.elevation()), gsPos.distance()),
                "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;
    }

luc · April 12, 2024, 1:48pm

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.

SUIFENGSK · April 12, 2024, 2:21pm

Thank you for your reply!

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.

luc · April 12, 2024, 2:28pm

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.

SUIFENGSK · April 12, 2024, 2:49pm

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)),
                new GeodeticPoint(Math.toRadians(gsPos.azimuth()), Math.toRadians(gsPos.elevation()), gsPos.distance()), "GS"));
        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.

luc · April 12, 2024, 9:25pm

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

SUIFENGSK · April 12, 2024, 9:41pm

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?

luc · April 12, 2024, 9:56pm

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)?

SUIFENGSK · April 12, 2024, 11:16pm

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)),
                new GeodeticPoint(Math.toRadians(latitude), Math.toRadians(longitude), altitude), "GS"));
        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;
    }

luc · April 13, 2024, 7:29am

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

SUIFENGSK · April 13, 2024, 9:54am

It works now! Thank you!

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

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?
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?
What is the difference between using rangeRate.estimateWithoutDerivatives and rangeRate.estimate to estimate Doppler shift?"

luc · April 13, 2024, 10:28am

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.

SUIFENGSK · April 13, 2024, 10:50am

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

I use the following equation:

luc · April 13, 2024, 11:48am

Yes, because what is computed is vrange.