diff --git a/doc/6-ImplementationStatus.md b/doc/6-ImplementationStatus.md
index 0c9444e88..50bf037d9 100644
--- a/doc/6-ImplementationStatus.md
+++ b/doc/6-ImplementationStatus.md
@@ -173,7 +173,7 @@ spellings.
 * `DateDash`
 * `DateDot`
 * `DateForeground`
-* `DateMenu`
+* `DatesMenu`
 * `DateSeparatorCommand`
 * `DateSlash`
 * `DateSpace`
@@ -1124,7 +1124,7 @@ The following is an extensive list of commands.
 * `DateDash`
 * `DateDot`
 * `DateForeground`
-* `DateMenu`
+* `DatesMenu`
 * `DateSeparatorCommand`
 * `DateSlash`
 * `DateSpace`
@@ -2073,7 +2073,7 @@ The following is an extensive list of commands.
 * `DateDash`
 * `DateDot`
 * `DateForeground`
-* `DateMenu`
+* `DatesMenu`
 * `DateSeparatorCommand`
 * `DateSlash`
 * `DateSpace`
diff --git a/doc/commands/system.md b/doc/commands/system.md
index 7afbb737b..dada4f638 100755
--- a/doc/commands/system.md
+++ b/doc/commands/system.md
@@ -1,6 +1,89 @@
-# Time, Date, Alarms and System Commands
 
-## DateMenu
+# Date and Time
+
+
+## Date entry / edit format
+
+A whole number with a `_date` unit attached is interpreted as a date according to the customary format `YYYYMMDD`.
+
+However, the year need not be a 4 digit number.
+It can a smaller, bigger, zero, even a negative number and leading zeros can be ommitted
+(as such it need not fit the `YYYY` format).
+In case of a negative number only the year is treated as such.
+Months and days, even in a negative year, run from low (positive numbers) to high within the year.
+
+Note: the date format is intentionally different from the format on the HP-48.
+
+Note: not all commands can handle a negative year.
+The DMCP platform for SwissMicros calculators only handles positive years.
+This restriction bubbles up to several DB48x commands.
+
+
+## Date display format
+
+By default a Date is displayed on the stack like `Sun 30/Nov/2025`.
+The format can be adjusted by changing the [Date, Time flags](#date-time-flags) or
+[Date, Time settings](#date-time-settings).
+
+
+## Date+Time entry / edit format
+
+A decimal number with a `_date` unit attached is interpreted as a Date + Time
+according to the format `YYYYMMDD.HHMMSS`.
+
+
+## Date+Time display format
+
+By default a Date+Time is displayed on the stack like `Sun 30/Nov/2025, 16:43:10`.
+The format can be adjusted by changing the [Date, Time flags](#date-time-flags) or
+[Date, Time settings](#date-time-settings).
+
+
+## Time entry format
+
+This is based on the DB48x feature to quickly enter an angle value,
+with repeated use of the _,_ (_._) key.
+For example to enter, `2°03′05″`, 
+you can enter the sequence _2.3.5_.
+
+After the minutes, the second '.' will change the display to `<hours>°<minutes>′_dms`,
+assuming you are entering an angle.
+Proceed with the seconds.
+Optionally you can give a third '.' followed by a fraction of a second.
+Enter the fraction as numerator '.' (yes, a fourth '.'; not a '/') denominator.
+
+Finally change the unit from `_dms` to `_hms`:
+- With the cursor positioned between the last digit and the unit separator (`_`)
+you can use the `_hms` menu button (it replaces the `_dms`).
+- With the cursor positioned after the the unit separator (`_`)
+you can use the _O_ (_x10n_, DM32/42/42n _E_) key (bound to the Cycle command)
+to alternate between `_dms` and `_hms`.
+- Or you can use the `→HMS` command.
+But, that also pushes the value on the stack.
+- Or, in alpha mode, you can change the 'd' in an 'h'.
+
+The `_hms` (InsertHms) and the `→HMS` (ToHMS) commands can be found in the [TimeMenu](#timemenu).
+
+A second way to enter a time is as `<hours>.<decimal fraction of an hour>_hms`.
+
+
+## Time edit format
+
+In edit mode a time wil show as `HH°MM′SS″<fraction of a second>_hms`.
+A leading zero of the hours will be left out.
+If the fraction of a second is zero it will be left out.
+
+
+## Time display format
+
+On the stack a time wil show as `HH:MM:SS.<decimal fraction of a second>`.
+A leading zero of the hours will be left out.
+Note the ':' as separator.
+If the fraction of a second is zero the '.' and the fraction will be left out.
+The `_hms` unit will not show.
+
+
+## DatesMenu
 
 Show a softkey menu for date-related commands, including:
 
@@ -12,9 +95,12 @@ Show a softkey menu for date-related commands, including:
 * `JulianDayNumber`
 * `DateFromJulianDayNumber`
 
+**Access:** [Time](#TimeMenu) menu.
+
+
 ## TimeMenu
 
-Show a softkey menu for time-related commands:
+Show a softkey menu for time-related commands, including:
 
 * `Time`
 * `ToHMS`
@@ -28,18 +114,72 @@ Show a softkey menu for time-related commands:
 * `TimedEval`
 * `SetTime`
 
+**Access:** 🟦 _V_ (_TIME_) key.
 
-## Date format
 
-The date format is `YYYYMMDD`, with an optional fractional part defining the
-time, as in `YYYYMMDD.HHMMSS`.
+## Date, Time flags
+
+* `HideDate` / `ShowDate`
+* `HideTime` / `ShowTime`
+* `ShowSeconds` / `HideSeconds`
+* `Time12H` / `Time24H`
+* `HideMonthName` / `ShowMonthName`
+* `MonthBeforeDay` / `DayBeforeMonth`
+* `YearFirst` / `YearLast`
+* `TwoDigitYear` / `FourDigitYear`
+* `HideDayOfWeek` / `ShowDayOfWeek`
+* `BCECEyearNumbering` / `AstronomicalYearNumbering`
+
+**Access:** 
+
+
+## Date, Time settings
+
+* `DateSlash` / `DateDash` / `DateDot` / `DateSpace`
+* `DateSeparatorCommand`
+
+**Access:** 
+
+
+## BCECEyearNumbering
+
+Activate the B(efore) C(ommon) E(ra) / C(ommon) E(ra) year numbering.
+
+Display years on the stack with a BCE (for negative years) or CE (for positive years) suffix.
+This uses a non uniform timeline without a year = 0.
+The year numbering jumps from 1 BCE to 1 CE. So negative years are shifted.
+See: [Proleptic Gregorian calendar](https://en.wikipedia.org/wiki/Proleptic_Gregorian_calendar)
+
+**Access:** 
 
-Note: the date format is intentionally different from the format on the HP-48.
 
-## Time format
+## AstronomicalYearNumbering
 
-The time format is `HH.MMSS` with optional hundredths of a second as
-in `HH.MMSSCC`.
+Activate the Astronomical year numbering.
+
+Display negative years on the stack with a minus sign.
+This uses a uniform timeline including a year = 0.
+See: [Astronomical year numbering](https://en.wikipedia.org/wiki/Astronomical_year_numbering)
+
+**Access:** 
+
+
+## DateSlash
+## DateDash
+## DateDot
+## DateSpace
+
+Activate the corresponding date separator.
+
+**Access:** trough the `Date separator` entry in the `Status bar` menu.
+
+
+## DateSeparatorCommand
+
+Return the most recent used or default date separator command.
+There is no DateSeparator variable to Recall.
+
+**Access:**
 
 
 ## SetDate
@@ -194,31 +334,44 @@ This displays on the stack as `HH:MM:SS`.
 
 ## JulianDayNumber
 
-Return the Julian day number for the given date and time.
+Convert a date to its Julian day number.
 
 For dates the Gregorian calendar is assumed.
 The Gregorian calendar jumps from 1582-10-04 to 1582-10-15.
-This command ignores that gap, so the Julian day number given by this command for dates on or before 1582-10-14 may deviate from other converters like
-[The NASA Julian Date/Time Converter](https://ssd.jpl.nasa.gov/tools/jdc).
+This command ignores that gap.
+The Julian day number given by this command for dates on or before 1582-10-14 may deviate from other converters
+like [The NASA Julian Date/Time Converter](https://ssd.jpl.nasa.gov/tools/jdc).
+
+The conversion uses software from the
+[Standards of Fundamental Astronomy of the International Astronomical Union](https://www.iausofa.org/).
+
+**Access:** [Date](#DatesMenu) menu.
 
 To compute the Julian Day Number for the first day of the millenium:
 ```rpl
 20000101 JDN
-@ Expecting 2 451 545
+@ Expecting 2 451 544 ¹/₂
 ```
 
-## Datefromjuliandaynumber
+## DateFromJulianDayNumber
 
-Return the date for a given Julian day number.
+Convert a Julian day number to its date.
 
 This command converts a Julian day number to a date in the Gregorian calendar.
 It is the opposite of the `JDN` command.
 
+The conversion uses software from the
+[Standards of Fundamental Astronomy of the International Astronomical Union](https://www.iausofa.org/).
+
+**Access:** [Date](#DatesMenu) menu.
+
 ```rpl
 2451545 JDN→
-@ Expecting Sat 1/Jan/2000
+@ Expecting Sat 1/Jan/2000, 12:00:00
 ```
 
+# Alarms and System Commands
+
 ## ACK
 Acknowledge oldest alarm (dismiss)
 
diff --git a/sim/db48x.pro b/sim/db48x.pro
index 7ef394655..4910104cf 100644
--- a/sim/db48x.pro
+++ b/sim/db48x.pro
@@ -92,7 +92,9 @@ SOURCES +=                                      \
         ../src/unit.cc                          \
         ../src/user_interface.cc                \
         ../src/util.cc                          \
-        ../src/variables.cc
+        ../src/variables.cc                     \
+        ../src/jd2cal.c                         \
+        ../src/cal2jd.c
 
 HEADERS +=                                      \
 	sim-window.h                            \
diff --git a/src/cal2jd.c b/src/cal2jd.c
new file mode 100644
index 000000000..a0411807f
--- /dev/null
+++ b/src/cal2jd.c
@@ -0,0 +1,192 @@
+#include "sofa.h"
+#include "sofam.h"
+
+int iauCal2jd(int iy, int im, int id, double *djm0, double *djm)
+/*
+**  - - - - - - - - - -
+**   i a u C a l 2 j d
+**  - - - - - - - - - -
+**
+**  Gregorian Calendar to Julian Date.
+**
+**  This function is part of the International Astronomical Union's
+**  SOFA (Standards of Fundamental Astronomy) software collection.
+**
+**  Status:  support function.
+**
+**  Given:
+**     iy,im,id  int     year, month, day in Gregorian calendar (Note 1)
+**
+**  Returned:
+**     djm0      double  MJD zero-point: always 2400000.5
+**     djm       double  Modified Julian Date for 0 hrs
+**
+**  Returned (function value):
+**               int     status:
+**                           0 = OK
+**                          -1 = bad year   (Note 3: JD not computed)
+**                          -2 = bad month  (JD not computed)
+**                          -3 = bad day    (JD computed)
+**
+**  Notes:
+**
+**  1) The algorithm used is valid from -4800 March 1, but this
+**     implementation rejects dates before -4799 January 1.
+**
+**  2) The Julian Date is returned in two pieces, in the usual SOFA
+**     manner, which is designed to preserve time resolution.  The
+**     Julian Date is available as a single number by adding djm0 and
+**     djm.
+**
+**  3) In early eras the conversion is from the "Proleptic Gregorian
+**     Calendar";  no account is taken of the date(s) of adoption of
+**     the Gregorian Calendar, nor is the AD/BC numbering convention
+**     observed.
+**
+**  Reference:
+**
+**     Explanatory Supplement to the Astronomical Almanac,
+**     P. Kenneth Seidelmann (ed), University Science Books (1992),
+**     Section 12.92 (p604).
+**
+**  This revision:  2021 May 11
+**
+**  SOFA release 2023-10-11
+**
+**  Copyright (C) 2023 IAU SOFA Board.  See notes at end.
+*/
+{
+   int j, ly, my;
+   long iypmy;
+
+/* Earliest year allowed (4800BC) */
+   const int IYMIN = -4799;
+
+/* Month lengths in days */
+   static const int mtab[]
+                     = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
+
+
+/* Preset status. */
+   j = 0;
+
+/* Validate year and month. */
+   if (iy < IYMIN) return -1;
+   if (im < 1 || im > 12) return -2;
+
+/* If February in a leap year, 1, otherwise 0. */
+   ly = ((im == 2) && !(iy%4) && (iy%100 || !(iy%400)));
+
+/* Validate day, taking into account leap years. */
+   if ( (id < 1) || (id > (mtab[im-1] + ly))) j = -3;
+
+/* Return result. */
+   my = (im - 14) / 12;
+   iypmy = (long) (iy + my);
+   *djm0 = DJM0;
+   *djm = (double)((1461L * (iypmy + 4800L)) / 4L
+                 + (367L * (long) (im - 2 - 12 * my)) / 12L
+                 - (3L * ((iypmy + 4900L) / 100L)) / 4L
+                 + (long) id - 2432076L);
+
+/* Return status. */
+   return j;
+
+/* Finished. */
+
+/*----------------------------------------------------------------------
+**
+**  Copyright (C) 2023
+**  Standards of Fundamental Astronomy Board
+**  of the International Astronomical Union.
+**
+**  =====================
+**  SOFA Software License
+**  =====================
+**
+**  NOTICE TO USER:
+**
+**  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING SIX TERMS AND
+**  CONDITIONS WHICH APPLY TO ITS USE.
+**
+**  1. The Software is owned by the IAU SOFA Board ("SOFA").
+**
+**  2. Permission is granted to anyone to use the SOFA software for any
+**     purpose, including commercial applications, free of charge and
+**     without payment of royalties, subject to the conditions and
+**     restrictions listed below.
+**
+**  3. You (the user) may copy and distribute SOFA source code to others,
+**     and use and adapt its code and algorithms in your own software,
+**     on a world-wide, royalty-free basis.  That portion of your
+**     distribution that does not consist of intact and unchanged copies
+**     of SOFA source code files is a "derived work" that must comply
+**     with the following requirements:
+**
+**     a) Your work shall be marked or carry a statement that it
+**        (i) uses routines and computations derived by you from
+**        software provided by SOFA under license to you; and
+**        (ii) does not itself constitute software provided by and/or
+**        endorsed by SOFA.
+**
+**     b) The source code of your derived work must contain descriptions
+**        of how the derived work is based upon, contains and/or differs
+**        from the original SOFA software.
+**
+**     c) The names of all routines in your derived work shall not
+**        include the prefix "iau" or "sofa" or trivial modifications
+**        thereof such as changes of case.
+**
+**     d) The origin of the SOFA components of your derived work must
+**        not be misrepresented;  you must not claim that you wrote the
+**        original software, nor file a patent application for SOFA
+**        software or algorithms embedded in the SOFA software.
+**
+**     e) These requirements must be reproduced intact in any source
+**        distribution and shall apply to anyone to whom you have
+**        granted a further right to modify the source code of your
+**        derived work.
+**
+**     Note that, as originally distributed, the SOFA software is
+**     intended to be a definitive implementation of the IAU standards,
+**     and consequently third-party modifications are discouraged.  All
+**     variations, no matter how minor, must be explicitly marked as
+**     such, as explained above.
+**
+**  4. You shall not cause the SOFA software to be brought into
+**     disrepute, either by misuse, or use for inappropriate tasks, or
+**     by inappropriate modification.
+**
+**  5. The SOFA software is provided "as is" and SOFA makes no warranty
+**     as to its use or performance.   SOFA does not and cannot warrant
+**     the performance or results which the user may obtain by using the
+**     SOFA software.  SOFA makes no warranties, express or implied, as
+**     to non-infringement of third party rights, merchantability, or
+**     fitness for any particular purpose.  In no event will SOFA be
+**     liable to the user for any consequential, incidental, or special
+**     damages, including any lost profits or lost savings, even if a
+**     SOFA representative has been advised of such damages, or for any
+**     claim by any third party.
+**
+**  6. The provision of any version of the SOFA software under the terms
+**     and conditions specified herein does not imply that future
+**     versions will also be made available under the same terms and
+**     conditions.
+*
+**  In any published work or commercial product which uses the SOFA
+**  software directly, acknowledgement (see www.iausofa.org) is
+**  appreciated.
+**
+**  Correspondence concerning SOFA software should be addressed as
+**  follows:
+**
+**      By email:  sofa@ukho.gov.uk
+**      By post:   IAU SOFA Center
+**                 HM Nautical Almanac Office
+**                 UK Hydrographic Office
+**                 Admiralty Way, Taunton
+**                 Somerset, TA1 2DN
+**                 United Kingdom
+**
+**--------------------------------------------------------------------*/
+}
diff --git a/src/datetime.cc b/src/datetime.cc
index d97998d44..63a0b5a50 100644
--- a/src/datetime.cc
+++ b/src/datetime.cc
@@ -36,6 +36,8 @@
 #include "tag.h"
 #include "unit.h"
 
+#include "sofa.h"
+#include "sofam.h"
 
 // ============================================================================
 //
@@ -132,7 +134,8 @@ uint to_date(object_p dtobj, dt_t &dt, tm_t &tm, bool error)
 
     const uint days[12] = { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };
     bool bisext = m == 2 && y % 4 == 0 && (y % 100 != 0 || y % 400 == 0);
-    if (m < 1 || m > 12 || d < 1 || d > days[m-1] + bisext || y == 0)
+    if (m < 1 || m > 12 || d < 1 || d > days[m-1] + bisext ||
+	(Settings.BCECEyearNumbering() && y == 0))
     {
         if (error)
             rt.invalid_date_error();
@@ -201,9 +204,14 @@ algebraic_p julian_day_number(const dt_t &dt, const tm_t &tm)
 //   Compute julian day number for a dt_t structure
 // ----------------------------------------------------------------------------
 {
+    double jdn_value = julian_day_number(dt.day, dt.month, dt.year);
+    long jdn_day = long(jdn_value);
+    long jdn_fraction_of_day = long((jdn_value - jdn_day) * 8640000);
+
     ularge csecs = (tm.hour * 3600 + tm.min * 60 + tm.sec) * 100 + tm.csec;
-    ularge jval = julian_day_number(dt.day, dt.month, dt.year);
-    algebraic_g jdn = integer::make(jval);
+    algebraic_g jdn = integer::make(jdn_day);
+
+    csecs += jdn_fraction_of_day;
     if (csecs)
     {
         algebraic_g frac = +fraction::make(integer::make(csecs),
@@ -227,27 +235,21 @@ algebraic_p julian_day_number(algebraic_p dtobj, bool error)
 }
 
 
-ularge julian_day_number(int d, int m, int y)
+double julian_day_number(int d, int m, int y)
 // ----------------------------------------------------------------------------
 //   Compute the Julian day number given day, month and year
 // ----------------------------------------------------------------------------
-// The algorithm below follows the description in
-// https://en.wikipedia.org/wiki/Julian_day#Converting_Gregorian_calendar_date_to_Julian_day_number
-// which itself seems to be based on the "Fliegel and van Flandern" algorithms:
-// https://aa.usno.navy.mil/faq/JD_formula --> JD.
-// With a little bit of algebra this function can be shown to be equivalent.
-//
-// According to https://en.wikipedia.org/wiki/Gregorian_calendar
-// "Thursday 4 October 1582 was followed by Friday 15 October 1582"
-// This function doesn't do the 10 day advance.
-{
-    int rm = (m-14)/12;
-    ularge jdn = ((1461 * (y + 4800 + rm)) / 4
-                  + (367 * (m - 2 - 12 * rm)) / 12
-                  - (3 * ((y + 4900 + rm) / 100)) / 4
-                  + d
-                  - 32075);
-    return jdn;
+{
+    int result;
+    double djm0, djm;
+
+    result = iauCal2jd(y, m, d, &djm0, &djm);
+    if (result < 0)
+    {
+        rt.invalid_date_error();
+        return 0.0;
+    }
+    return (djm0 + djm);
 }
 
 
@@ -255,63 +257,65 @@ algebraic_p date_from_julian_day(object_p jdn, bool error)
 // ----------------------------------------------------------------------------
 //   Create a date from a Julian day object
 // ----------------------------------------------------------------------------
-// See comments in julian_day_number
 {
     if (!jdn)
         return nullptr;
 
     if (algebraic_g jval = jdn->as_real())
     {
-        large jdn = jval->as_int64(0, error);
+        double dj1, dj2 = 0.0;
+        int year, month, day;
+        double jdn_fraction_of_day;
+        int result;
 
-        enum
+        // How on earth do I get the real or double value out of jval?
         {
-            y = 4716,
-            j = 1401,
-            m = 2,
-            n = 12,
-            r = 4,
-            p = 1461,
-            v = 3,
-            u = 5,
-            s = 153,
-            w = 2,
-            B = 274277,
-            C = -38
-        };
-
-        large f = jdn + j + (((4 * jdn + B) / 146097) * 3) / 4 + C;
-        large e = r * f + v;
-        large g = (e % p) / r;
-        large h = u * g + w;
-        uint day = (h % s) / u + 1;
-        uint month = (h / s + m ) % n + 1;
-        int year = e / p - y + (n + m - month) / n;
-
-	bool negativeYear = year < 0;
-	if (negativeYear)
+	    // Format hhmmsscc, thus cc is hundreds of a second.
+	    const long scale_to_cc = 100 * 100 * 100 * 100;
+            algebraic_g scale = integer::make(scale_to_cc);
+            algebraic_g j = jval * scale;
+            dj1 = double(j->as_uint64(0, false)) / scale_to_cc;
+        }
+
+        result = iauJd2cal(dj1, dj2, &year, &month, &day, &jdn_fraction_of_day);
+        if (result < 0)
+	{
+            if (error)
+	        rt.invalid_date_error();
+            return nullptr;
+	}
+
+        bool negative_year = year < 0;
+        if (negative_year)
             year = -year;
 
-        ularge dval = year * 10000 + month * 100 + day;
+        algebraic_g date_part = nullptr;
+        {
+          int date = (((year * 100) + month) * 100) + day;
+          date_part = integer::make(date);
+        }
 
-        algebraic_g date = integer::make(dval);
-        algebraic_g fp = integer::make(1);
-        fp = jval % fp;
-        if (!fp->is_zero())
+        algebraic_g time_part = nullptr;
         {
-            algebraic_g factor = integer::make(86400);
-            fp = fp * factor;
-            ularge hval = fp->as_uint64(0, false);
-            uint hour = hval / 3600;
-            uint min = (hval / 60) % 60;
-            uint sec = hval % 60;
-            hval = hour * 10000 + min * 100 + sec;
-            fp = +fraction::make(integer::make(hval),
-                                 integer::make(1000000));
-            date = date + fp;
+            const long scale_to_cc = 24 * 60 * 60 * 100;
+            const long scale_from_cc = 100 * 100 * 100 * 100;
+            int d = round(jdn_fraction_of_day * scale_to_cc);
+            int cc = d % 100;
+            d = d / 100;
+            int seconds = d % 60;
+            d = d / 60;
+            int minutes = d % 60;
+            d = d / 60;
+            int hours = d % 24;
+
+            int time = ((hours * 100 + minutes) * 100 + seconds) * 100 + cc;
+            time_part = +fraction::make(integer::make(time),
+                                        +integer::make(scale_from_cc));
         }
 
-	if (negativeYear)
+        algebraic_g date = date_part + time_part;
+
+        if (negative_year)
             date = -date;
 
         date = unit::make(date, +symbol::make("date"));
@@ -498,9 +502,6 @@ COMMAND_BODY(ChronoTime)
 }
 
 
-
-
-
 static bool setTime(object_p tobj)
 // ----------------------------------------------------------------------------
 //   Set the system date based on input
@@ -529,7 +530,6 @@ COMMAND_BODY(SetTime)
 }
 
 
-
 void render_time(renderer &r, algebraic_g &value,
                  cstring hrs, cstring min, cstring sec,
                  uint base, bool ampm)
@@ -609,8 +609,9 @@ size_t render_date(renderer &r, algebraic_g date)
 {
     if (!date || !date->is_real())
         return 0;
-    bool neg = date->is_negative();
-    if (neg)
+    
+    bool negative_date = date->is_negative();
+    if (negative_date)
         date = -date;
 
     algebraic_g factor = integer::make(100);
@@ -623,11 +624,14 @@ size_t render_date(renderer &r, algebraic_g date)
     date = date / factor;
     uint month = date->as_uint32(0, false) % 100;
     date = date / factor;
-    uint year = date->as_uint32(0, false);
-    if (year == 0)
+    int year = date->as_int32(0, false);
+
+    if (Settings.ShowDayOfWeek())
     {
-        rt.invalid_date_error();
-        return 0;
+        int syear = negative_date ? -int(year) : year;
+        ularge jdn = julian_day_number(day, month, syear);
+        uint dow = (jdn + 1) % 7;
+        r.printf("%s ", get_wday_shortcut(dow));
     }
 
     char mname[4];
@@ -636,19 +640,20 @@ size_t render_date(renderer &r, algebraic_g date)
     else
         snprintf(mname, 4, "%u", month);
 
+    bool zero_year = year == 0;
+    if (zero_year && !Settings.AstronomicalYearNumbering())
+        return 0;
+
+    if ((negative_date || zero_year) && Settings.BCECEyearNumbering())
+        year = year + 1;
+    if (negative_date && Settings.AstronomicalYearNumbering())
+        year = - year;
+    
     char ytext[16];
     if (Settings.TwoDigitYear())
-        snprintf(ytext, sizeof(ytext), "%02u", year % 100);
+        snprintf(ytext, sizeof(ytext), "%02i", year % 100);
     else
-        snprintf(ytext, sizeof(ytext), "%u", year);
-
-    if (Settings.ShowDayOfWeek())
-    {
-        int syear = neg ? -int(year) : year;
-        ularge jdn = julian_day_number(day, month, syear);
-        uint dow = jdn % 7;
-        r.printf("%s ", get_wday_shortcut(dow));
-    }
+        snprintf(ytext, sizeof(ytext), "%i", year);
 
     char sep   = Settings.DateSeparator();
     uint index = 2 * Settings.YearFirst() + Settings.MonthBeforeDay();
@@ -659,8 +664,13 @@ size_t render_date(renderer &r, algebraic_g date)
     case 2: r.printf("%s%c%u%c%s", ytext, sep, day,   sep, mname); break;
     case 3: r.printf("%s%c%s%c%u", ytext, sep, mname, sep, day);   break;
     }
-    if (neg)
-        r.printf(" BC");
+    if (Settings.BCECEyearNumbering())
+    {
+        r.printf(" ");
+        if (negative_date || zero_year)
+            r.printf("B");
+        r.printf("CE");
+    }
 
     if (time && !time->is_zero())
     {
@@ -673,7 +683,6 @@ size_t render_date(renderer &r, algebraic_g date)
 }
 
 
-
 // ============================================================================
 //
 //   HMS and DMS commands
diff --git a/src/datetime.h b/src/datetime.h
index dd97dc9a4..36a796f83 100644
--- a/src/datetime.h
+++ b/src/datetime.h
@@ -42,7 +42,7 @@ algebraic_p to_days(object_p days, bool error = true);
 
 // Convert date value to Julian day number, or 0 if fails
 algebraic_p julian_day_number(algebraic_p date, bool error = true);
-ularge julian_day_number(int d, int m, int y);
+double julian_day_number(int d, int m, int y);
 
 // Convert Julian day number to date
 algebraic_p date_from_julian_day(object_p jdn, bool error = true);
diff --git a/src/ids.tbl b/src/ids.tbl
index fb5c1529d..f58d52d41 100644
--- a/src/ids.tbl
+++ b/src/ids.tbl
@@ -1041,6 +1041,7 @@ FLAG(NumericalIntegration,      SymbolicIntegration)
 FLAG(TVMPayAtBeginningOfPeriod, TVMPayAtEndOfPeriod)
 FLAG(TruthLogicForIntegers,     BitwiseLogicForIntegers)
 FLAG(LaxArrayResizing,          StrictArrayResizing)
+FLAG(BCECEyearNumbering,        AstronomicalYearNumbering)
 
 ALIAS(HardwareFloatingPoint,    "HFP")
 ALIAS(HardwareFloatingPoint,    "HardFP")
@@ -1342,7 +1343,7 @@ MENU(PrintingMenu)
 MENU(IOMenu)
 MENU(FilesMenu)
 MENU(TimeMenu)
-MENU(DateMenu)
+MENU(DatesMenu)
 MENU(AlarmMenu)
 MENU(PlotMenu)
 MENU(GraphicsMenu)
diff --git a/src/jd2cal.c b/src/jd2cal.c
new file mode 100644
index 000000000..036c3cd89
--- /dev/null
+++ b/src/jd2cal.c
@@ -0,0 +1,251 @@
+#include "sofa.h"
+#include "sofam.h"
+#include <float.h>
+
+int iauJd2cal(double dj1, double dj2,
+              int *iy, int *im, int *id, double *fd)
+/*
+**  - - - - - - - - - -
+**   i a u J d 2 c a l
+**  - - - - - - - - - -
+**
+**  Julian Date to Gregorian year, month, day, and fraction of a day.
+**
+**  This function is part of the International Astronomical Union's
+**  SOFA (Standards of Fundamental Astronomy) software collection.
+**
+**  Status:  support function.
+**
+**  Given:
+**     dj1,dj2   double   Julian Date (Notes 1, 2)
+**
+**  Returned (arguments):
+**     iy        int      year
+**     im        int      month
+**     id        int      day
+**     fd        double   fraction of day
+**
+**  Returned (function value):
+**               int      status:
+**                           0 = OK
+**                          -1 = unacceptable date (Note 1)
+**
+**  Notes:
+**
+**  1) The earliest valid date is -68569.5 (-4900 March 1).  The
+**     largest value accepted is 1e9.
+**
+**  2) The Julian Date is apportioned in any convenient way between
+**     the arguments dj1 and dj2.  For example, JD=2450123.7 could
+**     be expressed in any of these ways, among others:
+**
+**            dj1             dj2
+**
+**         2450123.7           0.0       (JD method)
+**         2451545.0       -1421.3       (J2000 method)
+**         2400000.5       50123.2       (MJD method)
+**         2450123.5           0.2       (date & time method)
+**
+**     Separating integer and fraction uses the "compensated summation"
+**     algorithm of Kahan-Neumaier to preserve as much precision as
+**     possible irrespective of the jd1+jd2 apportionment.
+**
+**  3) In early eras the conversion is from the "proleptic Gregorian
+**     calendar";  no account is taken of the date(s) of adoption of
+**     the Gregorian calendar, nor is the AD/BC numbering convention
+**     observed.
+**
+**  References:
+**
+**     Explanatory Supplement to the Astronomical Almanac,
+**     P. Kenneth Seidelmann (ed), University Science Books (1992),
+**     Section 12.92 (p604).
+**
+**     Klein, A., A Generalized Kahan-Babuska-Summation-Algorithm.
+**     Computing, 76, 279-293 (2006), Section 3.
+**
+**  This revision:  2021 May 11
+**
+**  SOFA release 2023-10-11
+**
+**  Copyright (C) 2023 IAU SOFA Board.  See notes at end.
+*/
+{
+/* Minimum and maximum allowed JD */
+   const double DJMIN = -68569.5;
+   const double DJMAX = 1e9;
+
+   long jd, i, l, n, k;
+   double dj, f1, f2, d, s, cs, v[2], x, t, f;
+
+
+/* Verify date is acceptable. */
+   dj = dj1 + dj2;
+   if (dj < DJMIN || dj > DJMAX) return -1;
+
+/* Separate day and fraction (where -0.5 <= fraction < 0.5). */
+   d = dnint(dj1);
+   f1 = dj1 - d;
+   jd = (long) d;
+   d = dnint(dj2);
+   f2 = dj2 - d;
+   jd += (long) d;
+
+/* Compute f1+f2+0.5 using compensated summation (Klein 2006). */
+   s = 0.5;
+   cs = 0.0;
+   v[0] = f1;
+   v[1] = f2;
+   for ( i = 0; i < 2; i++ ) {
+      x = v[i];
+      t = s + x;
+      cs += fabs(s) >= fabs(x) ? (s-t) + x : (x-t) + s;
+      s = t;
+      if ( s >= 1.0 ) {
+         jd++;
+         s -= 1.0;
+      }
+   }
+   f = s + cs;
+   cs = f - s;
+
+/* Deal with negative f. */
+   if ( f < 0.0 ) {
+
+   /* Compensated summation: assume that |s| <= 1.0. */
+      f = s + 1.0;
+      cs += (1.0-f) + s;
+      s  = f;
+      f = s + cs;
+      cs = f - s;
+      jd--;
+   }
+
+/* Deal with f that is 1.0 or more (when rounded to double). */
+   if ( (f-1.0) >= -DBL_EPSILON/4.0 ) {
+
+   /* Compensated summation: assume that |s| <= 1.0. */
+      t = s - 1.0;
+      cs += (s-t) - 1.0;
+      s = t;
+      f = s + cs;
+      if ( -DBL_EPSILON/2.0 < f ) {
+         jd++;
+         f = gmax(f, 0.0);
+      }
+   }
+
+/* Express day in Gregorian calendar. */
+   l = jd + 68569L;
+   n = (4L * l) / 146097L;
+   l -= (146097L * n + 3L) / 4L;
+   i = (4000L * (l + 1L)) / 1461001L;
+   l -= (1461L * i) / 4L - 31L;
+   k = (80L * l) / 2447L;
+   *id = (int) (l - (2447L * k) / 80L);
+   l = k / 11L;
+   *im = (int) (k + 2L - 12L * l);
+   *iy = (int) (100L * (n - 49L) + i + l);
+   *fd = f;
+
+/* Success. */
+   return 0;
+
+/* Finished. */
+
+/*----------------------------------------------------------------------
+**
+**  Copyright (C) 2023
+**  Standards of Fundamental Astronomy Board
+**  of the International Astronomical Union.
+**
+**  =====================
+**  SOFA Software License
+**  =====================
+**
+**  NOTICE TO USER:
+**
+**  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING SIX TERMS AND
+**  CONDITIONS WHICH APPLY TO ITS USE.
+**
+**  1. The Software is owned by the IAU SOFA Board ("SOFA").
+**
+**  2. Permission is granted to anyone to use the SOFA software for any
+**     purpose, including commercial applications, free of charge and
+**     without payment of royalties, subject to the conditions and
+**     restrictions listed below.
+**
+**  3. You (the user) may copy and distribute SOFA source code to others,
+**     and use and adapt its code and algorithms in your own software,
+**     on a world-wide, royalty-free basis.  That portion of your
+**     distribution that does not consist of intact and unchanged copies
+**     of SOFA source code files is a "derived work" that must comply
+**     with the following requirements:
+**
+**     a) Your work shall be marked or carry a statement that it
+**        (i) uses routines and computations derived by you from
+**        software provided by SOFA under license to you; and
+**        (ii) does not itself constitute software provided by and/or
+**        endorsed by SOFA.
+**
+**     b) The source code of your derived work must contain descriptions
+**        of how the derived work is based upon, contains and/or differs
+**        from the original SOFA software.
+**
+**     c) The names of all routines in your derived work shall not
+**        include the prefix "iau" or "sofa" or trivial modifications
+**        thereof such as changes of case.
+**
+**     d) The origin of the SOFA components of your derived work must
+**        not be misrepresented;  you must not claim that you wrote the
+**        original software, nor file a patent application for SOFA
+**        software or algorithms embedded in the SOFA software.
+**
+**     e) These requirements must be reproduced intact in any source
+**        distribution and shall apply to anyone to whom you have
+**        granted a further right to modify the source code of your
+**        derived work.
+**
+**     Note that, as originally distributed, the SOFA software is
+**     intended to be a definitive implementation of the IAU standards,
+**     and consequently third-party modifications are discouraged.  All
+**     variations, no matter how minor, must be explicitly marked as
+**     such, as explained above.
+**
+**  4. You shall not cause the SOFA software to be brought into
+**     disrepute, either by misuse, or use for inappropriate tasks, or
+**     by inappropriate modification.
+**
+**  5. The SOFA software is provided "as is" and SOFA makes no warranty
+**     as to its use or performance.   SOFA does not and cannot warrant
+**     the performance or results which the user may obtain by using the
+**     SOFA software.  SOFA makes no warranties, express or implied, as
+**     to non-infringement of third party rights, merchantability, or
+**     fitness for any particular purpose.  In no event will SOFA be
+**     liable to the user for any consequential, incidental, or special
+**     damages, including any lost profits or lost savings, even if a
+**     SOFA representative has been advised of such damages, or for any
+**     claim by any third party.
+**
+**  6. The provision of any version of the SOFA software under the terms
+**     and conditions specified herein does not imply that future
+**     versions will also be made available under the same terms and
+**     conditions.
+*
+**  In any published work or commercial product which uses the SOFA
+**  software directly, acknowledgement (see www.iausofa.org) is
+**  appreciated.
+**
+**  Correspondence concerning SOFA software should be addressed as
+**  follows:
+**
+**      By email:  sofa@ukho.gov.uk
+**      By post:   IAU SOFA Center
+**                 HM Nautical Almanac Office
+**                 UK Hydrographic Office
+**                 Admiralty Way, Taunton
+**                 Somerset, TA1 2DN
+**                 United Kingdom
+**
+**--------------------------------------------------------------------*/
+}
diff --git a/src/menu.cc b/src/menu.cc
index 8cda94391..2b6e91435 100644
--- a/src/menu.cc
+++ b/src/menu.cc
@@ -198,7 +198,7 @@ static object::id unit_menu(unit_p u)
                 result = object::ID_TimeMenu;
             else if (sym->matches("date") || sym->matches("d") ||
                      sym->matches("yr"))
-                result = object::ID_DateMenu;
+                result = object::ID_DatesMenu;
         }
     }
     return result;
@@ -1773,11 +1773,11 @@ MENU(TimeMenu,
      "→Time",   ID_SetTime,
      "→Date",   ID_SetDate,
      "ClkAdj",  ID_Unimplemented,
-     "Dates",   ID_DateMenu,
+     "Dates",   ID_DatesMenu,
      "Alarms",  ID_AlarmMenu);
 
 
-MENU(DateMenu,
+MENU(DatesMenu,
 // ----------------------------------------------------------------------------
 //   Date operations
 // ----------------------------------------------------------------------------
@@ -1809,7 +1809,7 @@ MENU(AlarmMenu,
 
      "AckAll",  ID_Unimplemented,
      "Time",    ID_TimeMenu,
-     "Date",    ID_DateMenu);
+     "Dates",   ID_DatesMenu);
 
 
 MENU(TextMenu,
diff --git a/src/sofa.h b/src/sofa.h
new file mode 100644
index 000000000..b397769d3
--- /dev/null
+++ b/src/sofa.h
@@ -0,0 +1,627 @@
+#ifndef SOFAHDEF
+#define SOFAHDEF
+
+/*
+**  - - - - - - -
+**   s o f a . h
+**  - - - - - - -
+**
+**  Prototype function declarations for SOFA library.
+**
+**  This file is part of the International Astronomical Union's
+**  SOFA (Standards of Fundamental Astronomy) software collection.
+**
+**  This revision:   2023 April 16
+**
+**  SOFA release 2023-10-11
+**
+**  Copyright (C) 2023 IAU SOFA Board.  See notes at end.
+*/
+
+#include "math.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* Star-independent astrometry parameters */
+typedef struct {
+   double pmt;        /* PM time interval (SSB, Julian years) */
+   double eb[3];      /* SSB to observer (vector, au) */
+   double eh[3];      /* Sun to observer (unit vector) */
+   double em;         /* distance from Sun to observer (au) */
+   double v[3];       /* barycentric observer velocity (vector, c) */
+   double bm1;        /* sqrt(1-|v|^2): reciprocal of Lorenz factor */
+   double bpn[3][3];  /* bias-precession-nutation matrix */
+   double along;      /* longitude + s' + dERA(DUT) (radians) */
+   double phi;        /* geodetic latitude (radians) */
+   double xpl;        /* polar motion xp wrt local meridian (radians) */
+   double ypl;        /* polar motion yp wrt local meridian (radians) */
+   double sphi;       /* sine of geodetic latitude */
+   double cphi;       /* cosine of geodetic latitude */
+   double diurab;     /* magnitude of diurnal aberration vector */
+   double eral;       /* "local" Earth rotation angle (radians) */
+   double refa;       /* refraction constant A (radians) */
+   double refb;       /* refraction constant B (radians) */
+} iauASTROM;
+/* (Vectors eb, eh, em and v are all with respect to BCRS axes.) */
+
+/* Body parameters for light deflection */
+typedef struct {
+   double bm;         /* mass of the body (solar masses) */
+   double dl;         /* deflection limiter (radians^2/2) */
+   double pv[2][3];   /* barycentric PV of the body (au, au/day) */
+} iauLDBODY;
+
+/* Astronomy/Calendars */
+int iauCal2jd(int iy, int im, int id, double *djm0, double *djm);
+double iauEpb(double dj1, double dj2);
+void iauEpb2jd(double epb, double *djm0, double *djm);
+double iauEpj(double dj1, double dj2);
+void iauEpj2jd(double epj, double *djm0, double *djm);
+int iauJd2cal(double dj1, double dj2,
+                     int *iy, int *im, int *id, double *fd);
+int iauJdcalf(int ndp, double dj1, double dj2, int iymdf[4]);
+
+/* Astronomy/Astrometry */
+void iauAb(double pnat[3], double v[3], double s, double bm1,
+           double ppr[3]);
+void iauApcg(double date1, double date2,
+             double ebpv[2][3], double ehp[3],
+             iauASTROM *astrom);
+void iauApcg13(double date1, double date2, iauASTROM *astrom);
+void iauApci(double date1, double date2,
+             double ebpv[2][3], double ehp[3],
+             double x, double y, double s,
+             iauASTROM *astrom);
+void iauApci13(double date1, double date2,
+               iauASTROM *astrom, double *eo);
+void iauApco(double date1, double date2,
+             double ebpv[2][3], double ehp[3],
+             double x, double y, double s, double theta,
+             double elong, double phi, double hm,
+             double xp, double yp, double sp,
+             double refa, double refb,
+             iauASTROM *astrom);
+int iauApco13(double utc1, double utc2, double dut1,
+              double elong, double phi, double hm, double xp, double yp,
+              double phpa, double tc, double rh, double wl,
+              iauASTROM *astrom, double *eo);
+void iauApcs(double date1, double date2, double pv[2][3],
+             double ebpv[2][3], double ehp[3],
+             iauASTROM *astrom);
+void iauApcs13(double date1, double date2, double pv[2][3],
+               iauASTROM *astrom);
+void iauAper(double theta, iauASTROM *astrom);
+void iauAper13(double ut11, double ut12, iauASTROM *astrom);
+void iauApio(double sp, double theta,
+             double elong, double phi, double hm, double xp, double yp,
+             double refa, double refb,
+             iauASTROM *astrom);
+int iauApio13(double utc1, double utc2, double dut1,
+              double elong, double phi, double hm, double xp, double yp,
+              double phpa, double tc, double rh, double wl,
+              iauASTROM *astrom);
+void iauAtcc13(double rc, double dc,
+               double pr, double pd, double px, double rv,
+               double date1, double date2,
+               double *ra, double *da);
+void iauAtccq(double rc, double dc,
+              double pr, double pd, double px, double rv,
+              iauASTROM *astrom, double *ra, double *da);
+void iauAtci13(double rc, double dc,
+               double pr, double pd, double px, double rv,
+               double date1, double date2,
+               double *ri, double *di, double *eo);
+void iauAtciq(double rc, double dc, double pr, double pd,
+              double px, double rv, iauASTROM *astrom,
+              double *ri, double *di);
+void iauAtciqn(double rc, double dc, double pr, double pd,
+               double px, double rv, iauASTROM *astrom,
+               int n, iauLDBODY b[], double *ri, double *di);
+void iauAtciqz(double rc, double dc, iauASTROM *astrom,
+               double *ri, double *di);
+int iauAtco13(double rc, double dc,
+              double pr, double pd, double px, double rv,
+              double utc1, double utc2, double dut1,
+              double elong, double phi, double hm, double xp, double yp,
+              double phpa, double tc, double rh, double wl,
+              double *aob, double *zob, double *hob,
+              double *dob, double *rob, double *eo);
+void iauAtic13(double ri, double di,
+               double date1, double date2,
+               double *rc, double *dc, double *eo);
+void iauAticq(double ri, double di, iauASTROM *astrom,
+              double *rc, double *dc);
+void iauAticqn(double ri, double di, iauASTROM *astrom,
+               int n, iauLDBODY b[], double *rc, double *dc);
+int iauAtio13(double ri, double di,
+              double utc1, double utc2, double dut1,
+              double elong, double phi, double hm, double xp, double yp,
+              double phpa, double tc, double rh, double wl,
+              double *aob, double *zob, double *hob,
+              double *dob, double *rob);
+void iauAtioq(double ri, double di, iauASTROM *astrom,
+              double *aob, double *zob,
+              double *hob, double *dob, double *rob);
+int iauAtoc13(const char *type, double ob1, double ob2,
+              double utc1, double utc2, double dut1,
+              double elong, double phi, double hm, double xp, double yp,
+              double phpa, double tc, double rh, double wl,
+              double *rc, double *dc);
+int iauAtoi13(const char *type, double ob1, double ob2,
+              double utc1, double utc2, double dut1,
+              double elong, double phi, double hm, double xp, double yp,
+              double phpa, double tc, double rh, double wl,
+              double *ri, double *di);
+void iauAtoiq(const char *type,
+              double ob1, double ob2, iauASTROM *astrom,
+              double *ri, double *di);
+void iauLd(double bm, double p[3], double q[3], double e[3],
+           double em, double dlim, double p1[3]);
+void iauLdn(int n, iauLDBODY b[], double ob[3], double sc[3],
+            double sn[3]);
+void iauLdsun(double p[3], double e[3], double em, double p1[3]);
+void iauPmpx(double rc, double dc, double pr, double pd,
+             double px, double rv, double pmt, double pob[3],
+             double pco[3]);
+int iauPmsafe(double ra1, double dec1, double pmr1, double pmd1,
+              double px1, double rv1,
+              double ep1a, double ep1b, double ep2a, double ep2b,
+              double *ra2, double *dec2, double *pmr2, double *pmd2,
+              double *px2, double *rv2);
+void iauPvtob(double elong, double phi, double height, double xp,
+              double yp, double sp, double theta, double pv[2][3]);
+void iauRefco(double phpa, double tc, double rh, double wl,
+              double *refa, double *refb);
+
+/* Astronomy/Ephemerides */
+int iauEpv00(double date1, double date2,
+             double pvh[2][3], double pvb[2][3]);
+void iauMoon98(double date1, double date2, double pv[2][3]);
+int iauPlan94(double date1, double date2, int np, double pv[2][3]);
+
+/* Astronomy/FundamentalArgs */
+double iauFad03(double t);
+double iauFae03(double t);
+double iauFaf03(double t);
+double iauFaju03(double t);
+double iauFal03(double t);
+double iauFalp03(double t);
+double iauFama03(double t);
+double iauFame03(double t);
+double iauFane03(double t);
+double iauFaom03(double t);
+double iauFapa03(double t);
+double iauFasa03(double t);
+double iauFaur03(double t);
+double iauFave03(double t);
+
+/* Astronomy/PrecNutPolar */
+void iauBi00(double *dpsibi, double *depsbi, double *dra);
+void iauBp00(double date1, double date2,
+             double rb[3][3], double rp[3][3], double rbp[3][3]);
+void iauBp06(double date1, double date2,
+             double rb[3][3], double rp[3][3], double rbp[3][3]);
+void iauBpn2xy(double rbpn[3][3], double *x, double *y);
+void iauC2i00a(double date1, double date2, double rc2i[3][3]);
+void iauC2i00b(double date1, double date2, double rc2i[3][3]);
+void iauC2i06a(double date1, double date2, double rc2i[3][3]);
+void iauC2ibpn(double date1, double date2, double rbpn[3][3],
+               double rc2i[3][3]);
+void iauC2ixy(double date1, double date2, double x, double y,
+              double rc2i[3][3]);
+void iauC2ixys(double x, double y, double s, double rc2i[3][3]);
+void iauC2t00a(double tta, double ttb, double uta, double utb,
+               double xp, double yp, double rc2t[3][3]);
+void iauC2t00b(double tta, double ttb, double uta, double utb,
+               double xp, double yp, double rc2t[3][3]);
+void iauC2t06a(double tta, double ttb, double uta, double utb,
+               double xp, double yp, double rc2t[3][3]);
+void iauC2tcio(double rc2i[3][3], double era, double rpom[3][3],
+               double rc2t[3][3]);
+void iauC2teqx(double rbpn[3][3], double gst, double rpom[3][3],
+               double rc2t[3][3]);
+void iauC2tpe(double tta, double ttb, double uta, double utb,
+              double dpsi, double deps, double xp, double yp,
+              double rc2t[3][3]);
+void iauC2txy(double tta, double ttb, double uta, double utb,
+              double x, double y, double xp, double yp,
+              double rc2t[3][3]);
+double iauEo06a(double date1, double date2);
+double iauEors(double rnpb[3][3], double s);
+void iauFw2m(double gamb, double phib, double psi, double eps,
+             double r[3][3]);
+void iauFw2xy(double gamb, double phib, double psi, double eps,
+              double *x, double *y);
+void iauLtp(double epj, double rp[3][3]);
+void iauLtpb(double epj, double rpb[3][3]);
+void iauLtpecl(double epj, double vec[3]);
+void iauLtpequ(double epj, double veq[3]);
+void iauNum00a(double date1, double date2, double rmatn[3][3]);
+void iauNum00b(double date1, double date2, double rmatn[3][3]);
+void iauNum06a(double date1, double date2, double rmatn[3][3]);
+void iauNumat(double epsa, double dpsi, double deps, double rmatn[3][3]);
+void iauNut00a(double date1, double date2, double *dpsi, double *deps);
+void iauNut00b(double date1, double date2, double *dpsi, double *deps);
+void iauNut06a(double date1, double date2, double *dpsi, double *deps);
+void iauNut80(double date1, double date2, double *dpsi, double *deps);
+void iauNutm80(double date1, double date2, double rmatn[3][3]);
+double iauObl06(double date1, double date2);
+double iauObl80(double date1, double date2);
+void iauP06e(double date1, double date2,
+             double *eps0, double *psia, double *oma, double *bpa,
+             double *bqa, double *pia, double *bpia,
+             double *epsa, double *chia, double *za, double *zetaa,
+             double *thetaa, double *pa,
+             double *gam, double *phi, double *psi);
+void iauPb06(double date1, double date2,
+             double *bzeta, double *bz, double *btheta);
+void iauPfw06(double date1, double date2,
+              double *gamb, double *phib, double *psib, double *epsa);
+void iauPmat00(double date1, double date2, double rbp[3][3]);
+void iauPmat06(double date1, double date2, double rbp[3][3]);
+void iauPmat76(double date1, double date2, double rmatp[3][3]);
+void iauPn00(double date1, double date2, double dpsi, double deps,
+             double *epsa,
+             double rb[3][3], double rp[3][3], double rbp[3][3],
+             double rn[3][3], double rbpn[3][3]);
+void iauPn00a(double date1, double date2,
+              double *dpsi, double *deps, double *epsa,
+              double rb[3][3], double rp[3][3], double rbp[3][3],
+              double rn[3][3], double rbpn[3][3]);
+void iauPn00b(double date1, double date2,
+              double *dpsi, double *deps, double *epsa,
+              double rb[3][3], double rp[3][3], double rbp[3][3],
+              double rn[3][3], double rbpn[3][3]);
+void iauPn06(double date1, double date2, double dpsi, double deps,
+             double *epsa,
+             double rb[3][3], double rp[3][3], double rbp[3][3],
+             double rn[3][3], double rbpn[3][3]);
+void iauPn06a(double date1, double date2,
+              double *dpsi, double *deps, double *epsa,
+              double rb[3][3], double rp[3][3], double rbp[3][3],
+              double rn[3][3], double rbpn[3][3]);
+void iauPnm00a(double date1, double date2, double rbpn[3][3]);
+void iauPnm00b(double date1, double date2, double rbpn[3][3]);
+void iauPnm06a(double date1, double date2, double rnpb[3][3]);
+void iauPnm80(double date1, double date2, double rmatpn[3][3]);
+void iauPom00(double xp, double yp, double sp, double rpom[3][3]);
+void iauPr00(double date1, double date2,
+             double *dpsipr, double *depspr);
+void iauPrec76(double date01, double date02,
+               double date11, double date12,
+               double *zeta, double *z, double *theta);
+double iauS00(double date1, double date2, double x, double y);
+double iauS00a(double date1, double date2);
+double iauS00b(double date1, double date2);
+double iauS06(double date1, double date2, double x, double y);
+double iauS06a(double date1, double date2);
+double iauSp00(double date1, double date2);
+void iauXy06(double date1, double date2, double *x, double *y);
+void iauXys00a(double date1, double date2,
+               double *x, double *y, double *s);
+void iauXys00b(double date1, double date2,
+               double *x, double *y, double *s);
+void iauXys06a(double date1, double date2,
+               double *x, double *y, double *s);
+
+/* Astronomy/RotationAndTime */
+double iauEe00(double date1, double date2, double epsa, double dpsi);
+double iauEe00a(double date1, double date2);
+double iauEe00b(double date1, double date2);
+double iauEe06a(double date1, double date2);
+double iauEect00(double date1, double date2);
+double iauEqeq94(double date1, double date2);
+double iauEra00(double dj1, double dj2);
+double iauGmst00(double uta, double utb, double tta, double ttb);
+double iauGmst06(double uta, double utb, double tta, double ttb);
+double iauGmst82(double dj1, double dj2);
+double iauGst00a(double uta, double utb, double tta, double ttb);
+double iauGst00b(double uta, double utb);
+double iauGst06(double uta, double utb, double tta, double ttb,
+                double rnpb[3][3]);
+double iauGst06a(double uta, double utb, double tta, double ttb);
+double iauGst94(double uta, double utb);
+
+/* Astronomy/SpaceMotion */
+int iauPvstar(double pv[2][3], double *ra, double *dec,
+              double *pmr, double *pmd, double *px, double *rv);
+int iauStarpv(double ra, double dec,
+              double pmr, double pmd, double px, double rv,
+              double pv[2][3]);
+
+/* Astronomy/StarCatalogs */
+void iauFk425(double r1950, double d1950,
+              double dr1950, double dd1950,
+              double p1950, double v1950,
+              double *r2000, double *d2000,
+              double *dr2000, double *dd2000,
+              double *p2000, double *v2000);
+void iauFk45z(double r1950, double d1950, double bepoch,
+              double *r2000, double *d2000);
+void iauFk524(double r2000, double d2000,
+              double dr2000, double dd2000,
+              double p2000, double v2000,
+              double *r1950, double *d1950,
+              double *dr1950, double *dd1950,
+              double *p1950, double *v1950);
+void iauFk52h(double r5, double d5,
+              double dr5, double dd5, double px5, double rv5,
+              double *rh, double *dh,
+              double *drh, double *ddh, double *pxh, double *rvh);
+void iauFk54z(double r2000, double d2000, double bepoch,
+              double *r1950, double *d1950,
+              double *dr1950, double *dd1950);
+void iauFk5hip(double r5h[3][3], double s5h[3]);
+void iauFk5hz(double r5, double d5, double date1, double date2,
+              double *rh, double *dh);
+void iauH2fk5(double rh, double dh,
+              double drh, double ddh, double pxh, double rvh,
+              double *r5, double *d5,
+              double *dr5, double *dd5, double *px5, double *rv5);
+void iauHfk5z(double rh, double dh, double date1, double date2,
+              double *r5, double *d5, double *dr5, double *dd5);
+int iauStarpm(double ra1, double dec1,
+              double pmr1, double pmd1, double px1, double rv1,
+              double ep1a, double ep1b, double ep2a, double ep2b,
+              double *ra2, double *dec2,
+              double *pmr2, double *pmd2, double *px2, double *rv2);
+
+/* Astronomy/EclipticCoordinates */
+void iauEceq06(double date1, double date2, double dl, double db,
+               double *dr, double *dd);
+void iauEcm06(double date1, double date2, double rm[3][3]);
+void iauEqec06(double date1, double date2, double dr, double dd,
+               double *dl, double *db);
+void iauLteceq(double epj, double dl, double db, double *dr, double *dd);
+void iauLtecm(double epj, double rm[3][3]);
+void iauLteqec(double epj, double dr, double dd, double *dl, double *db);
+
+/* Astronomy/GalacticCoordinates */
+void iauG2icrs(double dl, double db, double *dr, double *dd);
+void iauIcrs2g(double dr, double dd, double *dl, double *db);
+
+/* Astronomy/GeodeticGeocentric */
+int iauEform(int n, double *a, double *f);
+int iauGc2gd(int n, double xyz[3],
+             double *elong, double *phi, double *height);
+int iauGc2gde(double a, double f, double xyz[3],
+              double *elong, double *phi, double *height);
+int iauGd2gc(int n, double elong, double phi, double height,
+             double xyz[3]);
+int iauGd2gce(double a, double f,
+              double elong, double phi, double height, double xyz[3]);
+
+/* Astronomy/Timescales */
+int iauD2dtf(const char *scale, int ndp, double d1, double d2,
+             int *iy, int *im, int *id, int ihmsf[4]);
+int iauDat(int iy, int im, int id, double fd, double *deltat);
+double iauDtdb(double date1, double date2,
+               double ut, double elong, double u, double v);
+int iauDtf2d(const char *scale, int iy, int im, int id,
+             int ihr, int imn, double sec, double *d1, double *d2);
+int iauTaitt(double tai1, double tai2, double *tt1, double *tt2);
+int iauTaiut1(double tai1, double tai2, double dta,
+              double *ut11, double *ut12);
+int iauTaiutc(double tai1, double tai2, double *utc1, double *utc2);
+int iauTcbtdb(double tcb1, double tcb2, double *tdb1, double *tdb2);
+int iauTcgtt(double tcg1, double tcg2, double *tt1, double *tt2);
+int iauTdbtcb(double tdb1, double tdb2, double *tcb1, double *tcb2);
+int iauTdbtt(double tdb1, double tdb2, double dtr,
+             double *tt1, double *tt2);
+int iauTttai(double tt1, double tt2, double *tai1, double *tai2);
+int iauTttcg(double tt1, double tt2, double *tcg1, double *tcg2);
+int iauTttdb(double tt1, double tt2, double dtr,
+             double *tdb1, double *tdb2);
+int iauTtut1(double tt1, double tt2, double dt,
+             double *ut11, double *ut12);
+int iauUt1tai(double ut11, double ut12, double dta,
+              double *tai1, double *tai2);
+int iauUt1tt(double ut11, double ut12, double dt,
+             double *tt1, double *tt2);
+int iauUt1utc(double ut11, double ut12, double dut1,
+              double *utc1, double *utc2);
+int iauUtctai(double utc1, double utc2, double *tai1, double *tai2);
+int iauUtcut1(double utc1, double utc2, double dut1,
+              double *ut11, double *ut12);
+
+/* Astronomy/HorizonEquatorial */
+void iauAe2hd(double az, double el, double phi,
+              double *ha, double *dec);
+void iauHd2ae(double ha, double dec, double phi,
+              double *az, double *el);
+double iauHd2pa(double ha, double dec, double phi);
+
+/* Astronomy/Gnomonic */
+int iauTpors(double xi, double eta, double a, double b,
+             double *a01, double *b01, double *a02, double *b02);
+int iauTporv(double xi, double eta, double v[3],
+             double v01[3], double v02[3]);
+void iauTpsts(double xi, double eta, double a0, double b0,
+              double *a, double *b);
+void iauTpstv(double xi, double eta, double v0[3], double v[3]);
+int iauTpxes(double a, double b, double a0, double b0,
+             double *xi, double *eta);
+int iauTpxev(double v[3], double v0[3], double *xi, double *eta);
+
+/* VectorMatrix/AngleOps */
+void iauA2af(int ndp, double angle, char *sign, int idmsf[4]);
+void iauA2tf(int ndp, double angle, char *sign, int ihmsf[4]);
+int iauAf2a(char s, int ideg, int iamin, double asec, double *rad);
+double iauAnp(double a);
+double iauAnpm(double a);
+void iauD2tf(int ndp, double days, char *sign, int ihmsf[4]);
+int iauTf2a(char s, int ihour, int imin, double sec, double *rad);
+int iauTf2d(char s, int ihour, int imin, double sec, double *days);
+
+/* VectorMatrix/BuildRotations */
+void iauRx(double phi, double r[3][3]);
+void iauRy(double theta, double r[3][3]);
+void iauRz(double psi, double r[3][3]);
+
+/* VectorMatrix/CopyExtendExtract */
+void iauCp(double p[3], double c[3]);
+void iauCpv(double pv[2][3], double c[2][3]);
+void iauCr(double r[3][3], double c[3][3]);
+void iauP2pv(double p[3], double pv[2][3]);
+void iauPv2p(double pv[2][3], double p[3]);
+
+/* VectorMatrix/Initialization */
+void iauIr(double r[3][3]);
+void iauZp(double p[3]);
+void iauZpv(double pv[2][3]);
+void iauZr(double r[3][3]);
+
+/* VectorMatrix/MatrixOps */
+void iauRxr(double a[3][3], double b[3][3], double atb[3][3]);
+void iauTr(double r[3][3], double rt[3][3]);
+
+/* VectorMatrix/MatrixVectorProducts */
+void iauRxp(double r[3][3], double p[3], double rp[3]);
+void iauRxpv(double r[3][3], double pv[2][3], double rpv[2][3]);
+void iauTrxp(double r[3][3], double p[3], double trp[3]);
+void iauTrxpv(double r[3][3], double pv[2][3], double trpv[2][3]);
+
+/* VectorMatrix/RotationVectors */
+void iauRm2v(double r[3][3], double w[3]);
+void iauRv2m(double w[3], double r[3][3]);
+
+/* VectorMatrix/SeparationAndAngle */
+double iauPap(double a[3], double b[3]);
+double iauPas(double al, double ap, double bl, double bp);
+double iauSepp(double a[3], double b[3]);
+double iauSeps(double al, double ap, double bl, double bp);
+
+/* VectorMatrix/SphericalCartesian */
+void iauC2s(double p[3], double *theta, double *phi);
+void iauP2s(double p[3], double *theta, double *phi, double *r);
+void iauPv2s(double pv[2][3],
+             double *theta, double *phi, double *r,
+             double *td, double *pd, double *rd);
+void iauS2c(double theta, double phi, double c[3]);
+void iauS2p(double theta, double phi, double r, double p[3]);
+void iauS2pv(double theta, double phi, double r,
+             double td, double pd, double rd,
+             double pv[2][3]);
+
+/* VectorMatrix/VectorOps */
+double iauPdp(double a[3], double b[3]);
+double iauPm(double p[3]);
+void iauPmp(double a[3], double b[3], double amb[3]);
+void iauPn(double p[3], double *r, double u[3]);
+void iauPpp(double a[3], double b[3], double apb[3]);
+void iauPpsp(double a[3], double s, double b[3], double apsb[3]);
+void iauPvdpv(double a[2][3], double b[2][3], double adb[2]);
+void iauPvm(double pv[2][3], double *r, double *s);
+void iauPvmpv(double a[2][3], double b[2][3], double amb[2][3]);
+void iauPvppv(double a[2][3], double b[2][3], double apb[2][3]);
+void iauPvu(double dt, double pv[2][3], double upv[2][3]);
+void iauPvup(double dt, double pv[2][3], double p[3]);
+void iauPvxpv(double a[2][3], double b[2][3], double axb[2][3]);
+void iauPxp(double a[3], double b[3], double axb[3]);
+void iauS2xpv(double s1, double s2, double pv[2][3], double spv[2][3]);
+void iauSxp(double s, double p[3], double sp[3]);
+void iauSxpv(double s, double pv[2][3], double spv[2][3]);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
+
+/*----------------------------------------------------------------------
+**
+**  Copyright (C) 2023
+**  Standards of Fundamental Astronomy Board
+**  of the International Astronomical Union.
+**
+**  =====================
+**  SOFA Software License
+**  =====================
+**
+**  NOTICE TO USER:
+**
+**  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING SIX TERMS AND
+**  CONDITIONS WHICH APPLY TO ITS USE.
+**
+**  1. The Software is owned by the IAU SOFA Board ("SOFA").
+**
+**  2. Permission is granted to anyone to use the SOFA software for any
+**     purpose, including commercial applications, free of charge and
+**     without payment of royalties, subject to the conditions and
+**     restrictions listed below.
+**
+**  3. You (the user) may copy and distribute SOFA source code to others,
+**     and use and adapt its code and algorithms in your own software,
+**     on a world-wide, royalty-free basis.  That portion of your
+**     distribution that does not consist of intact and unchanged copies
+**     of SOFA source code files is a "derived work" that must comply
+**     with the following requirements:
+**
+**     a) Your work shall be marked or carry a statement that it
+**        (i) uses routines and computations derived by you from
+**        software provided by SOFA under license to you; and
+**        (ii) does not itself constitute software provided by and/or
+**        endorsed by SOFA.
+**
+**     b) The source code of your derived work must contain descriptions
+**        of how the derived work is based upon, contains and/or differs
+**        from the original SOFA software.
+**
+**     c) The names of all routines in your derived work shall not
+**        include the prefix "iau" or "sofa" or trivial modifications
+**        thereof such as changes of case.
+**
+**     d) The origin of the SOFA components of your derived work must
+**        not be misrepresented;  you must not claim that you wrote the
+**        original software, nor file a patent application for SOFA
+**        software or algorithms embedded in the SOFA software.
+**
+**     e) These requirements must be reproduced intact in any source
+**        distribution and shall apply to anyone to whom you have
+**        granted a further right to modify the source code of your
+**        derived work.
+**
+**     Note that, as originally distributed, the SOFA software is
+**     intended to be a definitive implementation of the IAU standards,
+**     and consequently third-party modifications are discouraged.  All
+**     variations, no matter how minor, must be explicitly marked as
+**     such, as explained above.
+**
+**  4. You shall not cause the SOFA software to be brought into
+**     disrepute, either by misuse, or use for inappropriate tasks, or
+**     by inappropriate modification.
+**
+**  5. The SOFA software is provided "as is" and SOFA makes no warranty
+**     as to its use or performance.   SOFA does not and cannot warrant
+**     the performance or results which the user may obtain by using the
+**     SOFA software.  SOFA makes no warranties, express or implied, as
+**     to non-infringement of third party rights, merchantability, or
+**     fitness for any particular purpose.  In no event will SOFA be
+**     liable to the user for any consequential, incidental, or special
+**     damages, including any lost profits or lost savings, even if a
+**     SOFA representative has been advised of such damages, or for any
+**     claim by any third party.
+**
+**  6. The provision of any version of the SOFA software under the terms
+**     and conditions specified herein does not imply that future
+**     versions will also be made available under the same terms and
+**     conditions.
+*
+**  In any published work or commercial product which uses the SOFA
+**  software directly, acknowledgement (see www.iausofa.org) is
+**  appreciated.
+**
+**  Correspondence concerning SOFA software should be addressed as
+**  follows:
+**
+**      By email:  sofa@ukho.gov.uk
+**      By post:   IAU SOFA Center
+**                 HM Nautical Almanac Office
+**                 UK Hydrographic Office
+**                 Admiralty Way, Taunton
+**                 Somerset, TA1 2DN
+**                 United Kingdom
+**
+**--------------------------------------------------------------------*/
diff --git a/src/sofam.h b/src/sofam.h
new file mode 100644
index 000000000..9279f7277
--- /dev/null
+++ b/src/sofam.h
@@ -0,0 +1,222 @@
+#ifndef SOFAMHDEF
+#define SOFAMHDEF
+
+/*
+**  - - - - - - - -
+**   s o f a m . h
+**  - - - - - - - -
+**
+**  Macros used by SOFA library.
+**
+**  This file is part of the International Astronomical Union's
+**  SOFA (Standards of Fundamental Astronomy) software collection.
+**
+**  Please note that the constants defined below are to be used only in
+**  the context of the SOFA software, and have no other official IAU
+**  status.  In addition, self consistency is not guaranteed.
+**
+**  This revision:   2021 February 24
+**
+**  SOFA release 2023-10-11
+**
+**  Copyright (C) 2023 IAU SOFA Board.  See notes at end.
+*/
+
+/* Pi */
+#define DPI (3.141592653589793238462643)
+
+/* 2Pi */
+#define D2PI (6.283185307179586476925287)
+
+/* Radians to degrees */
+#define DR2D (57.29577951308232087679815)
+
+/* Degrees to radians */
+#define DD2R (1.745329251994329576923691e-2)
+
+/* Radians to arcseconds */
+#define DR2AS (206264.8062470963551564734)
+
+/* Arcseconds to radians */
+#define DAS2R (4.848136811095359935899141e-6)
+
+/* Seconds of time to radians */
+#define DS2R (7.272205216643039903848712e-5)
+
+/* Arcseconds in a full circle */
+#define TURNAS (1296000.0)
+
+/* Milliarcseconds to radians */
+#define DMAS2R (DAS2R / 1e3)
+
+/* Length of tropical year B1900 (days) */
+#define DTY (365.242198781)
+
+/* Seconds per day. */
+#define DAYSEC (86400.0)
+
+/* Days per Julian year */
+#define DJY (365.25)
+
+/* Days per Julian century */
+#define DJC (36525.0)
+
+/* Days per Julian millennium */
+#define DJM (365250.0)
+
+/* Reference epoch (J2000.0), Julian Date */
+#define DJ00 (2451545.0)
+
+/* Julian Date of Modified Julian Date zero */
+#define DJM0 (2400000.5)
+
+/* Reference epoch (J2000.0), Modified Julian Date */
+#define DJM00 (51544.5)
+
+/* 1977 Jan 1.0 as MJD */
+#define DJM77 (43144.0)
+
+/* TT minus TAI (s) */
+#define TTMTAI (32.184)
+
+/* Astronomical unit (m, IAU 2012) */
+#define DAU (149597870.7e3)
+
+/* Speed of light (m/s) */
+#define CMPS 299792458.0
+
+/* Light time for 1 au (s) */
+#define AULT (DAU/CMPS)
+
+/* Speed of light (au per day) */
+#define DC (DAYSEC/AULT)
+
+/* L_G = 1 - d(TT)/d(TCG) */
+#define ELG (6.969290134e-10)
+
+/* L_B = 1 - d(TDB)/d(TCB), and TDB (s) at TAI 1977/1/1.0 */
+#define ELB (1.550519768e-8)
+#define TDB0 (-6.55e-5)
+
+/* Schwarzschild radius of the Sun (au) */
+/* = 2 * 1.32712440041e20 / (2.99792458e8)^2 / 1.49597870700e11 */
+#define SRS 1.97412574336e-8
+
+/* dint(A) - truncate to nearest whole number towards zero (double) */
+#define dint(A) ((A)<0.0?ceil(A):floor(A))
+
+/* dnint(A) - round to nearest whole number (double) */
+#define dnint(A) (fabs(A)<0.5?0.0\
+                                :((A)<0.0?ceil((A)-0.5):floor((A)+0.5)))
+
+/* dsign(A,B) - magnitude of A with sign of B (double) */
+#define dsign(A,B) ((B)<0.0?-fabs(A):fabs(A))
+
+/* max(A,B) - larger (most +ve) of two numbers (generic) */
+#define gmax(A,B) (((A)>(B))?(A):(B))
+
+/* min(A,B) - smaller (least +ve) of two numbers (generic) */
+#define gmin(A,B) (((A)<(B))?(A):(B))
+
+/* Reference ellipsoids */
+#define WGS84 1
+#define GRS80 2
+#define WGS72 3
+
+#endif
+
+/*----------------------------------------------------------------------
+**
+**  Copyright (C) 2023
+**  Standards of Fundamental Astronomy Board
+**  of the International Astronomical Union.
+**
+**  =====================
+**  SOFA Software License
+**  =====================
+**
+**  NOTICE TO USER:
+**
+**  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING SIX TERMS AND
+**  CONDITIONS WHICH APPLY TO ITS USE.
+**
+**  1. The Software is owned by the IAU SOFA Board ("SOFA").
+**
+**  2. Permission is granted to anyone to use the SOFA software for any
+**     purpose, including commercial applications, free of charge and
+**     without payment of royalties, subject to the conditions and
+**     restrictions listed below.
+**
+**  3. You (the user) may copy and distribute SOFA source code to others,
+**     and use and adapt its code and algorithms in your own software,
+**     on a world-wide, royalty-free basis.  That portion of your
+**     distribution that does not consist of intact and unchanged copies
+**     of SOFA source code files is a "derived work" that must comply
+**     with the following requirements:
+**
+**     a) Your work shall be marked or carry a statement that it
+**        (i) uses routines and computations derived by you from
+**        software provided by SOFA under license to you; and
+**        (ii) does not itself constitute software provided by and/or
+**        endorsed by SOFA.
+**
+**     b) The source code of your derived work must contain descriptions
+**        of how the derived work is based upon, contains and/or differs
+**        from the original SOFA software.
+**
+**     c) The names of all routines in your derived work shall not
+**        include the prefix "iau" or "sofa" or trivial modifications
+**        thereof such as changes of case.
+**
+**     d) The origin of the SOFA components of your derived work must
+**        not be misrepresented;  you must not claim that you wrote the
+**        original software, nor file a patent application for SOFA
+**        software or algorithms embedded in the SOFA software.
+**
+**     e) These requirements must be reproduced intact in any source
+**        distribution and shall apply to anyone to whom you have
+**        granted a further right to modify the source code of your
+**        derived work.
+**
+**     Note that, as originally distributed, the SOFA software is
+**     intended to be a definitive implementation of the IAU standards,
+**     and consequently third-party modifications are discouraged.  All
+**     variations, no matter how minor, must be explicitly marked as
+**     such, as explained above.
+**
+**  4. You shall not cause the SOFA software to be brought into
+**     disrepute, either by misuse, or use for inappropriate tasks, or
+**     by inappropriate modification.
+**
+**  5. The SOFA software is provided "as is" and SOFA makes no warranty
+**     as to its use or performance.   SOFA does not and cannot warrant
+**     the performance or results which the user may obtain by using the
+**     SOFA software.  SOFA makes no warranties, express or implied, as
+**     to non-infringement of third party rights, merchantability, or
+**     fitness for any particular purpose.  In no event will SOFA be
+**     liable to the user for any consequential, incidental, or special
+**     damages, including any lost profits or lost savings, even if a
+**     SOFA representative has been advised of such damages, or for any
+**     claim by any third party.
+**
+**  6. The provision of any version of the SOFA software under the terms
+**     and conditions specified herein does not imply that future
+**     versions will also be made available under the same terms and
+**     conditions.
+*
+**  In any published work or commercial product which uses the SOFA
+**  software directly, acknowledgement (see www.iausofa.org) is
+**  appreciated.
+**
+**  Correspondence concerning SOFA software should be addressed as
+**  follows:
+**
+**      By email:  sofa@ukho.gov.uk
+**      By post:   IAU SOFA Center
+**                 HM Nautical Almanac Office
+**                 UK Hydrographic Office
+**                 Admiralty Way, Taunton
+**                 Somerset, TA1 2DN
+**                 United Kingdom
+**
+**--------------------------------------------------------------------*/
diff --git a/src/tests.cc b/src/tests.cc
index 7c31c3f19..74116d5fe 100644
--- a/src/tests.cc
+++ b/src/tests.cc
@@ -10118,8 +10118,10 @@ void tests::date_operations()
         .test(CLEAR, "19690217.03550197_date", ENTER)
         .expect("Mon 17/Feb/1969, 3:55:01.97");
     step("Displaying invalid date and time")
+        .test(ID_HideDayOfWeek, ENTER)
         .test(CLEAR, "999999999.99999999_date", ENTER)
-        .expect("Sat 99/99/99999, 99:99:99.99");
+        .expect("99/99/99999, 99:99:99.99")
+        .test(ID_ShowDayOfWeek, ENTER);
 
     step("Difference between two dates using DDays")
         .test(CLEAR, "20230908", ENTER)
@@ -10171,10 +10173,13 @@ void tests::date_operations()
 
     step("Converting from a date to a Julian day number")
         .test(CLEAR, "20250919_date", ENTER, ID_JulianDayNumber)
-        .expect("2 460 938");
+        .expect("2 460 937 ¹/₂");
     step("Converting from a Julian day number to a date")
         .test(CLEAR, "1 000 000", ENTER, ID_DateFromJulianDayNumber)
-        .expect("Tue 21/Oct/1975 BC");
+        .expect("Tue 21/Oct/-1975, 12:00:00");
+    step("Change to BCECE representation")
+        .test(ID_BCECEyearNumbering)
+        .expect("Tue 21/Oct/1976 BCE, 12:00:00");
     step("Adding invalid dates")
         .test(CLEAR, "1 2 DATE+", ENTER)
         .error("Invalid date");