Advanced Question Computing RA/Dec from LOS vector

[Gondos]

This is not strictly an SDK question since I'm working in the core but here it is anyway :
I'd like to compute the right ascension/declination of a radio signal received during a timestep, with sub timestep accuracy.
I have something that seems to work for Az/El, but I tried a couple of things for RA/Dec and cannot land anywhere close to the values shown in the Info dialog for planets.
It could be a reference frame issue or something like that, it also doesn't help that most formulas I found are not using the same model as Orbiter for planetary rotation.
It would be nice to have something that works on any planet if possible.
Here is what I've got for Az/El :

static Vector AddAberration(const Vector& n, const Vector& observerVel)
{
    Vector beta = observerVel / C0;
    double b2 = dotp(beta, beta);

    if (b2 < 1e-16)
        return n;

    DoubleDouble gamma = 1.0 / (dd_one - b2).sqrt();
    double ndotb = dotp(n, beta);
    Vector result = (n / (double)gamma + beta * (double)(1.0 + ((gamma - 1.0) / b2) * ndotb)) / (1.0 + ndotb);
    return result.unit();
}


struct AzEl {
    double azimuth;
    double elevation;
};

static AzEl GetAzimuthElevation(const RadioEvent *evt, const LightTimeSolution &sol, const CelestialBody *planet)
{
    AzEl azel;
    Vector photonDir = (sol.recvPos - evt->emitPos).unit();
    Vector apparentLOS = AddAberration(-photonDir, sol.recvVel);
    StateVectors sv = planet->InterpolateState(sol.u);
    Vector up = (sol.recvPos - sv.pos).unit();
    Vector pole = {sv.R.m12, sv.R.m22, sv.R.m32};  // planet rotation axis in global frame
    Vector east  = crossp(up, pole).unit();
    Vector north = crossp(east, up).unit(); // left-handed

    double z = std::clamp(dotp(apparentLOS, up), -1.0, 1.0);
    azel.elevation = asin(z);

    Vector horiz = apparentLOS - up * dotp(apparentLOS, up);
    double hlen = horiz.length();
    if (hlen < 1e-12) {
        azel.azimuth = 0.0;
        return azel;
    }
    horiz /= hlen;
    azel.azimuth = atan2(dotp(horiz, east), dotp(horiz, north));

    if (azel.azimuth < 0)
        azel.azimuth += PI2;
    
    return azel;
}

The goal is to be able to compute observables that could then potentially be fed into an orbit determination program (I also have 2-way range/range rate available), to check the consistency of a radio propagation model.
Thanks

 

[jarmonik]

I haven't run any tests but this is what I came up quickly.
Should precession be included ?

double z = dotp(apparentLOS, pole);
double Dec = asin(z);
 
Vector Aries = _V(1,0,0)
Vector X = crossp(pole, Aries).unit(); 
 
Vector horiz = apparentLOS - pole* dotp(apparentLOS, pole);
double hlen = horiz.length();
horiz /= hlen;
RA = atan2(dotp(horiz, X), dotp(horiz, Aries));

 

[Gondos]

I haven't run any tests but this is what I came up quickly.
Should precession be included ?

double z = dotp(apparentLOS, pole);
double Dec = asin(z);
 
Vector Aries = _V(1,0,0)
Vector X = crossp(pole, Aries).unit();
 
Vector horiz = apparentLOS - pole* dotp(apparentLOS, pole);
double hlen = horiz.length();
horiz /= hlen;
RA = atan2(dotp(horiz, X), dotp(horiz, Aries));

Thanks, the dec is somewhat close but the RA is off.
I get for example Venus from Earth :
23h 09m 12.53s and +13° 53' 12.81"
In the info window at the same time :
00h 50m 47.34s and +13° 51' 21.65"

For Mars :
07h 10m 57.85s -21° 45' 58.82" vs 16h 49m 2.78s -21° 45' 51.48"

The sum of the RA between the 2 solutions is close to 24h so there must be an inversion somewhere.
I'll try to see the impact of the retarded time if it can explain some of the difference.
Thanks!

 

[jarmonik]

Thanks, the dec is somewhat close but the RA is off.
I get for example Venus from Earth :
23h 09m 12.53s and +13° 53' 12.81"
In the info window at the same time :
00h 50m 47.34s and +13° 51' 21.65"

For Mars :
07h 10m 57.85s -21° 45' 58.82" vs 16h 49m 2.78s -21° 45' 51.48"

The sum of the RA between the 2 solutions is close to 24h so there must be an inversion somewhere.
I'll try to see the impact of the retarded time if it can explain some of the difference.
Thanks!

Yeah, it looks like RA is counted in opposite direction.
And there is then of-course the Earth Mean Equator and the current actual pole from oapiGetPlanetRotationMatrix() which might explain the difference in Dec. In which case the pole should be replaced by earth mean pole, obliquity of ecliptic is given in orbiter documentation. There must be constant defined for it too.

 

[n72.75]

This is not strictly an SDK question since I'm working in the core but here it is anyway :
I'd like to compute the right ascension/declination of a radio signal received during a timestep, with sub timestep accuracy.

Are we able to make API calls at the substep level?

 

[Gondos]

Are we able to make API calls at the substep level?

No, the integration substeps are not accessible. However at the end of state update (between pre and post callbacks, but after integration), there are some interpolation functions available inside the Orbiter core to get an estimation of the state at a time between the start and the end of the timestep.