shashwatak / satellite-js

Modular set of functions for SGP4 and SDP4 propagation of TLEs.
MIT License
896 stars 144 forks source link

Error in calculating the Doppler Factor. #126

Open studentkra opened 7 months ago

studentkra commented 7 months ago

I am using the example from the Readme to calculate the Doppler Factor, slightly changing it because the observerCoordsEcf is not defined. I was guided by this theme:

https://github.com/shashwatak/satellite-js/issues/24

//dopplerFactor = satellite.dopplerFactor(observerCoordsEcf, positionEcf, velocityEcf);
velocityEcf = satellite.eciToEcf(velocityEci, gmst);
dopplerFactor = satellite.dopplerFactor(observerEcf, positionEcf, velocityEcf);
console.log(Math.round(dopplerFactor * centerFreq));

Everything seems to be working fine as long as the satellite is approaching. The frequency is decreasing, as it should be, and fully corresponds to the frequency shown by the Look4Sat program. But when the Doppler Factor is equal to one and the satellite begins to move away, the frequency begins to increase, it should continue to decrease. What is my mistake?

Update: I think there is a mistake in this function of the source code. An error in mathematical signs. The function never returns a value less than 1.

function dopplerFactor(location, position, velocity) {
    var mfactor = 7.292115E-5;
    var c = 299792.458; // Speed of light in km/s

    var range = {
      x: position.x - location.x,
      y: position.y - location.y,
      z: position.z - location.z
    };
    range.w = Math.sqrt(Math.pow(range.x, 2) + Math.pow(range.y, 2) + Math.pow(range.z, 2));
    var rangeVel = {
      x: velocity.x + mfactor * location.y,
      y: velocity.y - mfactor * location.x,
      z: velocity.z
    };
    function sign(value) {
      return value >= 0 ? 1 : -1;
    }
    var rangeRate = (range.x * rangeVel.x + range.y * rangeVel.y + range.z * rangeVel.z) / range.w;
    return 1 + rangeRate / c * sign(rangeRate);
  }
thkruz commented 7 months ago

@studentkra, looking through the history of the file, @Alesha72003 had recommended some changes (see https://github.com/shashwatak/satellite-js/pull/79) 3 years ago because they believed the previous logic was incorrect.

Looking through the tests that were created for this, it is solely testing that it can generate a doppler factor of 1. So your theory is entirely possible.

@Alesha72003 introduced an "mfactor" that I am not familiar with and do not see referenced anywhere on google. It is possible that m is a reference to a Russian word/letter which is why I am not finding it.

I am happy to find the issue and make the PR for you, but it would really help if you provided an (or ideally 3-5) example set(s) of input data and your expected output. Look4Sat is an awesome library btw (never seen it before), but it has been almost a decade since I have done anything with java so I wasn't able to solve it by just referencing their code.

I did find:

fun getDownlinkFreq(freq: Long): Long {
    return (freq.toDouble() * (SPEED_OF_LIGHT - distanceRate * 1000.0) / SPEED_OF_LIGHT).toLong()
}

fun getUplinkFreq(freq: Long): Long {
    return (freq.toDouble() * (SPEED_OF_LIGHT + distanceRate * 1000.0) / SPEED_OF_LIGHT).toLong()
}

Given that we are just calculating the dopplerFactor in satellite.js it would seem that the output you are expecting is:

...
return c + sign(rangeRate) / c;

Assuming that sign returns the correct + or -.

Alternatively we could just refactor this whole function to do a basic eci to range/az/el and then calculate doppler based on the the difference in range. It would be a lot simpler.

@studentkra let me know if you are willing to provide some examples so I can troubleshoot it for you. @ezze are you still around for this if I do a PR or should I direct him to ootk-core and make the change over there?

ezze commented 7 months ago

@thkruz Hi, mate! I still exist somewhere in the Universe so there are some chances to merge your PR. :)

But I believe that we should ask @shashwatak to give you all permissions to manage this repo and npm package. Unfortunally I can't do it by myself. You're definitely the best person who can maintain this one.

studentkra commented 7 months ago

Thanks for the answer! I replaced the function return with this:

if (rangeRate < 0) return 1 + rangeRate / c * sign(rangeRate);
if (rangeRate >= 0) return 1 - rangeRate / c * sign(rangeRate);

Everything seems to be working correctly now.

shashwatak commented 7 months ago

@thkruz collaborator invite sent!

thkruz commented 7 months ago

@studentkra I am digging into this code and I am having a hard time making sense of your solution.

Edit: It turns out the double negatives were tricking me. I think we can refactor your solution to ditch that sign function all together. Here is an example of two satellites, one that is moving towards us at 7km/s and one that is moving away at the same speed.

if (rangeRate < 0) return 1 - rangeRate / c; // 1 - (-7 / 299792.48) = ~1.00002334
if (rangeRate >= 0) return 1 - rangeRate / c; // 1 - (7 / 299792.48) = ~0.99997665

// So the return value could just be:
return 1 - rangeRate / c;

Does that make sense to you (and anyone else reading) and feel like an adequate solution?

As a side note, I found the mFactor by searching for the constant. It is the rotation speed of the earth and is being used to account for the observer moving while the earth spins.