Single-precision floating-point level rounding errors when using pmsafe about erfa

[PTW19]

I've experimented with various cut-off values but have yet to get convincing results. Although this fix may indeed end up being adopted, I'll need to understand precisely what's going on first.

[erykoff]

Now that I understand the issues, I'm happy to wait for the upcoming release (which sounds like fortuitous timing!). Also, it's my understanding that ERFA was supposed to be precisely matched to SOFA (warts and all). I certainly don't want to step into that!

Also, I understand that my patches are not to be used as-is, I just wanted to make sure that it was recorded precisely what code I tested. And I'm also assuming there will need to be corresponding updates to the SOFA FORTRAN code.

[PTW19]

I can't influence the date of the SOFA release much, but I should think it will be within two months.

The new versions, both C and Fortran, are done, and have passed the t_sofa_c/t_sofa_f tests. The Board will need to see them. If ERFA could use a sneak preview, best to contact me by email.

[mhvk]

Interesting. For further debugging, here's a more minimal version that does not use astropy at all (but uses pyerfa):

import matplotlib.pyplot as plt
import numpy as np
import erfa

pm_ra = np.arange(-5e-5, 5e-5, 1e-6)
starpm = erfa.pmsafe(np.pi, 0., pm_ra, 0., 0.1, 0.0, 2451545.0, 0., 2451545.0, 365.25)
plt.plot(pm_ra, starpm[2] - pm_ra)
plt.show()

[mhvk]

Looking through the actual starpm code, I think the reason for the fractional errors of order 1e-6 is that internally, all the positions are turned into double precision X, Y, Z coordinates, to which the velocities over a year are added. Now that itself should only give fractional uncertainties of 1e-11, but then there is the light travel time correction which does seem to be second order and thus possibly responsible: see

erfa/src/starpm.c

Line 141 in dcef368

tl2 = (-rdv + sqrt(rdv*rdv + c2mv2*r2)) / c2mv2;

(I have a feeling this has come up before and someone suggested a better approximation, but am not sure...)

[mhvk]

Separately, if we do find a better approximation, we'd need to bring this up with SOFA - in the end, ERFA is just a copy of that with a different license...

[cmsaunders]

@erykoff and I spent some time digging further into this yesterday and figured out that the problem behavior seems to be triggered here:

erfa/src/starpv.c

Line 203 in dcef368

w = (betsr != 0.0) ? d + del / betsr : 1.0;

depending on whether betsr evaluates to zero. We found that if you replace betsr != 0 with betsr > 1e-17, the RA proper motion becomes smooth and consistent. However, we aren't confident whether the values are now all correct or all incorrect with that change.

[mhvk]

@PTW19 - sorry to bother you again, but here is what would seem like another small buglet in SOFA. A possible solution at #95 (comment)

[PTW19]

Is it convenient to give me a program like the one in the original post but written in C and calling SOFA functions (or the rebadged ERFA equivalents)?

[cmsaunders]

This is what @erykoff put together yesterday:

#include <stdio.h>
#include "sofa.h"
​
​
int main(int argc, char **argv) {
    double lon[101];
    double lat[101];
    double pmra[101];
    double pmdec[101];
    double plx[101];
    double rv[101];
    double t1jd1, t1jd2, t2jd1, t2jd2;
​
    int i;
​
    double ra2[101];
    double dec2[101];
    double pmr2[101];
    double pmd2[101];
    double px2[101];
    double rv2[101];
​
    t1jd1 = 2451546.0;
    t1jd2 = -0.06145166083470002;
    t2jd1 = 2451911.0;
    t2jd2 = 0.18854833931867876;
​
    for (i=0; i<101; i++) {
        lon[i] = 3.14159265;
        lat[i] = 0.0;
        pmdec[i] = 0.0;
        plx[i] = 0.1;
        rv[i] = 0.0;
​
        pmra[i] = 1e-6*i + 5e-5;
    }
​
    for (i=0; i<101; i++) {
        iauPmsafe(lon[i],
                  lat[i],
                  pmra[i],
                  pmdec[i],
                  plx[i],
                  rv[i],
                  t1jd1,
                  t1jd2,
                  t2jd1,
                  t2jd2,
                  &ra2[i],
                  &dec2[i],
                  &pmr2[i],
                  &pmd2[i],
                  &px2[i],
                  &rv2[i]);
    }
​
    for (i=0; i<101; i++) {
        printf("Delta (%.10lf) = %.18lf\n", pmra[i], pmr2[i] - pmra[i]);
    }
​
    return 0;
}
​
/*
​
$CC sofa_bug.c -I./sofa/20210512/c/src/ -lm -L./sofa/20210512/c/src/ -lsofa_c
​
 */

[PTW19]

When I run that code I get output that starts:

Delta (0.0000500000) = -0.000000000000125000
Delta (0.0000510000) = -0.000000000000132651
Delta (0.0000520000) = -0.000000000000140608
Delta (0.0000530000) = -0.000000000000148877
Delta (0.0000540000) = -0.000000000000157464
Delta (0.0000550000) = -0.000000000000166375
Delta (0.0000560000) = -0.000000000000175616
Delta (0.0000570000) = -0.000000000000185193

If I plot the right hand column with Excel I get this:

What am I supposed to be seeing?

[erykoff]

That is different from what I see on both Linux x86 and macos arm64. To be clear I'm using gcc 11 or clang 14. I've downloaded sofa_c from https://www.iausofa.org/current_C.html#Downloads (sofa_c-20210512.tar.gz). I go into the src directory and type make. Then I go to the directory with this program, compile it and run ./a.out and get:

Delta (0.0000500000) = -0.000000000070688088
Delta (0.0000510000) = -0.000000000075014760
Delta (0.0000520000) = -0.000000000079514476
Delta (0.0000530000) = -0.000000000084190629
Delta (0.0000540000) = -0.000000000000157464
Delta (0.0000550000) = -0.000000000094085816
Delta (0.0000560000) = -0.000000000099311637
Delta (0.0000570000) = -0.000000000104727466
Delta (0.0000580000) = -0.000000000000195112
Delta (0.0000590000) = -0.000000000000205379
...

on linux or macos.
And interestingly the 5th element matches your output, but the first 4 do not.

[PTW19]

In brief as I have a broadband outage (thank you again Virgin Media) and having to use my iPhone.

All Windows10 + gfortran + gcc.

I have done a Fortran version of the program and it gives identical results to the C program (all on my platform). But when I compile it quad precision I get the phenomenon described in the open post:

I will investigate further.

[erykoff]

Thank you! Very interesting. The problem that I have found is (perhaps) around starpv.c where the line w = (betsr != 0.0) ? d + del / betsr : 1.0; depends on whether betsr compares to zero exactly or not. What is even more peculiar is that when it evaluates to zero is when we get the wrong answer. And at that point I've lost the thread because I don't fully understand what is happening in this function.

I'd be curious to hear if on windows (double or quad) if it's evaluating to zero or not. As I said I only tested linux and macos (but different architectures).

[cmsaunders]

I think the issue is that if betsr compares to zero, you end up doing the relativistic correction in the transverse direction, but not in the radial direction. I played around with forcing either both w and d to be one, and that makes the discrepancy go away (but I think results in a different wrong answer).

[PTW19]

I'll return to this tomorrow.

[PTW19]

I can at this stage confirm that there is a bug in iauStarpv, causing radically different behavior when the variable betsr (the radial component of the space velocity in units of c) is zero or non-zero. The reason the effect comes and goes is that this test example has zero radial velocity (a scalar - I do wish people would say recession speed); the radial component is calculated using the dot product of the direction and velocity vectors, and this is vulnerable to rounding errors, so doesn't come out to exactly zero every time and on every platform. Whether it does or not decides which of two expressions is used to calculate the scratch variable w; the two are radically different, and I'm beginning to suspect neither is correct.

Working out the correct coding is what Sherlock Holmes would have called a three-pipe problem, and I'm working on it. If there's anyone reading this who is an expert on special relativity, the algorithm is in Stumpff, P., Astron.Astrophys. 144, 232-240 (1985), and I'm trying to match the code and the paper without understanding either.

More in due course.

[erykoff]

Looking at the Stumpff paper, equations 59-60, it looks to me like the implementation is iauStarpv is incorrect. Specifically, in the paper delta_0 = (1 - \beta_0^2)^(1/2) - 1, which I think is incorrectly encoded in the line del = - w / (sqrt(1.0 - w) + 1.0);. Replacing this line with del = sqrt(1.0 - w) - 1.0; (assuming starpv del should be the same as Stumpff delta_0, and also observing that the temporary work variable w is the approximation of beta_0^2 in the iterative solution) is one part of my tentative fix. And the second part is replacing w = (betsr != 0.0) ? d + del / betsr : 1.0; with the much simpler w = d + del; (in line with the promise from Stumpff that this can be done in "only a few lines of FORTRAN code"). These changes seem to produce values that (at least in my short tests) "pass the smell test". That is, the final version of w (which according to the comments is the approximation to the inertial approximation of the radial velocity) does not noticeably depend on floating-point errors, smoothly varies with pmr, and comes close to but does not exceed 1.0.

With the caveats that I am no expert on special relativity, and I haven't smoked a pipe in over 10 years 😄.

[erykoff]

I may have spoken too soon. Or at least there needs to be a corresponding fix in iauPvstar that I haven't looked at in detail yet. (And it may be that I'm completely wrong!)

[mhvk]

I have a vague recollection this came up previously - that the current equation is there to avoid relying on 1-beta^2, since beta is generally very small.

[PTW19]

Specifically, in the paper delta_0 = (1 - \beta_0^2)^(1/2) - 1, which I think is incorrectly encoded in the line del = - w / (sqrt(1.0 - w) + 1.0);. Replacing this line with del = sqrt(1.0 - w) - 1.0; (assuming starpv del should be the same as Stumpff delta_0, and also observing that the temporary work variable w is the approximation of beta_0^2 in the iterative solution) is one part of my tentative fix.

That had me puzzled, too, but turns out to be correct (as I see mhvk has just pointed out). The expression coded is Stumpff's formula rearranged to minimize rounding errors: for small w the original formula involves subtracting 1 from a number close to one. Rather than trying to set this out as math, just run this sanity-check...

#include <stdio.h>
#include <math.h>
int main ()
{
double beta, beta2, a, b;
beta = 0.12345; /* or whatever /
beta2 = betabeta;
a = sqrt ( 1.0 - beta2 ) - 1.0; /* Stumpff /
b = - beta2 / ( sqrt(1.0-beta2) + 1.0 ); / as coded */
printf ( "%g\n", a-b );
return 0;
}

And the second part is replacing w = (betsr != 0.0) ? d + del / betsr : 1.0; with the much simpler w = d + del;

I had reached the same conclusion - that the problem is this line of code - but had yet to satisfy myself that the fix is identical to the Stumpff formula. Still working on that (but busy today).

I found old tests that demonstrated excellent round-trip agreement between STARPV and PVSTAR, which must mean both of these have the bug.

[erykoff]

Okay, here's a simple round-trip test that I wrote that touches on these issues. Note that my original code that @cmsaunders posted had a bug in that the tested pm values started at the maximum 5e-5 and went up from there (stupid sign error), but the same point still holds. In this simple test I'm just round-tripping iauStarpv/iauPvstar and for some values it works, and for other values the pmr and rv values go a little 🤪 (on my computer at least!).

#include <stdio.h>
#include "sofa.h"

int main(int argc, char **argv) {
    double ra[101];
    double dec[101];
    double pmr[101];
    double pmd[101];
    double px[101];
    double rv[101];

    int i;

    double ra2[101];
    double dec2[101];
    double pmr2[101];
    double pmd2[101];
    double px2[101];
    double rv2[101];

    double pv[2][3];

    for (i=0; i<101; i++) {
        ra[i] = 3.14159265;
        dec[i] = 0.0;
        pmd[i] = 0.0;
        px[i] = 0.1;
        rv[i] = 0.0;
        pmr[i] = 1e-6*i - 5e-5;
    }

    for (i=0; i<101; i++) {
        iauStarpv(ra[i], dec[i], pmr[i], pmd[i], px[i], rv[i], pv);
        iauPvstar(pv, &ra2[i], &dec2[i], &pmr2[i], &pmd2[i], &px2[i], &rv2[i]);
    }

    for (i=0; i<101; i++) {
        printf("%d: delta-lon = %.20lf\n", i, ra2[i] - ra[i]);
        printf("%d: delta-lat = %.20lf\n", i, dec2[i] - dec[i]);
        printf("%d: delta-pmra = %.20lf\n", i, pmr2[i] - pmr[i]);
        printf("%d: delta-pmdec = %.20lf\n", i, pmd2[i] - pmd[i]);
        printf("%d: delta-px = %.20lf\n", i, px2[i] - px[i]);
        printf("%d: delta-rv = %.20lf\n", i, rv2[i] - rv[i]);
    }

}

/*

$CC sofa_pvstar_roundtrip.c -I./sofa/20210512/c/src/ -lm -L./sofa/20210512/c/src/ -lsofa_c

*/

[erykoff]

Oh, and if I set the rv to something not zero (say 10 km/s) then the pmra round-trips okay, but there are residual errors on the returned rv.

[erykoff]

Oooh, this is potentially interesting. The forward transformation in starpv does the iterative solution to compute d and del. But the reverse transformation doesn't. And in the forward transformation d ends up as something not quite 1.0, but in the reverse it ends up exactly zero (no iteration) when betsr=0. So it doesn't transform back the same way (and this happens through the tangential component). I don't know if the iteration is coming to the wrong conclusion, or if it's supposed to iterate the same way in reverse, but I think this is what is causing my round-trip tests to fail.

[PTW19]

I had been assuming that you only need the iteration in one direction, IOW there's an uphill and a downhill transformation. But right now all bets are off. I've found some test programs from nearly 20 years ago (including comparisons with USNO NOVAS-C) that may have some clues. I'll report back in due course.

[PTW19]

My best shot at present is the following:

:

if (i >= IMAX) iwarn += 4;

/* Replace observed tangential velocity with inertial value. */
iauSxp(d, ust, ut);

/* Replace observed radial velocity with inertial value. */
w = (vst != 0.0) ? vst : 1.0;
iauSxp(d+del/w, usr, ur);

/* Combine the two to obtain the inertial space velocity. */
iauPpp(ur, ut, pv[1]);
:

I have still not mastered the Stumpff paper, partly because I keep getting velocity and speed mixed up. His expressions (59) for example are all scalar whereas iauStarpv is all about generating vectors.

See what you think. I'll need to make a similar modification in iauPvstar.

The numbers that the test program now produces start like this:

Delta (0.0000500000) = -0.000000000070688088
Delta (0.0000510000) = -0.000000000075014760
Delta (0.0000520000) = -0.000000000079514476
Delta (0.0000530000) = -0.000000000084190629
Delta (0.0000540000) = -0.000000000089046611
Delta (0.0000550000) = -0.000000000094085816
Delta (0.0000560000) = -0.000000000099311637
Delta (0.0000570000) = -0.000000000104727466
Delta (0.0000580000) = -0.000000000110336697
Delta (0.0000590000) = -0.000000000116142723

[erykoff]

Hmm. I'm not sure about the line w = (vst != 0.0) ? vst : 1.0; because this will evaluate to vst if it is very, very small, but to 1.0 if it is exactly zero. And that discontinuity worries me.

Meanwhile, for the reverse transform. For the forward transform I noticed that betsr starts out exactly 0, but betr is not because of the del factor which depends on the full speed (bett*bett + betr*betr). So therefore d != 1.0 when we get to the iauSxp(d, ust, ut) call. But in the reverse transform, betr is coming out exactly 0 (is this a bug? I'm not sure) and therefore d = 1.0 exactly, and we don't properly reverse the operation.

[PTW19]

Hmm. I'm not sure about the line w = (vst != 0.0) ? vst : 1.0; because this will evaluate to vst if it is very, very small, but to 1.0 if it is exactly zero. And that discontinuity worries me.

My suggested change is indeed incorrect (for a start vst is transverse not radial), but let me try to explain where the ugly 1.0 value came from.

Stumpff's expression for inertial radial speed (in units of c) is beta(inertial) = d*beta(observed)+delta. However, iauStarpv doesn't want adjusted radial speed (a scalar): it wants adjusted radial velocity (a vector}. The function works by resolving the original velocity vector into a radial component and a transverse component that are then each scaled before being added back together to reconstitute the final velocity vector. So what we are looking for is not Stumpff's speeds but the scaling factors we need to apply to our vector components. Thus setting a scaling factor to 1.0 to avoid a divide by zero would simply mean "leave the original vector unchanged", which sounds benign.

My current experimental coding looks similar to the released algorithm, which sounds reassuring except I am back to all the former numerical oddities for zero radial speed. Watch this space.

In case my contributions to this topic look shambolic, they are. Fortran code for the released algorithm existed before my electronic email records begin in 1999, so the early provenance is unknown. I may have to resort to understanding enough of the Stumpff paper to start again from scratch. If anyone reading this thread would like the challenge please pipe up!

[erykoff]

Thanks so much for the explanation, I'm starting to understand what is going on (though I don't know how to fix it).

On the forward starPv side, we have from Stumpff \beta_r^{inertial} = d*\beta_r^{obs} + \delta for the scalar. And then we multiply this by \vec{\beta_r^{obs}}/\beta_r^{obs} and that's where we get \vec{\beta_r^{inertial}} = (d + \delta/\beta_r^{obs})*\vec{\beta_r^{obs}} which is what is used. And where this goes wrong is when \beta_r^{obs} -> 0. At that point this goes off the rails.

I was thinking about just using \vec{\beta_r^{inertial}} = d*\vec{\beta_r^{obs}} + \delta but this doesn't seem to work (well, it certainly gives different answers). So there's much about the scalar/vector differences that I certainly don't understand!

On the reverse pvStar side, it seems to be doing something slightly different. At least, in the forward side then if \beta_r^{obs} == 0 (or approximately so) then after the iteration \beta_r^{inertial} != 0 because of the modification by delta which includes the full vector magnitude (both beta_r and beta_t). But in the reverse application, maybe due to rounding (?) then I see that \beta_r^{inertial} == 0 at the start, and therefore d = 1 + betr is identically 1.0 (which is different from the forward case because \beta_r^{inertial} != 0) and therefore the transverse direction doesn't properly round-trip.

(Hopefully these observations are closer to correct and therefore helpful in some way ...)

[PTW19]

Thanks.

I have a growing conviction that the STARPV algorithm is basically correct and that the problems are numerical. We need to bear in mind that the test program is actually displaying the difference between two very similar numbers, so the rounding errors are under a microscope.

Somewhere in the calculation - for example I've been looking at the separation of the space motion into transverse and radial components - there is an undue sensitivity to calculation precision that I hope a reformulation can cure. I'm comparing internal variables from both a C and a Fortran implemention, the latter in both double and quad compilations, but have yet to hit pay dirt. I don't even know which of several different sets of results is closest to the truth.

I've been referring both to the Stumpff A&A paper and this more concise one:

https://adsabs.harvard.edu/full/1986IAUS..114..193S

[erykoff]

You may be right. Upon further inspection, my simple round-trip tests only fail when rv=0.0. A possible fix, such as it is, is to demand that rv be greater than some epsilon value and add a new Note of warning when that happens (as currently happens with the parallax).

[PTW19]

I think I have diagnosed the problem.

As I said earlier, the principle STARPV and PVSTAR use is that the transformation between observed and inertial velocity vectors is a scaling factor close to unity. This scaling factor is applied separately to the tangential and radial components of the velocity vector before the two are recombined. The sensitivity of this procedure to rounding errors happens in cases where the radial component is very close to zero, and an enormous scaling factor has to be applied to reach the adjusted result. The original, apparently random, precaution of setting the scaling factor to unity to avoid a divide-by-zero was sound enough, because 1.0 times [0,0,0] is no different from (say) 1e40 times [0,0,0].

The answer is to perform the calculation in a far more obvious way. The Stumpff expression for beta_r computes the inertial radial speed directly and without any numerical complications; simply multiplying the unit vector towards the star by this quantity delivers the required adjusted radial velocity vector directly, ready to be combined with the adjusted tangential velocity vector.

I am working on revised coding for the STARPV and PVSTAR routines; in the meantime, these are the results I am getting from preliminary experiments:

DELTA (0.0000500000) = -0.000000000070280550
DELTA (0.0000510000) = -0.000000000074582278
DELTA (0.0000520000) = -0.000000000079056051
DELTA (0.0000530000) = -0.000000000083705245
DELTA (0.0000540000) = -0.000000000088533231
DELTA (0.0000550000) = -0.000000000093543383

[PTW19]

This is my current recoding of iauStarpv:

if (i >= IMAX) iwarn += 4;

/* Scale observed tangential velocity into inertial. */
iauSxp(d, ust, ut);

/* Compute inertial radial velocity vector. */
iauSxp((betsr-del)/d, x, ur);

/* Combine the two to obtain the inertial space velocity vector. */
iauPpp(ur, ut, pv[1]);

...and this is iauPvstar:

del = - w / (sqrt(1.0-w) + 1.0);

/* Compute inertial radial velocity vector. /
iauSxp ( dbetr+del, x, usr );

/* Scale inertial tangential velocity vector into observed */
iauSxp(1.0/d, ut, ust);

/* Combine the two to obtain the observed velocity vector. */
iauPpp(usr, ust, pv[1]);

My current results from the test program are as follows:

DELTA (0.0000500000) = -0.000000000071095627
DELTA (0.0000510000) = -0.000000000075447243
DELTA (0.0000520000) = -0.000000000079972901
DELTA (0.0000530000) = -0.000000000084676013
DELTA (0.0000540000) = -0.000000000089559992
DELTA (0.0000550000) = -0.000000000094628249

I will be interested to know how these new versions perform in your round-trip tests.

[PTW19]

Health warning - I think there's a missing factor of c in the radial velocity calculation. I'll be back.

[PTW19]

Latest versions:

STARPV

if (i >= IMAX) iwarn += 4;

/* Scale observed tangential velocity vector into inertial (au/d). */
iauSxp(d, ust, ut);

/* Compute inertial radial velocity vector (au/d). /
iauSxp(DC(d*betsr+del), x, ur);

/* Combine the two to obtain the inertial space velocity vector. */
iauPpp(ur, ut, pv[1]);

PVSTAR

del = - w / (sqrt(1.0-w) + 1.0);

/* Scale inertial tangential velocity vector into observed */
iauSxp(1.0/d, ut, ust);

/* Compute observed radial velocity vector. /
iauSxp(DC(betr-del)/d, x, usr);

/* Combine the two to obtain the observed velocity vector. */
iauPpp(usr, ust, pv[1]);

The results of running the test program are now:

Delta (0.0000500000) = -0.000000000000125000
Delta (0.0000510000) = -0.000000000000132651
Delta (0.0000520000) = -0.000000000000140608
Delta (0.0000530000) = -0.000000000000148877
Delta (0.0000540000) = -0.000000000000157464
Delta (0.0000550000) = -0.000000000000166375

IIRC, these are the same as what I get from the currently released SOFA code.

I have tested the two functions against Stumpff's Barnard's Star data and get his inertial values for proper motion and radial velocity, and an accurate round-trip:

10.3062860385891004917613308551090542
107.993949591739319780211341702375798
10.3100000000000004031609462578437016
108.000000000000000000000000000000074

[erykoff]

Utterly brilliant! That seems to have done it. I'm getting consistent results in all my tests, and they are invariant for small and 0 values around rv=0 km/s.

I note that github had some formatting problems with mistaking asterisks for markdown formatting.

For the record (since this is an erfa thread) the patches that I applied to erfa starpv.c and pvstar.c are here:
0001_relativity.patch
0002_relativity2.patch

Similarly, the patches for sofa starpv.c and pvstar.c are here:
0001_relativity_sofa.patch
0002_relativity2_sofa.patch

NB: I have only tested this for our original problem. I haven't tested for all possible values, but it makes sense to me and avoids the nasty divide-by-zero.

Thanks so much!

[PTW19]

I'm glad you like the fix. I'm pleased it came out so neatly in the end. Reassuringly, the new code still passes the existing t_sofa_c test.

We have a SOFA release imminent, which has up to now been almost entirely documentation enhancements; the new STARPV and PVSTAR will be an exception. Be aware that when the revised code appears it will contain other minor changes, so will differ from your patched version. Just a change of variable name and improvements to the commenting.

[mhvk]

Very nice indeed! If a new SOFA version is to be released soon, then simplest for ERFA would seem to just wait for that and release a new version based on it. Would that be OK for you, @erykoff?

[mhvk]

@PTW19 - thanks so much! We have an automatic process for turning SOFA into ERFA, so it should all be OK. Happy that this solution will be in a new version so soon!

[erykoff]

Any updates on a possible SOFA release to get this fixed? Or can we patch erfa? Thanks!

[PTW19]

Any updates on a possible SOFA release to get this fixed?

There is a SOFA release in the pipeline. I was surprised to find that my last record of the topic is an email dated 21 June asking when the release would occur. I will make enquiries. The version of starpv.c in the new release is attached. Patrick Wallace ----------------------------------------------------------------- TPOINT CONSULTING [phone] Sutton Wick Lane, Drayton [phone] Abingdon Oxfordshire OX14 4HH United Kingdom www.tpointsw.uk

…

_________________________________________________________________ Except where indicated, the contents of this email are for the named recipient(s) only and must not be made public or otherwise disclosed to third parties. ----------------------------------------------------------------- #include "sofa.h" #include "sofam.h" int iauStarpv(double ra, double dec, double pmr, double pmd, double px, double rv, double pv[2][3]) /* ** - - - - - - - - - - ** i a u S t a r p v ** - - - - - - - - - - ** ** Convert star catalog coordinates to position+velocity vector. ** ** This function is part of the International Astronomical Union's ** SOFA (Standards of Fundamental Astronomy) software collection. ** ** Status: support function. ** ** Given (Note 1): ** ra double right ascension (radians) ** dec double declination (radians) ** pmr double RA proper motion (radians/year) ** pmd double Dec proper motion (radians/year) ** px double parallax (arcseconds) ** rv double radial velocity (km/s, positive = receding) ** ** Returned (Note 2): ** pv double[2][3] pv-vector (au, au/day) ** ** Returned (function value): ** int status: ** 0 = no warnings ** 1 = distance overridden (Note 6) ** 2 = excessive speed (Note 7) ** 4 = solution didn't converge (Note 8) ** else = binary logical OR of the above ** ** Notes: ** ** 1) The star data accepted by this function are "observables" for an ** imaginary observer at the solar-system barycenter. Proper motion ** and radial velocity are, strictly, in terms of barycentric ** coordinate time, TCB. For most practical applications, it is ** permissible to neglect the distinction between TCB and ordinary ** "proper" time on Earth (TT/TAI). The result will, as a rule, be ** limited by the intrinsic accuracy of the proper-motion and ** radial-velocity data; moreover, the pv-vector is likely to be ** merely an intermediate result, so that a change of time unit ** would cancel out overall. ** ** In accordance with normal star-catalog conventions, the object's ** right ascension and declination are freed from the effects of ** secular aberration. The frame, which is aligned to the catalog ** equator and equinox, is Lorentzian and centered on the SSB. ** ** 2) The resulting position and velocity pv-vector is with respect to ** the same frame and, like the catalog coordinates, is freed from ** the effects of secular aberration. Should the "coordinate ** direction", where the object was located at the catalog epoch, be ** required, it may be obtained by calculating the magnitude of the ** position vector pv[0][0-2] dividing by the speed of light in ** au/day to give the light-time, and then multiplying the space ** velocity pv[1][0-2] by this light-time and adding the result to ** pv[0][0-2]. ** ** Summarizing, the pv-vector returned is for most stars almost ** identical to the result of applying the standard geometrical ** "space motion" transformation. The differences, which are the ** subject of the Stumpff paper referenced below, are: ** ** (i) In stars with significant radial velocity and proper motion, ** the constantly changing light-time distorts the apparent proper ** motion. Note that this is a classical, not a relativistic, ** effect. ** ** (ii) The transformation complies with special relativity. ** ** 3) Care is needed with units. The star coordinates are in radians ** and the proper motions in radians per Julian year, but the ** parallax is in arcseconds; the radial velocity is in km/s, but ** the pv-vector result is in au and au/day. ** ** 4) The RA proper motion is in terms of coordinate angle, not true ** angle. If the catalog uses arcseconds for both RA and Dec proper ** motions, the RA proper motion will need to be divided by cos(Dec) ** before use. ** ** 5) Straight-line motion at constant speed, in the inertial frame, ** is assumed. ** ** 6) An extremely small (or zero or negative) parallax is interpreted ** to mean that the object is on the "celestial sphere", the radius ** of which is an arbitrary (large) value (see the constant PXMIN). ** When the distance is overridden in this way, the status, ** initially zero, has 1 added to it. ** ** 7) If the space velocity is a significant fraction of c (see the ** constant VMAX), it is arbitrarily set to zero. When this action ** occurs, 2 is added to the status. ** ** 8) The relativistic adjustment involves an iterative calculation. ** If the process fails to converge within a set number (IMAX) of ** iterations, 4 is added to the status. ** ** 9) The inverse transformation is performed by the function ** iauPvstar. ** ** Called: ** iauS2pv spherical coordinates to pv-vector ** iauPm modulus of p-vector ** iauZp zero p-vector ** iauPn decompose p-vector into modulus and direction ** iauPdp scalar product of two p-vectors ** iauSxp multiply p-vector by scalar ** iauPmp p-vector minus p-vector ** iauPpp p-vector plus p-vector ** ** Reference: ** ** Stumpff, P., 1985, Astron.Astrophys. 144, 232-240. ** ** This revision: 2023 May 4 ** ** Copyright1 */ { /* Smallest allowed parallax */ static const double PXMIN = 1e-7; /* Largest allowed speed (fraction of c) */ static const double VMAX = 0.5; /* Maximum number of iterations for relativistic solution */ static const int IMAX = 100; int i, iwarn; double w, r, rd, rad, decd, v, pu[3], usr[3], ust[3], vsr, vst, betst, betsr, bett, betr, dd, ddel, ur[3], ut[3], d = 0.0, del = 0.0, /* to prevent */ odd = 0.0, oddel = 0.0, /* compiler */ od = 0.0, odel = 0.0; /* warnings */ /* Distance (au). */ if (px >= PXMIN) { w = px; iwarn = 0; } else { w = PXMIN; iwarn = 1; } r = DR2AS / w; /* Radial speed (au/day). */ rd = DAYSEC * rv * 1e3 / DAU; /* Proper motion (radian/day). */ rad = pmr / DJY; decd = pmd / DJY; /* To pv-vector (au,au/day). */ iauS2pv(ra, dec, r, rad, decd, rd, pv); /* If excessive velocity, arbitrarily set it to zero. */ v = iauPm(pv[1]); if (v / DC > VMAX) { iauZp(pv[1]); iwarn += 2; } /* Isolate the radial component of the velocity (au/day). */ iauPn(pv[0], &w, pu); vsr = iauPdp(pu, pv[1]); iauSxp(vsr, pu, usr); /* Isolate the transverse component of the velocity (au/day). */ iauPmp(pv[1], usr, ust); vst = iauPm(ust); /* Special-relativity dimensionless parameters. */ betsr = vsr / DC; betst = vst / DC; /* Determine the observed-to-inertial correction terms. */ bett = betst; betr = betsr; for (i = 0; i < IMAX; i++) { d = 1.0 + betr; w = betr*betr + bett*bett; del = - w / (sqrt(1.0 - w) + 1.0); betr = d * betsr + del; bett = d * betst; if (i > 0) { dd = fabs(d - od); ddel = fabs(del - odel); if ((i > 1) && (dd >= odd) && (ddel >= oddel)) break; odd = dd; oddel = ddel; } od = d; odel = del; } if (i >= IMAX) iwarn += 4; /* Scale observed tangential velocity vector into inertial (au/d). */ iauSxp(d, ust, ut); /* Compute inertial radial velocity vector (au/d). */ iauSxp(DC*(d*betsr+del), pu, ur); /* Combine the two to obtain the inertial space velocity vector. */ iauPpp(ur, ut, pv[1]); /* Return the status. */ return iwarn; /* Finished. */ /* Copyright2 */ }

[timj]

Looks like there was a release of SOFA in 2023-10-11. I have not confirmed that this included the fix.

[erykoff]

Ah! And the patch is included in the latest erfa release. Hooray!

[mhvk]

Ha, thanks for checking! And of course for reporting the original.