[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Orekit Users] Right ascension and declination




Philippe Demoulin <Philippe.Demoulin@aeronomie.be> a écrit :

Thank you Luc for your fast response !

Hi Philippe,

Reading again what I wrote, I found that there is something
slightly wrong. In fact, there are two errors that compensate
each other in the code update I proposed.

I have computed the position of Jupiter and Earth in ICRF,
then subtracted them and implicitly used the fact ICRF and
EME2000 are aligned to declare the result of this subtraction
a position in EME2000 frame.

This is wrong.

ICRF and EME2000 are not really aligned. It seems we forgot
to update the javadoc comments in FramesFactory when we fixed
this in Orekit. There is a so called "frame bias" between
ICRF and EME2000, resulting from IAU 2000 resolution B1.6 which
can be found for example in IERS conventions 2010 on page 44.
The end of this resolution reads:

  In addition to that model are frame bias values between the
  mean pole and equinox at J2000.0 and the GCRS.

This frame bias is a constant rotation between GCRF and EME200,
with a total angle about 23 milliarcseconds (it is really a
combination of 3 small rotations, see EME2000Provider).

So IF we really wanted to compute Jupiter position in EME2000,
we should really do :

 Vector3D jupiterGCRF = body.
                        getPVCoordinates(date.shiftedBy(-delay), icrf).
                        getPosition().
                        subtract(earth);
Transform gcrf2EME2000 = FramesFactory.getGCRF().getTransformTo(eme2000, date);
 Vector3D jupiterEME2000 = gcrf2EME2000.transformPosition(jupiterGCRF);


HOWEVER, fixing this error and using the code above would not provide the same
results as Horizon. The reason is that there is a second error (but not in
Orekit) that compensate the first one. Despite Horizons (and also Spice/Naif) claim to compute things in frames aligned with EME2000, they really use frames alignes with ICRF. At some point in spice documentation (if I remember well), they say they use frames aligned with ICRF or EME2000 because they have the same orientation. This is probably because the bias is small. This puzzled us for a long time when doing
validation with Naif (see for example PredefinedIAUPolesTest.testNaif() in
Orekit unit tests). We finally asked to Naif developers and they confirmed
that what was really used was aligned with ICRF, but they never changed the
documentation and many people do not care about the small bias.

Hence, the updated code from my previous mail indeed does compute the same
thing as Horizon does, but the variable bodyEME2000 should really be renamed
bodyGCRF.

I will update Orekit javadoc to make it clear that solar system barycenter
(and Earth-Moon barycenter are really GCRF aligned).

best regards,
Luc


Best regards,

Philippe
-----------------------------------------------

On 07-Dec-17 14:08, Luc Maisonobe wrote:
Le 07/12/2017 à 12:00, Philippe Demoulin a écrit :
Hi everybody,

Hi Philippe,


I am completely new to Orekit, and I followed the 3-day course at CS
last week (thanks Luc !).

You're welcome.

I just started with a small program to
compute the right ascension and declination of Jupiter :

           AbsoluteDate date = new AbsoluteDate(2017, 12, 6, 10, 17, 0.0,
                   TimeScalesFactory.getUTC());

           CelestialBody body   = CelestialBodyFactory.getJupiter() ;
           Frame eme2000 = FramesFactory.getEME2000();

           Vector3D bodyInEME2000  = body.getPVCoordinates(date,
eme2000).getPosition();

           double rightAscension = bodyInEME2000.getAlpha() ;
           if (rightAscension < 0) rightAscension = rightAscension + 2.
* FastMath.PI ;
           double declination = bodyInEME2000.getDelta() ;
           System.out.format(Locale.US, "At %s UTC, %s RA = %10.5f°, DEC
= %10.5f° %n",
                   date, body.getName(),
FastMath.toDegrees(rightAscension),
                   FastMath.toDegrees(declination)) ;

The output of this program is :
"At 2017-12-06T10:17:00.000 UTC, Jupiter RA =  219.81633°, DEC =
-14.44226°"

Computing the same parameters with the JPL Horizons system
(https://ssd.jpl.nasa.gov/horizons.cgi) gives RA = 219.81392° and DEC =
-14.44186°, i.e. a difference of 8.7" in RA and -1.4" in DEC.

I tried to run the same from Horizon, ang I got the following:

 2017-Dec-06 10:17     14 39 15.36 -14 26 29.4  -1.72   5.51
6.24293497155197 -15.3706304  32.3247 /L   5.5594

So it is rather RA = 219.81400° (14H 39min 15.36sec) and DEC =
-14.4415°. My input on the HORIZONS system was:

 Ephemeris Type [change] :     OBSERVER
 Target Body [change] :        Jupiter [599]
 Observer Location [change] :  Geocentric [500]
 Time Span [change] :          Start=2017-12-06 10:10,
                               Stop=2017-12-06 10:20,
                               Step=1 m
Table Settings [change] :      defaults
Display/Output [change] :      download/save (plain text file)



Do you have any idea of the reason of these different answers ?

The reason for the discrepancy is that Horizons computes what the
observer will "see". As Jupiter is far from Earth, the light
coming from it takes almost 52 minutes to arrive at Earth center.
So what we see at 10:17 corresponds to a Jupiter position at about
09:25. So we have to take into account Jupiter motion.

The simplest way to do that is to compute Jupiter position at photon
travel start and Earth position at photo travel end, both in the Solar
System barycenter frame, and then to compute the difference of these
two positions. Here is how I did it:


 CelestialBody body  = CelestialBodyFactory.getJupiter() ;
 Frame         icrf  = FramesFactory.getICRF();
 Vector3D      earth = CelestialBodyFactory.
                       getEarth().
                       getPVCoordinates(date, icrf).
                       getPosition();

 double delay = 0;

 // search signal transit date,
 // computing the signal travel in inertial frame
 final double cReciprocal = 1.0 / Constants.SPEED_OF_LIGHT;
 double delta;
 int count = 0;
 do {
     double previous = delay;
     Vector3D bodyP  = body.
                       getPVCoordinates(date.shiftedBy(-delay), icrf).
                       getPosition();
     delay           = Vector3D.distance(bodyP, earth) * cReciprocal;
     delta           = FastMath.abs(delay - previous);
 } while (count++ < 10 && delta >= 2 * FastMath.ulp(delay));

 Vector3D bodyEME2000  = body.
                         getPVCoordinates(date.shiftedBy(-delay), icrf).
                         getPosition().
                         subtract(earth);

 double alphaDeg =
   FastMath.toDegrees(MathUtils.normalizeAngle(bodyInEME2000.getAlpha(),
                                               FastMath.PI));
 double deltaDeg = FastMath.toDegrees(bodyInEME2000.getDelta());
 System.out.format(Locale.US,
                   "At %s UTC arrival time, %s light time = %11.6fs," +
                   " RA = %s, DEC = %s%n",
                   date, body.getName(), delay,
                   toHMinSec(alphaDeg), toDegMinSec(deltaDeg));


with the two helper functions

    private String toHMinSec(double aDeg) {
        double s = FastMath.abs(aDeg / 15);
        int h = (int) FastMath.copySign(FastMath.floor(s), aDeg);
        s = (s - FastMath.abs(h)) * 60;
        int m = (int) FastMath.floor(s);
        s = (s - m) * 60;
        StringBuilder sb = new StringBuilder();
        Formatter f = new Formatter(sb);
        f.format(Locale.US, "%02dH %02dmin %6.3fsec", h, m, s);
        f.close();
        return sb.toString();
    }

    private String toDegMinSec(double aDeg) {
        double s = FastMath.abs(aDeg);
        int d = (int) FastMath.copySign(FastMath.floor(s), aDeg);
        s = (s - FastMath.abs(d)) * 60;
        int m = (int) FastMath.floor(s);
        s = (s - m) * 60;
        StringBuilder sb = new StringBuilder();
        Formatter f = new Formatter(sb);
        f.format(Locale.US, "%d° %02d' %6.3f''", d, m, s);
        f.close();
        return sb.toString();
    }


Notice the use of MathUtils.normalizeAngle (from Hipparchus) to ensure
the right ascension is centered around PI (i.e. it lies between 0 and
2PI).

An important part is that Jupiter position and earth position are *not*
computed at the same time, but they are nevertheless subtracted to
compute the light travel.

With this setup, I get:

At 2017-12-06T10:17:00.000 UTC arrival time, Jupiter light time =
  3115.254693s, RA = 14H 39min 15.359sec, DEC = -14° 26' 29.450''

So it is consistent with the results I get from Horizons.


For a real accurate observation, we should in fact also take into
account the aberration of light (i.e. the fact that Earth moves with
a non-negligible speed in the ICRF frame, so there is a velocity
composition to compute at photons arrival). You can see an example
of such computation in the Rugged library (which is built on top
of Orekit). The effect is however much smaller than the effect of
light time (order of magnitude of the milliarcsecond, whereas as
you noticed light time is several arc seconds).

best regards,
Luc



Thank you for your help.

Philippe Demoulin
Institut royal d'Aéronomie Spatiale de Belgique
Avenue circulaire 3, 1180 Uccle