liberfa / erfa Goto Github PK

The pre-1972 _changes are approximated to $10^{-7}\kern{3 mu}\text{s}$, leading to strange results.

Example (After #91 fixed)

>>> from astropy.time import Time, TimeDelta
>>> t = Time('1971-12-31 23:59:60.107757999', format = 'iso', scale = 'utc', precision = 9)
>>> t
<Time object: scale='utc' format='iso' value=1971-12-31 23:59:60.107757999>
>>> t.tai
<Time object: scale='tai' format='iso' value=1972-01-01 00:00:10.000000002>

In theory, it should be earlier than 1972-01-01 00:00:10, for 1972-01-01 00:00:00 UTC = 1972-01-01 00:00:10 TAI and the step was -0.1077580 s.

And the method of detecting jump is also inexact.

/* Separate TAI-UTC change into per-day (DLOD) and any jump (DLEAP). */
   dlod = 2.0 * (dat12 - dat0);
   dleap = dat24 - (dat0 + dlod);

The end of a day may not be 24:00:00 but 24:00:00+dleap.

So $\text{dleap}=\text{dat24}-\left(\text{dat0}+\dfrac{24\kern{3 mu}\text{h}+\text{dleap}}{24\kern{3 mu}\text{h}}\text{dlod}\right)$, that is, $\text{dleap}=\dfrac{\text{ERFA\_DAYSEC}}{\text{ERFA\_DAYSEC}+\text{dlod}}(\text{dat24}-(\text{dat0}+\text{dlod}))$.

Interestingly, this inexact method just matches $10^{-7}\kern{3 mu}\text{s}$ precision, dleap is calculated correctly (just a coincidence).

So the table _changes and the method of detecting jump may need to be modified at the same time.

/* Dates and Delta(AT)s */
   static const eraLEAPSECOND _changes[] = {
      { 1960,  1,  1.41781799014226 },
      { 1961,  1,  1.42281799021726 },
      { 1961,  8,  1.37281798946726 },
      { 1962,  1,  1.84585798946726 },
      { 1963, 11,  1.94585799076726 },
      { 1964,  1,  3.24012999076726 },
      { 1964,  4,  3.34012999226726 },
      { 1964,  9,  3.44012999376726 },
      { 1965,  1,  3.54012999526726 },
      { 1965,  3,  3.64012999676726 },
      { 1965,  7,  3.74012999826726 },
      { 1965,  9,  3.84012999976726 },
      { 1966,  1,  4.31316999976726 },
      { 1968,  2,  4.21316999676726 },
      { 1972,  1, 10.0              },
      { 1972,  7, 11.0              },
      { 1973,  1, 12.0              },
      { 1974,  1, 13.0              },
      { 1975,  1, 14.0              },
      { 1976,  1, 15.0              },
      { 1977,  1, 16.0              },
      { 1978,  1, 17.0              },
      { 1979,  1, 18.0              },
      { 1980,  1, 19.0              },
      { 1981,  7, 20.0              },
      { 1982,  7, 21.0              },
      { 1983,  7, 22.0              },
      { 1985,  7, 23.0              },
      { 1988,  1, 24.0              },
      { 1990,  1, 25.0              },
      { 1991,  1, 26.0              },
      { 1992,  7, 27.0              },
      { 1993,  7, 28.0              },
      { 1994,  7, 29.0              },
      { 1996,  1, 30.0              },
      { 1997,  7, 31.0              },
      { 1999,  1, 32.0              },
      { 2006,  1, 33.0              },
      { 2009,  1, 34.0              },
      { 2012,  7, 35.0              },
      { 2015,  7, 36.0              },
      { 2017,  1, 37.0              }
   };

See astropy/astropy#10092 for details of this issue emerging in astropy.

If you imagine trying to update a source with proper motion of order 100 mas/yr, but with no distance information. erfa.starpm tries to treat the object as if it were at "infinite distance". In practice this means setting the parallax to 1e-7 arcseconds (10 Mpc).

In an intermediate step, the erfa.starpv routine calculates the cartesian position and velocity vectors of your star, but at this large distance the space velocity is >1/2 the speed of light, which erfa.starpv doesn't like. As a result it sets the velocity vector to 0.

Hence erfa.starpm doesn't update the position at all.

Perhaps objects with no distance information should be special cased in erfa.starpm?

Specifically, could link to the Julia implementation at JuliaAstro mentioned in #67 (comment)

Can anyone let me know if ERFA and SOFA are thread safe?
I could not find this information while googling.

Thank you,

Currently leap seconds are hard-coded into the C, as SOFA does. However, there are use cases (e.g. astropy/astropy#5783) where it would be useful to have the option for leap seconds to be added at runtime. This option should be added.

The main debate here is whether or not this is acceptable, as it's a break from how SOFA does it. I think yes, because it's a pretty important use case - users should be able to get leap seconds from online sources and update at runtime without having to recompile...

erfam.h offers a series of macros, basically astronomical constants, that are part of the public API of erfa (#11). Some of them have very short names. To avoid collisions I propose to prefix all the macros in erfam.h with ERFA_

DAU -> ERFA_DAU
WGS84 -> ERFA_WGS84

and so on.

That's what C libraries usually do. Macros in gsl start with GSL_, macros in math.h start with M_, in zlib start with Z_ etc

At https://github.com/liberfa/erfa/blob/master/src/atco13.c#L152
** eraAtioq quick ICRS to observed
should read
** eraAtioq quick CIRS to observed

Hello,

SOFA file cr.c for function iauCr (eraCr in ERFA) has an error in its documentation. This error is present in 20160503 and the previous version as well.

Line 19-20 is:

** Returned:
** char[] double[3][3] copy

It should be:

** Returned:
** c double[3][3] copy

This error was not caught by erfa_generator in astropy._erfa because this function is in VectorMatrix category and so is not processed.

I have notified SOFA at [email protected].

Regards,
Prasanth

See https://github.com/liberfa/erfa/blob/master/src/erfaversion.c#L10

@eteq - since you did the release, it's probably easiest for you to fix... I have a PR ready for the astropy side (indeed, that is how I found out...)

We should update to 12c; this supposedly includes #40 and #41.

Hi All,

I have a request from the SOFA board, of which I am a member. A user of SOFA has requested a change in the way that the SOFA header files are used to make SOFA easier to use in C++ code. I'm including the messages (from emails, with names removed) below.

Please let me know if you think there might be issues with these changes. Our initial looks suggest that we can fairly easily make the changes and that they would affect only a small number of users.

Thanks,

Scott Ransom [email protected]

We are proposing a change in the way SOFA manages its header files.

The present method is as follows:

• sofa.h contains function prototypes and #includes sofam.h.
• sofam.h contains typedefs plus macro definitions.
• Individual functions #include only sofa.h.
• User applications #include only sofa.h.

This turns out to cause problems for C++ developers, and the user
proposes the following (see note, appended):

• sofa.h contains function prototypes plus typedefs; it does not
  #include sofam.h.
• sofam.h contains only macro definitions.
• Individual functions #include sofa.h and, in individual cases,
  sofam.h as well.
• User applications usually #include only sofa.h, but if they have
  exploited the definitions in sofa.h (which nowhere in the
  documentation are they encouraged to do, and indeed in the AST
  cookbook are specifically told not to) must #include that as well.

Extracts from user's message:

I recently ran into a problem #including sofa.h into some code that
also uses a C++ template library. The problem was easily resolved
but I’ve included that actual error at the end of this message...it
might be somewhat intimidating for the novice programmer.

The root of the problem is that it is conventional to use short
upper case names for template parameters – in this particular case
“DC”, and sofam.c (#included by sofa.h) #defines DC to be a number
which causes the havoc below. As far as I can tell, the macro
definitions are not required to define the API, so could the
typedefs be moved to sofa.h and sofam.h not included? I can see
that this could break existing programs but the manual says that
the constants are defined in sofam.h so programs that use them
should be including that as well as sofa.h.

                 from FocalPlane.cpp:5:
/usr/include/eigen3/Eigen/src/LU/arch/Inverse_SSE.h: In static member
function ‘static void Eigen::internal::compute_inverse_size4<1, float,
MatrixType, ResultType>::run(const MatrixType&, ResultType&)’:
/home/dlt87/include/sofam.h:85:17: error: expected unqualified-id before
numeric constant
#define DAYSEC (86400.0)
                 ^
/home/dlt87/include/sofam.h:85:17: error: expected ‘)’ before numeric
constant
/home/dlt87/include/sofam.h:85:17: error: expected ‘)’ before numeric
constant
In file included from /usr/include/eigen3/Eigen/LU:40:0,
                 from /usr/include/eigen3/Eigen/Dense:2,
                 from AffineTransform.h:4,
                 from CoordinateFrames.h:16,
                 from CelestialPosition.h:9,
                 from FocalPlane.h:8,
                 from FocalPlane.cpp:7:
/usr/include/eigen3/Eigen/src/LU/arch/Inverse_SSE.h:84:5: error: ‘AB’ was
not declared in this scope
     AB = _mm_mul_ps(_mm_shuffle_ps(A,A,0x0F), B);
     ^~
/usr/include/eigen3/Eigen/src/LU/arch/Inverse_SSE.h:84:5: note: suggested
alternative: ‘dB’
     AB = _mm_mul_ps(_mm_shuffle_ps(A,A,0x0F), B);
     ^~
     dB
/usr/include/eigen3/Eigen/src/LU/arch/Inverse_SSE.h:88:86: error: cannot
convert ‘double’ to ‘__m128 {aka __vector(4) float}’ for argument ‘1’ to
‘__m128 _mm_sub_ps(__m128, __m128)’
  _mm_sub_ps(DC,_mm_mul_ps(_mm_shuffle_ps(D,D,0xA5),
_mm_shuffle_ps(C,C,0x4E)));
                                                                            
  ^
/usr/include/eigen3/Eigen/src/LU/arch/Inverse_SSE.h:105:45: error: cannot
convert ‘double’ to ‘__m128 {aka __vector(4) float}’ for argument ‘1’ to
‘__m128 _mm_shuffle_ps(__m128, __m128, int)’
     d = _mm_mul_ps(_mm_shuffle_ps(DC,DC,0xD8),AB);
                                             ^
/usr/include/eigen3/Eigen/src/LU/arch/Inverse_SSE.h:111:66: error: cannot
convert ‘double’ to ‘__m128 {aka __vector(4) float}’ for argument ‘1’ to
‘__m128 _mm_movelh_ps(__m128, __m128)’
     iA = _mm_mul_ps(_mm_shuffle_ps(B,B,0xA0), _mm_movelh_ps(DC,DC));
                                                                  ^
/usr/include/eigen3/Eigen/src/LU/arch/Inverse_SSE.h:112:80: error: cannot
convert ‘double’ to ‘__m128 {aka __vector(4) float}’ for argument ‘1’ to
‘__m128 _mm_movehl_ps(__m128, __m128)’
iA = _mm_add_ps(iA,_mm_mul_ps(_mm_shuffle_ps(B,B,0xF5),
_mm_movehl_ps(DC,DC)));
                                                                           
^
/usr/include/eigen3/Eigen/src/LU/arch/Inverse_SSE.h:138:49: error: cannot
convert ‘double’ to ‘__m128 {aka __vector(4) float}’ for argument ‘1’ to
‘__m128 _mm_shuffle_ps(__m128, __m128, int)’
     iC = _mm_mul_ps(A, _mm_shuffle_ps(DC,DC,0x33));
                                                 ^
/usr/include/eigen3/Eigen/src/LU/arch/Inverse_SSE.h:139:87: error: cannot
convert ‘double’ to ‘__m128 {aka __vector(4) float}’ for argument ‘1’ to
‘__m128 _mm_shuffle_ps(__m128, __m128, int)’
m_sub_ps(iC, _mm_mul_ps(_mm_shuffle_ps(A,A,0xB1),
_mm_shuffle_ps(DC,DC,0x66)));
                                         
                                  ^
Makefile.Release:653: recipe for target 'release/FocalPlane.o' failed

We will want a Makefile to compile ERFA and also allow make test to run the tests.

As detailed further in astropy/astropy#6239 and astropy/astropy#6240, it turns out to be rather complicated to adapt the changes from #42 into the python wrapper used in astropy (which we hope to someday make into a standalone Python library instead of only being accessible through astropy). So the question is: should the changes in #42 be altered to be more "SOFA-like" ? Primarily this would mean having the return values be arguments instead of actual return values (this is how SOFA "returns" multiple values from one function).

Also, it's not clear how best to copy over the config.h macros with the version info into astropy. It would be better if those could somehow get injected into the source code (although I'm not sure if that's possible with the C build machinery we're using here... I'm not familiar enough with autoconf to be sure).

Or if it turns out to be easier to solve astropy/astropy#6240 by making the generator machinery more flexible than to make the change here and deal with backwards-compatibility, that's fine too!

cc @olebole

As discussed in #49 (comment) and implemented in #58, there is now an API for updating the leap seconds. However, some compromises were made in this API, as it was basically done using Astropy as a "guinea pig". So this issue is a space to discuss practical drawbacks of the current API and would might be done to either improve or replace it.

There's not a clear end-point here, but tentatively I'd say if in 2 years time (i.e., the end of the LTS cycle of the Astropy version that this is used in), if we are happy with it we should declare the API "good". So by end of 2021 this can be closed if not before.

Given that astropy is going to update its copy of erfa to include 2015-Jun-30 leap second (astropy/astropy/pull/3794), we need a release of erfa with (at least) the same change. Package managers that use system erfa will need it.

It seems that the after the last update of the library (#26 ) it was never released

I think it would be good to rename the README.rst to README.

RST is quite Python specific, and more traditional people may get confused by this name. However, the content can/should be rst formatted.

This will fix #76

There is a new SOFA release 2012-12-02 available at www.iausofa.org/current_C.html
so, erfa should be upgraded as well.

It contains a number of new functions, summarized by the upstream-tracker

Related: #76 #82

#76 was fixed in SOFA 17a, however, it has brought another problem.

      leap = (fabs(dleap) > 0.5);

The threshold is 0.5 s, which means that mini-leaps, of which the maximum absolute value is 0.1077580 s, cannot be processed properly.

>>> from astropy.time import Time
>>> Time('1965-06-30 12:00:00', format = 'iso', scale = 'utc')
<Time object: scale='utc' format='iso' value=1965-06-30 11:59:59.950>

Considering that the minimum absolute value of UTC mini-leaps is 0.005 s (is there such a serious rounding error?), it is possible to determine an appropriate threshold.

According to Vallado's book (Fundamentals of Astrodynamics and Applications, 4th Ed.), the Fundamental arguments can be calculated using equation 3-82 (page 225).

In this equation, the second item is 1325*r*T, 99*r*T..., r = 360°,
The following code is how erfa calculate Nutation, IAU 1980 model, in function eraNut80:

/* Mean longitude of Moon minus mean longitude of Moon's perigee. */
   el = eraAnpm(
        (485866.733 + (715922.633 + (31.310 + 0.064 * t) * t) * t)
        * ERFA_DAS2R + fmod(1325.0 * t, 1.0) * ERFA_D2PI);

/* Mean longitude of Sun minus mean longitude of Sun's perigee. */
   elp = eraAnpm(
         (1287099.804 + (1292581.224 + (-0.577 - 0.012 * t) * t) * t)
         * ERFA_DAS2R + fmod(99.0 * t, 1.0) * ERFA_D2PI);

/* Mean longitude of Moon minus mean longitude of Moon's node. */
   f = eraAnpm(
       (335778.877 + (295263.137 + (-13.257 + 0.011 * t) * t) * t)
       * ERFA_DAS2R + fmod(1342.0 * t, 1.0) * ERFA_D2PI);

/* Mean elongation of Moon from Sun. */
   d = eraAnpm(
       (1072261.307 + (1105601.328 + (-6.891 + 0.019 * t) * t) * t)
       * ERFA_DAS2R + fmod(1236.0 * t, 1.0) * ERFA_D2PI);

/* Longitude of the mean ascending node of the lunar orbit on the */
/* ecliptic, measured from the mean equinox of date. */
   om = eraAnpm(
        (450160.280 + (-482890.539 + (7.455 + 0.008 * t) * t) * t)
        * ERFA_DAS2R + fmod(-5.0 * t, 1.0) * ERFA_D2PI);

My question is why use fmod for the term 1325.0 * t , 99.0 * t , ... ?
According to my knowledge fo c++, fmod function is "Computes the floating-point remainder of the division operation x/y.", so the result is less than 1.0, inconsistent with the equation above.

SOFA is very unusual for a C project in that it includes its documentation in the .c implementation files, not in the .h header files. This is the opposite of standard procedure in C, where generally a system library may have the header files (which also serve as usage documentation) and compiled binaries, but the original .c files are not available.

• IDEs/editors that provide pop-up help about functions are unable to do so.
• Tools that want to generate wrappers and include the documentation, or use the extra information about in/out variables contained therein, can not do so without the full source.

At the end of the day, this is probably a minor point relative to the effort it would take to fix it, but its rather irksome.

Whith this PR /pull/15 we are providing sonames that are wrong and break applications using erfa shared library. With the PR the soname is 1.0.0 and it should be 0.0.1.

The derived shared library is liberfa.so.1 instead liberfa.so.0, and everything linked against liberfa.so.0 gets broken.

@eteq I don' think you can get the soname version from the package version

Mostly procedural since #52 actually has the update and I'm about to do the release with it merged... But at least this serves as an FYI for @timj @david-terrett and @sergiopasra that I'm doing the release!

As discussed in #35, there's currently no check that the src/Makefile.am is kept up-to-date if new files are added.

Probably the best/easiest thing (as suggested in #35 (comment)) is to write a script that checks this. This should also get run as part of the travis tests that run on new PR's to make sure they add anything that might be needed.

README.rst still references version 1.6.0. It would be great to know what 1.7.0 -- the actual current release -- has changed.

For dates in the range where the UTC scale was variable (around 1961? to 1972), the SOFA routines iauTaiutc and iauUtctai do not round trip. In other words after converting from TAI => UTC => TAI, the end result is not identical to the input. This seems like an inconsistency, but perhaps I'm not understanding some subtlety in this situation?

Below is a demonstration of the issue.

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

int main()
{
    int j, iy, im, id, ihmsf[4];
    double a1, a2, u1, u2, a1_rt, a2_rt;
    /* Encode TAI date and time into internal format. */
    j = iauDtf2d ( "TAI", 1965, 1, 1, 0, 0, 33.7, &a1, &a2 );
    if ( j ) return 1;

    /* Transform TAI to UTC. */
    j = iauTaiutc ( a1, a2, &u1, &u2 );
    if ( j ) return 1;

    /* Transform UTC back to TAI */
    j = iauUtctai ( u1, u2, &a1_rt, &a2_rt );
    if ( j ) return 1;

    /* Print round trip difference */
    printf("TAI -> UTC -> TAI round trip difference: %f seconds\n", 
           ((a1_rt - a1) + (a2_rt - a2))*86400);

    /* Output: TAI -> UTC -> TAI round trip difference: -0.002592 seconds */

    return 0;
}

I'm trying to compile erfa for new Fedora that comes with gcc 6. This is causing a problem in one of the tests. In particular I see this problem in i386 (in x86_64 the test passes)

eraFk52h failed: rv want -7.6000000940000251859 got -7.6000000939851357629 (1/5.1e+11)
t_erfa_c validation failed!

The radial velocity is computed here

When compiling the current version, the SONAME is built from the package version, not from the LIB_VERSION:

$ readelf -d ./src/.libs/liberfa.so |grep SONAME
 0x000000000000000e (SONAME)             Library soname: [liberfa.so.1]

Looking into configure.ac however gives:

So it should actually be liberfa.so.7.
While this is not an issue now, we may face backward incompatible changes in future, #58.
(Just curious: where does the current number come from?)

With the current src/Makefile.am, erfam.h is installed alongside erfa.h when the users does make install

Is this correct or erfam.h is an internal header, used for compilation only?

@eteq - I see that on zenodo this is still at 1.4. Would be good to auto-update it! Unfortunately, I think you'll have to adjust the settings for that...

As discussed while updating erfa to 1.3 in astropy, it would be good if we could expose the version number in the erfa library, as then we could have workarounds be used only when actually needed. See astropy/astropy#5418 (comment)

Once #7 and #14 have been either merged or pushed to later, we can do a release that integrates the new SOFA version 20131202_b

cc @sergiopasra

The only version of erfa in Zenodo is 1.4. We should upload the next releases or enable the link between Zenodo and github so that every release in github gets an entry in Zenodo

As noted by @rmathar (#44 (comment)), SOFA version 13 is out, and we should update ERFA to it, since it solves #40 and #41. It is also a good opportunity to finally get in #42...

This is the last step discussed in #49 (comment) . To summarize, now that #58 has provide a possible leap second interface (although emphasis on "possible" because of #61), we should figure out a way to implement @olebole's idea in #49 (comment), which is also connected with @nikkelj and @rmathar's ideas in #43.

Perhaps @rmathar's code in #43 (comment) can be recycled here?

This needs to be populated with "derived" SOFA code from the liberfa/erfa-fetch repository. It should probably go in a sec directory, not straight into root.

Is there any hope of building and distributing Windows shared library binaries along with the source tarballs for each release? This would make it easier to support Windows in ERFA.jl. Ref: JuliaAstro/ERFA.jl#17

There were some confusions about the SONAME version in erfa-1.0.0, and we are now going to have Fedora (@sergiopasra) released with SONAME 0.0.0, and Debian with 1.0.0.

Please release a version 1.0.1 ASAP with (just) the correct number so that we can update our packages and have them using the same SONAME on both distributions.

See: http://www.iausofa.org

Once we have the finalized notice, we should add a file containining the letter from the SOFA board that indicates we have the right to make this derived version.

Over at JuliaAstro/AstroBase.jl we are currently in the process of porting several ERFA functions to Julia. The rest of AstroBase.jl is published under the MIT license. We are now asking ourselves how we can acknowledge the ERFA heritage. According to the internet, it is not really clear to what degree a port is a derivative work as long as it is not an automatic or line-by-line translation.

My new-jerk reaction would be to adapt our LICENSE file like shown below and put a reference to the respective ERFA source code in the function docstring:

$ASTRO_BASE_LICENSE (MIT)
--
Some functions in AstroBase.jl are based on ERFA which is published under the following license:
$ERFA_LICENSE

But IANAL 🤷‍♂ and thus I would like to ask you folks how you want this to be handled.
Many thanks in advance!

Nothing super-urgent in it, but should stay up to date.

As I understand it, erfa's time routines currently represent time as a pair of double-precision floating point numbers. The first is the day part, and the second is a fractional sub-day part, of a Julian date.

I think that using two 64-bit integers would be superior. The first integer could be the Julian date kiloseconds, and the second integer could be the femtoseconds of remainder (a femtosecond is 1E-15 seconds, or 1 million nanoseconds).

This system:

• has a superior range for the values that can be expressed at both the top end (10x more years covered) and bottom end (10,000x finer precision).
• would use exact arithmetic, avoiding an entire class of nasty bugs, and making numerical algorithms much simpler
• would probably be slightly faster
• would permit stable conversion algorithms, since round-tripping from (say) UTC to TAI and back again (as described here) would not incur any floating point rounding error
• would prevent user-side confusion like that evinced in astropy/astropy#14455

I realize that this is not a simple change to implement! It would require a rewrite of every aspect of erfa's time-handling code. But I believe it could be a superior system, and since erfa is so fundamental to astronomy, this improvement could have a very wide impact.

Ideally, this change in internal representation would not cause much turmoil for users. Most astronomers interact with erfa through other tools like Astropy which could perhaps manage a transition.

An example value

As an example, consider the timestamp "2023-03-19 13:52:14.53526." In erfa, that gets represented as two 64-bit floats:

2460023.0            # Julian days
0.077982964907407414 # Julian fractional days

Note that due to floating point error, the fractional day is ± 4.9e-18 days. The whole day part is exact, however.

I propose it be represented as two 64-bit integers; the first signed, and the second unsigned:

212545993          # Julian kiloseconds
937728179931640576 # Julian fractional kiloseconds, in femtoseconds

Using Python as a shorthand, here is an example of conversion from floating point to integer:

def convert(jd, jd_fraction):
     ksec1, remainder1 = divmod(jd * 86400, 1000)
     ksec2, remainder2 = divmod(jd_fraction * 86400, 1000)
     jd_kilosec = int(ksec1) + int(ksec2)
     jd_kilosec_fraction = int((remainder1 + remainder2) * 1e15)
     return jd_kilosec, jd_kilosec_fraction

Range of values calculations

64-bit floating point numbers have a maximum whole-number value of 2^53+1; beyond that, the value will no longer be a whole number. So the maximum Julian Day which can be represented is 9,007,199,254,740,993, which is about 2.47E13 years.

A signed 64-bit integer for kiloseconds could represent whole numbers up to 2^63 kiloseconds, which is 2.92E14 years, so this permits about 10x larger range at the top end.

The fractional floating point value for days is precise down to 2^-53 days, which is 1.11E-16 days, or 0.009 nanoseconds.

An unsigned 64-bit integer for the fractional part would need to represent values in [0, 1) kiloseconds, or [0, 1000) seconds. 1 kilosecond is 1e18 femtoseconds, which is less than 2^64, so each integer value can be as single femtosecond and the whole fraction's range is covered. This means the precision floor is 1 femtosecond which is 0.000001 nanoseconds, about 10,000x more precise than erfa's current implementation.

For SOFA 18/erfa 2.0, there was a change in how the header files were dealt with. In particular, the structs defined in sofam.h were moved to sofa.h, and sofam.h was no longer automatically included in sofa.h. We concluded this was an API change based on our own definition (see #84 (comment)),

However, as noted in #84 (comment), it does not see SOFA agrees, as http://www.iausofa.org/current_changes.html states:

The consequence of the change for developers is that any application that uses constants from sofa.h will now need an explicit #include <sofam.h>. It should, however, be noted that the contents of sofam.h are not formally part of the SOFA API, and so such an application is vulnerable to future changes. Rather than adding an #include <sofam.h> to applications, it would be better to replicate the needed macros in the user code.

So, we now have a choice: either we follow suit and no longer count the macros in erfam.h as part of the API, or we slightly diverge from SOFA and continue to count them. Obviously, in principle best to stay as close as possible, though frankly the SOFA rationale surprises me. Surely, one would want compiled code that uses SOFA to use constants like the tropical year instead of copying them, so that if there is a correction (of a typo, say), that automatically gets propagated. Indeed, for this very reason, we export the values in erfam.h in pyerfa. (And even for the more proper macros, such as DNINT, it seems correct to count them as API; after all, those were changed because they were not quite right initially). cc @scottransom who may have more insight in the SOFA thinking.

p.s. Note that I think in practice this is unlikely to matter much - changes to the macros are few and far in between.

There is a May 2016 update to SOFA. From the release notes, this seems to include new routines to go from equatorial to ecliptic, so quite possibly relevant. cc @eteq.

In astropy/astropy#5425 and astropy/astropy#5466 it was noted that for i386, the astropy tests of erfa fail; this seems to be due to a compiler optimization problem (see astropy/astropy#5425 (comment)), which may be solved either by not including -O2 or by adding some other options (see #33 (comment)). I don't know to what extent we can control that here, but it might at least be good to add a test case which ensures that output is correct. Possibly this should be brought to the attention of SOFA as well.

It seems that for a time just before the day boundary, it's possible to get the wrong day / hour / min / second representation:

In [11]: t = Time('1960-01-02') + 0.0

In [12]: t.jd1
Out[12]: 2436935.0

In [13]: t.jd2
Out[13]: 0.4999999999999998

In [14]: erfa.d2dtf('TAI', 6,  t.jd1, t.jd2)
Out[14]: (1960, 1, 2, (0, 0, 0, 0))

In [15]: erfa.d2dtf('UTC', 6,  t.jd1, t.jd2)
Out[15]: (1960, 1, 1, (23, 59, 60, 0))

This doesn't seem to occur for more "modern" dates like 1972-01-02.

As @taldcroft mentioned in #7, we may be nearly for a 1.0 release of ERFA. However, one important question needs to be resolved before we can do a release: where will we host the release tarballs? A few options:

1. Create github branches for each release. Each branch would have a single commit with just the version in question. That way github automatically hosts via the "Download ZIP" button. I don't like this because it's sort of a misuse of git, and the number of branches will decrease without bound.
2. Same as 1, but use tags that are detached from a branch instead of branches. I'm not sure it's actually possible to push detached commits to github, but I could experiment with it if this is judged the best way.
3. Host on sourceforge/google code/etc.
4. Find some generous soul or institution to host the tarball.

Any possibilities I missed? Obviously 4 would be the most straightforward, so any ideas on that front?

cc @taldcroft @astrofrog @olebole @timj @RobertLuptonTheGood @perrygreenfield

Just what the title says - this is particularly relevant once #5 is merged, as that needs to know the scheme to produce the correct files.

We were having problems testing the Rubin Science Pipelines PM implementation with astropy/erfa and we traced it down to what we believe are floating-point errors in erfa. The simplest repro that I have come up with is here. It may be that we (or astropy) are misunderstanding how the routines are supposed to be used, but the oscillation between very small (double-precision) residuals and larger (single-precision) residuals is surprising. I note that this gets worse the smaller the parallax is (which can produce a warning at very small parallax). cc @cmsaunders, @timj .

from astropy.coordinates import SkyCoord
import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt
import astropy.time
from astropy.coordinates.representation import (
    RadialDifferential,
    SphericalDifferential,
    SphericalRepresentation,
)
from astropy.coordinates.builtin_frames.icrs import ICRS
import erfa
from astropy.time import Time


testPM_RA = np.arange(-5e-5, 5e-5, 1e-6) * (u.radian/u.yr)
testPM_Dec = np.zeros(len(testPM_RA)) * (u.radian/u.yr)
testRA = 180 * np.ones(len(testPM_RA)) * u.degree
testDec = 0 * np.ones(len(testPM_RA)) * u.degree

cat = SkyCoord(testRA, testDec,
               pm_ra_cosdec=testPM_RA, pm_dec=testPM_Dec)

t1 = Time("J2000")
t2 = t1 + 1 * u.yr

t1 = t1.tdb
t2 = t2.tdb

icrsrep = cat.icrs.represent_as(SphericalRepresentation, SphericalDifferential)
icrsvel = icrsrep.differentials["s"]

starpm = erfa.pmsafe(
    icrsrep.lon.radian,
    icrsrep.lat.radian,
    icrsvel.d_lon.to_value(u.radian / u.yr),
    icrsvel.d_lat.to_value(u.radian / u.yr),
    np.zeros(len(cat)) + 0.1,
    0.0,
    t1.jd1,
    t1.jd2,
    t2.jd1,
    t2.jd2,
)

print(starpm[2] - icrsvel.d_lon.value)

plt.plot(icrsvel.d_lon.value, starpm[2] - icrsvel.d_lon.value, 'r-')
plt.show()

Prints:

[ 1.25000000e-13  1.17648987e-13  6.25403038e-11  1.03822983e-13
  9.73359848e-14  5.15316310e-11  4.81719691e-11  4.49616002e-11
  4.18971311e-11  3.89751689e-11  6.40000002e-14  3.35451933e-11
  3.10303938e-11  2.86445292e-11  4.66559978e-14  2.42460328e-11
  2.22266148e-11  2.03225598e-11  3.27679981e-14  1.68469662e-11
  1.52686417e-11  2.43890025e-14  1.24139720e-11  1.11308409e-11
  9.93932156e-12  1.56249999e-14  7.81754604e-12  6.88050386e-12
  1.06479970e-14  9.26100085e-15  4.52404333e-12  6.85900176e-15
  3.29802774e-12  2.77832829e-12  4.09599722e-15  1.90858098e-12
  1.55174705e-12  2.19699679e-15  1.72800312e-15  7.52687831e-13
  5.65505516e-13  4.12253527e-13  5.11999229e-16  1.93968399e-13
  2.15999331e-16  1.24999193e-16  6.40009625e-17  2.70008693e-17
  7.99980267e-18  9.99922394e-19  0.00000000e+00 -9.99922394e-19
 -4.52404448e-15 -2.70000222e-17 -3.61923550e-14 -7.06881927e-14
 -2.15998484e-16 -1.93968398e-13 -2.89538830e-13 -4.12253527e-13
 -5.65505518e-13 -7.52687834e-13 -1.72800312e-15 -2.19699848e-15
 -1.55174705e-12 -1.90858098e-12 -4.09599722e-15 -2.77832829e-12
 -3.29802774e-12 -6.85900176e-15 -4.52404333e-12 -9.26100085e-15
 -1.06479970e-14 -6.88050386e-12 -7.81754604e-12 -1.56249999e-14
 -9.93932156e-12 -1.11308409e-11 -1.24139720e-11 -2.43890025e-14
 -1.52686417e-11 -1.68469662e-11 -3.27680049e-14 -2.03225598e-11
 -2.22266148e-11 -2.42460328e-11 -4.66560045e-14 -2.86445292e-11
 -3.10303938e-11 -5.93190048e-14 -3.61923207e-11 -3.89751689e-11
 -4.18971311e-11 -4.49616002e-11 -4.81719691e-11 -5.15316310e-11
 -5.50439787e-11 -1.03822997e-13 -6.25403038e-11 -6.65310671e-11
 -1.25000006e-13]

And the plot is: