Release 3

dcf21 · Aug 23, 2024 · 82ea99a · 82ea99a
1 parent c1c2b10
commit 82ea99a
Showing 41 changed files with 291 additions and 125 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -1,6 +1,8 @@
 cmake_minimum_required(VERSION 3.3)
 project(ephemeris_compute)
 
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O0 -Wall -std=c99 -D SRCDIR='\"${CMAKE_SOURCE_DIR}/src/\"' -D DCFVERSION='\"x\"' -D DEBUG=0 -D MEMDEBUG1=0 -D MEMDEBUG2=0")
+
 include_directories(src)
 include_directories(src/argparse)
 include_directories(src/coreUtils)

diff --git a/Dockerfile b/Dockerfile
@@ -1,16 +1,20 @@
-# Use Python 3.10 running on Debian Bullseye
-FROM python:3.10-bullseye
+# Use Python 3.12 running on Ubuntu 24.04
+FROM ubuntu:24.04
 
 # Install required libraries
 RUN apt-get update
-RUN apt-get install -y wget gzip libgsl-dev ; apt-get clean
+RUN apt-get install -y apt-utils dialog file git vim python3 python3-dev \
+                       build-essential make gcc wget gzip libgsl-dev \
+                       pkg-config \
+                       ; apt-get clean
 
 # Copy code into container
 WORKDIR /
 ADD . ephemeris-compute
 
 # Fetch data
 WORKDIR /ephemeris-compute
+RUN /ephemeris-compute/prettymake clean
 RUN /ephemeris-compute/setup.sh
 
 # Check that binary quick-lookup files are generated

diff --git a/Makefile b/Makefile
@@ -1,7 +1,7 @@
 # Makefile for EphemerisCompute
 #
 # -------------------------------------------------
-# Copyright 2015-2020 Dominic Ford
+# Copyright 2015-2024 Dominic Ford
 #
 # This file is part of EphemerisCompute.
 #
@@ -21,8 +21,8 @@
 
 CWD=$(shell pwd)
 
-VERSION = 2.0
-DATE    = 09/06/2019
+VERSION = 3.0
+DATE    = 16/08/2024
 PATHLINK= /
 
 WARNINGS    = -Wall -Wno-format-truncation -Wno-unknown-pragmas

diff --git a/README.md b/README.md
@@ -29,7 +29,7 @@ Unix-like operating systems.
 ### License
 
 This code is distributed under the Gnu General Public License. It is (C)
-Dominic Ford 2010 - 2022.
+Dominic Ford 2010 - 2024.
 
 ### Set up
 
@@ -47,14 +47,14 @@ provided to build a selection of example starcharts:
 
 ```
 docker compose build
-docker compose run ephmeris-compute-de430
+docker compose run ephemeris-compute-de430
 ```
 
 This produces a single demo ephemeris. To make other ephemerides, open a shell
 within the Docker container as follows:
 
 ```
-docker run -it ephmeris-compute-de430:v1 /bin/bash
+docker run -it ephemeris-compute-de430:v3 /bin/bash
 ```
 
 ### Producing an ephemeris
@@ -124,6 +124,12 @@ This section lists the names which are recognised by the `--objects` command-lin
 * `CJ95O010`. Comets may be referred to by their Minor Planet Center designations
 * `C<n>`: Comer number `n`. `n` is the line number within the file [Soft00Cmt.txt](http://www.minorplanetcenter.net/iau/Ephemerides/Comets/Soft00Cmt.txt), downloaded from the Minor Planet Center.
 
+### Change history
+
+**Version 3.0** (23 Aug 2024) - Added corrections for light travel time and annual aberration. Improvements to numerical stability when calculating hyperbolic orbits.
+
+**Version 2.0** (16 Oct 2022) - Initial public release.
+
 ## Author
 
 This code was developed by Dominic Ford <https://dcford.org.uk>. It is distributed under the Gnu General Public License V3.

diff --git a/constellations/constellation_names.dat b/constellations/constellation_names.dat
@@ -1,7 +1,7 @@
 # constellation_names.dat
 #
 # -------------------------------------------------
-# Copyright 2015-2019 Dominic Ford
+# Copyright 2015-2024 Dominic Ford
 #
 # This file is part of EphemerisCompute.
 #

diff --git a/docker-compose.yml b/docker-compose.yml
@@ -1,6 +1,5 @@
-version: '3'
 services:
-  ephmeris-compute-de430:
+  ephemeris-compute-de430:
     build: .
-    image: ephmeris-compute-de430:v1
+    image: ephemeris-compute-de430:v3
     command: ./bin/ephem.bin
diff --git a/src/asteroids.c b/src/asteroids.c
@@ -2,7 +2,7 @@
 // Dominic Ford
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/coreUtils/asciiDouble.c b/src/coreUtils/asciiDouble.c
@@ -1,7 +1,7 @@
 // asciiDouble.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/coreUtils/asciiDouble.h b/src/coreUtils/asciiDouble.h
@@ -1,7 +1,7 @@
 // asciiDouble.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/coreUtils/errorReport.c b/src/coreUtils/errorReport.c
@@ -1,7 +1,7 @@
 // errorReport.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/coreUtils/errorReport.h b/src/coreUtils/errorReport.h
@@ -1,7 +1,7 @@
 // errorReport.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/coreUtils/makeRasters.c b/src/coreUtils/makeRasters.c
@@ -1,7 +1,7 @@
 // makeRasters.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/coreUtils/makeRasters.h b/src/coreUtils/makeRasters.h
@@ -1,7 +1,7 @@
 // makeRasters.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/coreUtils/strConstants.h b/src/coreUtils/strConstants.h
@@ -1,7 +1,7 @@
 // strConstants.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/ephemCalc/constellations.c b/src/ephemCalc/constellations.c
@@ -1,7 +1,7 @@
 // constellations.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //
@@ -58,7 +58,8 @@ typedef struct {
 } constel_desc;
 
 //! constel_desc constel_data - Storage for the 88 constellations which make up the night sky
-static constel_desc constel_data[90];
+#define MAX_CONSTELLATIONS 90
+static constel_desc constel_data[MAX_CONSTELLATIONS];
 
 //! Nconstel - A counter for the number of constellations we have loaded so far
 static int Nconstel = 0;
@@ -94,7 +95,8 @@ static double dWind(double RA, double Dec, double RA0, double Dec0, double RA1,
     double xC1 = xB1;
     double yC0 = yB0 * cos(-a) + zB0 * sin(-a);
     double yC1 = yB1 * cos(-a) + zB1 * sin(-a);
-//double zC0 = yB0*sin(  a) + zB0*cos(- a);      double zC1 = yB1*sin(  a) + zB1*cos(- a);
+    // double zC0 = yB0*sin(  a) + zB0*cos(- a);
+    // double zC1 = yB1*sin(  a) + zB1*cos(- a);
 
     double dW;
 
@@ -108,15 +110,25 @@ static double dWind(double RA, double Dec, double RA0, double Dec0, double RA1,
 
 void constellations_init() {
     FILE *file;
-    int i;
-    char line[FNAME_LENGTH], *scan, constellation[6] = "@@@";
-
-    Nconstel = 0;
-
-    // Scan through catalogue of bound
-    file = fopen(SRCDIR "../constellations/bound_20.dat", "r");
-    if (file == NULL) ephem_fatal(__FILE__, __LINE__, "Could not open constellation boundary data");
+    char line[FNAME_LENGTH], *scan;
+    char constellation[6] = "@@@"; // Name of constellation we're currently reading in
+
+    Nconstel = 0; // Constellation counter. NB: This is 1 when we're reading in constel_data[0]
+
+    // Scan through catalogue of constellation boundaries
+    {
+        char in_path[FNAME_LENGTH];
+        snprintf(in_path, FNAME_LENGTH, "%s", SRCDIR "../constellations/bound_20.dat");
+        file = fopen(in_path, "r");
+        if (file == NULL) {
+            char buffer[LSTR_LENGTH];
+            snprintf(buffer, LSTR_LENGTH, "Could not open constellation boundary data from <%s>.", in_path);
+            ephem_fatal(__FILE__, __LINE__, buffer);
+            exit(1);
+        }
+    }
 
+    // Loop over lines of constellation boundary data
     while ((!feof(file)) && (!ferror(file))) {
         double ra, dec;
         file_readline(file, line);
@@ -129,15 +141,25 @@ void constellations_init() {
         dec = get_float(scan, NULL);
         if (line[11] == '-') dec = -dec;
 
+        // If constellation name has changed since previous entry, we have a new constellation.
         if (strncmp(line + 23, constellation, 4) != 0) {
-            Nconstel++;
+            Nconstel++; // Advance constellation counter
+
+            // Check for data structure overflow
+            if (Nconstel >= MAX_CONSTELLATIONS) {
+                ephem_fatal(__FILE__, __LINE__, "Too many constellations; bound_20.dat is probably corrupted.");
+            }
+
+            // Initialise new constellation data structure
+            // Note -1 because Nconstel is 1 when reading first record, etc
             constel_data[Nconstel - 1].ShortNameSet = 1;
             constel_data[Nconstel - 1].LongNameSet = 0;
             constel_data[Nconstel - 1].Npoints = 0;
-            strncpy(constel_data[Nconstel - 1].ShortName, line + 23, 4);
-            constel_data[Nconstel - 1].ShortName[5] = '\0';
-            strncpy(constellation, line + 23, 4);
-            constellation[5] = '\0';
+
+            // shortName is 4-letter constellation abbrev
+            // Note, "SER1" has four letter abbrevs. Most are filed as "AND " with trailing space.
+            snprintf(constel_data[Nconstel - 1].ShortName, 5, "%s", line + 23);
+            snprintf(constellation, 5, "%s", line + 23);  // Most are filed as "AND " with trailing space
         }
 
         constel_data[Nconstel - 1].point[constel_data[Nconstel - 1].Npoints].RA = ra / 12. * M_PI;
@@ -148,32 +170,49 @@ void constellations_init() {
     fclose(file);
 
     // Scan through list of full names
-    file = fopen(SRCDIR "../constellations/constellation_names.dat", "r");
-    if (file == NULL) ephem_fatal(__FILE__, __LINE__, "Could not open constellation name data");
+    {
+        char in_path[FNAME_LENGTH];
+        snprintf(in_path, FNAME_LENGTH, "%s", SRCDIR "../constellations/constellation_names.dat");
+        file = fopen(in_path, "r");
+        if (file == NULL) {
+            char buffer[LSTR_LENGTH];
+            snprintf(buffer, LSTR_LENGTH, "Could not open constellation name data from <%s>.", in_path);
+            ephem_fatal(__FILE__, __LINE__, buffer);
+            exit(1);
+        }
+    }
 
     while ((!feof(file)) && (!ferror(file))) {
         char ShortName[6], LongName[64];
         int i, done;
 
-        // Read short/long name pair
+        // Lines take the form:
+        // short_name  long_name
         file_readline(file, line);
         if ((line[0] == '#') || (strlen(line) < 4)) continue; // Comment line
+
+        // Read short name
         scan = line + 0;
         while ((scan[0] > '\0') && (scan[0] <= ' ')) scan++;
         for (i = 0; scan[0] > ' '; i++, scan++) ShortName[i] = *scan;
+
+        // Pad short name to ensure it's 4 characters long to match abbreviation above
         while (i < 4) ShortName[i++] = ' ';
         ShortName[i] = '\0';
+
+        // Read long name
         while ((scan[0] > '\0') && (scan[0] <= ' ')) scan++;
         for (i = 0; scan[0] > ' '; i++, scan++) {
-            if ((*scan < 'a') && (i > 0)) LongName[i++] = '@';
+            if ((*scan < 'a') && (i > 0)) LongName[i++] = ' ';
             LongName[i] = *scan;
         }
         LongName[i] = '\0';
 
-        // Search for entry in constel_data
+        // Search for entry in constel_data with matching abbreviated name
         for (done = i = 0; i < Nconstel; i++) {
             if (!constel_data[i].ShortNameSet) ephem_fatal(__FILE__, __LINE__, "ShortName of constellation not set");
             if (strcmp(ShortName, constel_data[i].ShortName) == 0) {
+                // If we've found constellation with matching abbreviation, insert it's long name into record
                 strcpy(constel_data[i].LongName, LongName);
                 constel_data[i].LongNameSet = 1;
                 done = 1;
@@ -186,7 +225,7 @@ void constellations_init() {
         }
     }
 
-    for (i = 0; i < Nconstel; i++)
+    for (int i = 0; i < Nconstel; i++)
         if (!constel_data[i].LongNameSet) {
             sprintf(temp_err_string, "Could not find long name for constellation'%s'", constel_data[i].ShortName);
             ephem_fatal(__FILE__, __LINE__, temp_err_string);

diff --git a/src/ephemCalc/constellations.h b/src/ephemCalc/constellations.h
@@ -1,7 +1,7 @@
 // constellations.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/ephemCalc/jpl.c b/src/ephemCalc/jpl.c
@@ -1,7 +1,7 @@
 // jpl.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //
@@ -25,6 +25,7 @@
 #include <string.h>
 
 #include <gsl/gsl_math.h>
+#include <gsl/gsl_const_mksa.h>
 
 #include "coreUtils/asciiDouble.h"
 #include "coreUtils/errorReport.h"
@@ -33,8 +34,6 @@
 #include "listTools/ltDict.h"
 #include "listTools/ltMemory.h"
 
-#include "settings/settings.h"
-
 #include "jpl.h"
 #include "orbitalElements.h"
 #include "magnitudeEstimate.h"
@@ -128,7 +127,7 @@ int JPL_ReadBinaryData() {
     dcffread((void *) JPL_ShapeData, sizeof(int), 13 * 3, JPL_EphemFile);
 
     // We have now reached the actual ephemeris data. We don't load this into RAM since it is large and this would
-    // take time. Instead store a pointer to the offset of the start of the ephemeris from the beginning of file.
+    // take time. Instead, store a pointer to the offset of the start of the ephemeris from the beginning of file.
     JPL_EphemData_offset = (int) ftell(JPL_EphemFile);
 
     // Allocate memory to use to store ephemeris, as we load it
@@ -667,6 +666,10 @@ void jpl_computeEphemeris(settings *i, int bodyId, double jd, double *x, double
     // Boolean flags indicating whether this is the Earth, Sun or Moon (which need special treatment)
     int is_moon = 0, is_earth = 0, is_sun = 0;
 
+    double EMX_future, EMY_future, EMZ_future; // Position of the Earth-Moon centre of mass
+    double moon_pos_x_future, moon_pos_y_future, moon_pos_z_future;
+    double earth_pos_x_future, earth_pos_y_future, earth_pos_z_future;
+
     // Body 19 is the Earth.
     // DE430 gives us the Earth/Moon barycentre (body 2), from which we subtract a small fraction of the Moon's
     // offset (body 9) to get the Earth's centre of mass
@@ -708,9 +711,6 @@ void jpl_computeEphemeris(settings *i, int bodyId, double jd, double *x, double
     const double moon_mass = 0.1093189565989898e-10;
     const double moon_earth_mass_ratio = moon_mass / (moon_mass + earth_mass);
 
-    // Look up the Sun's position
-    jpl_computeXYZ(10, jd, &sun_pos_x, &sun_pos_y, &sun_pos_z);
-
     // Look up the Earth-Moon centre of mass position
     jpl_computeXYZ(2, jd, &EMX, &EMY, &EMZ);
 
@@ -722,39 +722,98 @@ void jpl_computeEphemeris(settings *i, int bodyId, double jd, double *x, double
     earth_pos_y = EMY - moon_earth_mass_ratio * moon_pos_y;
     earth_pos_z = EMZ - moon_earth_mass_ratio * moon_pos_z;
 
+    // Look up the Sun's position, taking light travel time into account
+    {
+        jpl_computeXYZ(10, jd, &sun_pos_x, &sun_pos_y, &sun_pos_z);
+
+        // Calculate light travel time
+        const double distance = gsl_hypot3(sun_pos_x - earth_pos_x,
+                                           sun_pos_y - earth_pos_y,
+                                           sun_pos_z - earth_pos_z);  // AU
+        const double light_travel_time = distance * GSL_CONST_MKSA_ASTRONOMICAL_UNIT / GSL_CONST_MKSA_SPEED_OF_LIGHT;
+
+        // Look up position of requested object at the time the light left the object
+        jpl_computeXYZ(10, jd - light_travel_time / 86400, &sun_pos_x, &sun_pos_y, &sun_pos_z);
+    }
+
+    // If the user's query was about the Earth, we already know its position
     if (is_earth) {
-        // If the user's query was about the Earth, we already know its position
         *x = earth_pos_x;
         *y = earth_pos_y;
         *z = earth_pos_z;
     }
+
+        // If the user's query was about the Sun, we already know that position too
     else if (is_sun) {
-        // If the user's query was about the Sun, we already know that position too!
         *x = sun_pos_x;
         *y = sun_pos_y;
         *z = sun_pos_z;
     }
+
+        // If the user's query was about the Moon, we already know that position too
     else if (is_moon) {
-        // If the user's query was about the Moon, we already know that position too!
-        *x = moon_pos_x;
-        *y = moon_pos_y;
-        *z = moon_pos_z;
+        *x = moon_pos_x + earth_pos_x;
+        *y = moon_pos_y + earth_pos_y;
+        *z = moon_pos_z + earth_pos_z;
     }
+
+        // Otherwise we need to query DE430 for the particular object the user was looking for,
+        // taking light travel time into account
     else {
-        // ... otherwise we need to query DE430 for the particular object the user was looking for
+        // Calculate position of requested object at specified time
         jpl_computeXYZ(bodyId, jd, x, y, z);
+
+        // Calculate light travel time
+        const double distance = gsl_hypot3(*x - earth_pos_x, *y - earth_pos_y, *z - earth_pos_z);  // AU
+        const double light_travel_time = distance * GSL_CONST_MKSA_ASTRONOMICAL_UNIT / GSL_CONST_MKSA_SPEED_OF_LIGHT;
+
+        // Look up position of requested object at the time the light left the object
+        jpl_computeXYZ(bodyId, jd - light_travel_time / 86400, x, y, z);
     }
 
-    // If the body was the Moon, don't forget the add the position of the Earth-Moon centre of mass
-    if (is_moon) {
-        *x += earth_pos_x;
-        *y += earth_pos_y;
-        *z += earth_pos_z;
+    // Look up the Earth-Moon centre of mass position, a short time in the future
+    // We use this to calculate the Earth's velocity vector, which is needed to correct for aberration
+    // (see eqn 7.119 of the Explanatory Supplement)
+    const double eb_dot_timestep = 1e-6; // days
+    const double eb_dot_timestep_sec = eb_dot_timestep * 86400;
+    jpl_computeXYZ(2, jd + eb_dot_timestep, &EMX_future, &EMY_future, &EMZ_future);
+    jpl_computeXYZ(9, jd + eb_dot_timestep, &moon_pos_x_future, &moon_pos_y_future, &moon_pos_z_future);
+    earth_pos_x_future = EMX_future - moon_earth_mass_ratio * moon_pos_x_future;
+    earth_pos_y_future = EMY_future - moon_earth_mass_ratio * moon_pos_y_future;
+    earth_pos_z_future = EMZ_future - moon_earth_mass_ratio * moon_pos_z_future;
+
+    // Equation (7.118) of the Explanatory Supplement - correct for aberration
+    if (!is_earth) {
+        const double u1[3] = {
+                *x - earth_pos_x,
+                *y - earth_pos_y,
+                *z - earth_pos_z
+        };
+        const double u1_mag = gsl_hypot3(u1[0], u1[1], u1[2]);
+        const double u[3] = {u1[0] / u1_mag, u1[1] / u1_mag, u1[2] / u1_mag};
+        const double eb_dot[3] = {
+                earth_pos_x_future - earth_pos_x,
+                earth_pos_y_future - earth_pos_y,
+                earth_pos_z_future - earth_pos_z
+        };
+
+        // Speed of light in AU per time step
+        const double c = GSL_CONST_MKSA_SPEED_OF_LIGHT / GSL_CONST_MKSA_ASTRONOMICAL_UNIT * eb_dot_timestep_sec;
+        const double V[3] = {eb_dot[0] / c, eb_dot[1] / c, eb_dot[2] / c};
+        const double V_mag = gsl_hypot3(V[0], V[1], V[2]);
+        const double beta = sqrt(1 - gsl_pow_2(V_mag));
+        const double f1 = u[0] * V[0] + u[1] * V[1] + u[2] * V[2];
+        const double f2 = 1 + f1 / (1 + beta);
+
+        // Correct for aberration
+        *x = earth_pos_x + (beta * u1[0] + f2 * u1_mag * V[0]) / (1 + f1);
+        *y = earth_pos_y + (beta * u1[1] + f2 * u1_mag * V[1]) / (1 + f1);
+        *z = earth_pos_z + (beta * u1[2] + f2 * u1_mag * V[2]) / (1 + f1);
     }
 
     // Populate other quantities, like the brightness, RA and Dec of the object, based on its XYZ position
-    magnitudeEstimate(bodyId, *x, *y, *z, earth_pos_x, earth_pos_y, earth_pos_z,
-                      sun_pos_x, sun_pos_y, sun_pos_z, ra, dec, mag, phase, angSize, phySize,
+    magnitudeEstimate(bodyId, *x, *y, *z, earth_pos_x, earth_pos_y, earth_pos_z, sun_pos_x, sun_pos_y, sun_pos_z, ra,
+                      dec, mag, phase, angSize, phySize,
                       albedo, sunDist, earthDist, sunAngDist, theta_ESO, eclipticLongitude, eclipticLatitude,
                       eclipticDistance, i);
 }
diff --git a/src/ephemCalc/jpl.h b/src/ephemCalc/jpl.h
@@ -1,7 +1,7 @@
 // jpl.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/ephemCalc/magnitudeEstimate.c b/src/ephemCalc/magnitudeEstimate.c
@@ -1,7 +1,7 @@
 // magnitudeEstimate.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/ephemCalc/magnitudeEstimate.h b/src/ephemCalc/magnitudeEstimate.h
@@ -1,7 +1,7 @@
 // magnitudeEstimate.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/ephemCalc/meeus.c b/src/ephemCalc/meeus.c
@@ -1,7 +1,7 @@
 // meeus.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/ephemCalc/meeus.h b/src/ephemCalc/meeus.h
@@ -1,7 +1,7 @@
 // meeus.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/ephemCalc/orbitalElements.c b/src/ephemCalc/orbitalElements.c
@@ -1,7 +1,7 @@
 // orbitalElements.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //
@@ -483,22 +483,22 @@ void orbitalElements_asteroids_readAsciiData() {
                                                           ((int) floor(tmp)) % 100, 0, 0, 0, &i, temp_err_string);
         line_ptr = next_word(line_ptr);
 
-        // radians; J2000.0
+        // Read mean anomaly -- radians; J2000.0
         asteroid_database[n].meanAnomaly = get_float(line_ptr, NULL) * M_PI / 180;
         line_ptr = next_word(line_ptr);
-        // radians; J2000.0
+        // Read argument of perihelion -- radians; J2000.0
         asteroid_database[n].argumentPerihelion = get_float(line_ptr, NULL) * M_PI / 180;
         line_ptr = next_word(line_ptr);
-        // radians; J2000.0
+        // Read longitude of ascending node -- radians; J2000.0
         asteroid_database[n].longAscNode = get_float(line_ptr, NULL) * M_PI / 180;
         line_ptr = next_word(line_ptr);
-        // radians; J2000.0
+        // Read inclination of orbit -- radians; J2000.0
         asteroid_database[n].inclination = get_float(line_ptr, NULL) * M_PI / 180;
         line_ptr = next_word(line_ptr);
-        // dimensionless
+        // Read eccentricity of orbit -- dimensionless
         asteroid_database[n].eccentricity = get_float(line_ptr, NULL);
         line_ptr = next_word(line_ptr);
-        // AU
+        // Read semi-major axis of orbit -- AU
         asteroid_database[n].semiMajorAxis = get_float(line_ptr, NULL);
     }
     fclose(input);
@@ -621,7 +621,6 @@ void orbitalElements_comets_readAsciiData() {
         while ((line[j] > ' ') && (k < 23)) comet_database[comet_count].name2[k++] = line[j++];
         comet_database[comet_count].name2[k] = '\0';
 
-
         // Read perihelion distance
         perihelion_dist = get_float(line + 31, NULL);
 
@@ -706,7 +705,7 @@ void orbitalElements_comets_readAsciiData() {
     OrbitalElements_DumpBinaryData("dcfbinary.cmt", comet_database, &comet_database_offset,
                                    comet_count, comet_secure_count);
 
-    // Make table indicating that we have loaded all of the orbital elements in this table
+    // Make table indicating that we have loaded all the orbital elements in this table
     comet_database_items_loaded = (unsigned char *) lt_malloc(comet_count * sizeof(unsigned char));
     memset(comet_database_items_loaded, 1, comet_count);
 
@@ -827,8 +826,10 @@ void orbitalElements_computeXYZ(int body_id, double jd, double *x, double *y, do
 
     double v, r;
 
+    // const double epsilon = (23.4393 - 3.563E-7 * (jd - 2451544.5)) * M_PI / 180;
+
+    // Case 1: Object is a planet
     if (body_id < 1000000) {
-        // Case 1: Object is a planet
         // Planets occupy body numbers 1-19
         const int index = body_id;
 
@@ -844,8 +845,8 @@ void orbitalElements_computeXYZ(int body_id, double jd, double *x, double *y, do
         orbital_elements = orbitalElements_planets_fetch(index);
     }
 
-    else if (body_id < 2000000) {
         // Case 2: Object is an asteroid
+    else if (body_id < 2000000) {
         // Asteroids occupy body numbers 1e6 - 2e6
         const int index = body_id - 1000000;
 
@@ -861,8 +862,8 @@ void orbitalElements_computeXYZ(int body_id, double jd, double *x, double *y, do
         orbital_elements = orbitalElements_asteroids_fetch(index);
     }
 
-    else {
         // Case 3: Object is a comet
+    else {
         // Comets occupy body numbers 2e6 - 3e6
         const int index = body_id - 2000000;
 
@@ -1054,6 +1055,10 @@ void orbitalElements_computeEphemeris(settings *i, int bodyId, double jd, double
     // Boolean flags indicating whether this is the Earth, Sun or Moon (which need special treatment)
     int is_moon = 0, is_earth = 0, is_sun = 0;
 
+    double EMX_future, EMY_future, EMZ_future; // Position of the Earth-Moon centre of mass
+    double moon_pos_x_future, moon_pos_y_future, moon_pos_z_future;
+    double earth_pos_x_future, earth_pos_y_future, earth_pos_z_future;
+
     // Earth: Need to convert from Earth/Moon barycentre to geocentre
     if (bodyId == 19) {
         bodyId = 2;
@@ -1071,14 +1076,13 @@ void orbitalElements_computeEphemeris(settings *i, int bodyId, double jd, double
     }
 
     // Look up position of the Earth at this JD, so that we can convert XYZ coordinates relative to Sun into
-    // RA and Dec as observed from the Earth
-    const double earth_mass = 5.9736e24;
-    const double moon_mass = 7.3477e22;
+    // RA and Dec as observed from the Earth.
+    // Below are values of GM3 and GMM from DE405. See
+    // <https://web.archive.org/web/20120220062549/http://iau-comm4.jpl.nasa.gov/de405iom/de405iom.pdf>
+    const double earth_mass = 0.8887692390113509e-9;
+    const double moon_mass = 0.1093189565989898e-10;
     const double moon_earth_mass_ratio = moon_mass / (moon_mass + earth_mass);
 
-    // Look up the Sun's position
-    jpl_computeXYZ(10, jd, &sun_pos_x, &sun_pos_y, &sun_pos_z);
-
     // Look up the Earth-Moon centre of mass position
     jpl_computeXYZ(2, jd, &EMX, &EMY, &EMZ);
 
@@ -1090,45 +1094,98 @@ void orbitalElements_computeEphemeris(settings *i, int bodyId, double jd, double
     earth_pos_y = EMY - moon_earth_mass_ratio * moon_pos_y;
     earth_pos_z = EMZ - moon_earth_mass_ratio * moon_pos_z;
 
+    // Look up the Sun's position, taking light travel time into account
+    {
+        jpl_computeXYZ(10, jd, &sun_pos_x, &sun_pos_y, &sun_pos_z);
+
+        // Calculate light travel time
+        const double distance = gsl_hypot3(sun_pos_x - earth_pos_x,
+                                           sun_pos_y - earth_pos_y,
+                                           sun_pos_z - earth_pos_z);  // AU
+        const double light_travel_time = distance * GSL_CONST_MKSA_ASTRONOMICAL_UNIT / GSL_CONST_MKSA_SPEED_OF_LIGHT;
+
+        // Look up position of requested object at the time the light left the object
+        jpl_computeXYZ(10, jd - light_travel_time / 86400, &sun_pos_x, &sun_pos_y, &sun_pos_z);
+    }
+
+    // If the user's query was about the Earth, we already know its position
     if (is_earth) {
-        // If the user's query was about the Earth, we already know its position
         *x = earth_pos_x;
         *y = earth_pos_y;
         *z = earth_pos_z;
     }
+
+        // If the user's query was about the Sun, we already know that position too
     else if (is_sun) {
-        // If the user's query was about the Sun, we already know that position too!
         *x = sun_pos_x;
         *y = sun_pos_y;
         *z = sun_pos_z;
     }
+
+        // If the user's query was about the Moon, we already know that position too
     else if (is_moon) {
-        // If the user's query was about the Moon, we already know that position too!
-        *x = moon_pos_x;
-        *y = moon_pos_y;
-        *z = moon_pos_z;
+        *x = moon_pos_x + earth_pos_x;
+        *y = moon_pos_y + earth_pos_y;
+        *z = moon_pos_z + earth_pos_z;
     }
+
+        // Otherwise we need to use the orbital elements for the particular object the user was looking for,
+        // taking light travel time into account
     else {
-        // ... otherwise we need to use the orbital elements for the particular object the user was looking for
+        // Calculate position of requested object at specified time
         orbitalElements_computeXYZ(bodyId, jd, x, y, z);
 
-        // Orbital elements give the position of objects relative to the Sun.
-        // Switch coordinates to be relative to the solar system barycentre
-        *x += sun_pos_x;
-        *y += sun_pos_y;
-        *z += sun_pos_z;
+        // Calculate light travel time
+        const double distance = gsl_hypot3(*x - earth_pos_x, *y - earth_pos_y, *z - earth_pos_z);  // AU
+        const double light_travel_time = distance * GSL_CONST_MKSA_ASTRONOMICAL_UNIT / GSL_CONST_MKSA_SPEED_OF_LIGHT;
+
+        // Look up position of requested object at the time the light left the object
+        orbitalElements_computeXYZ(bodyId, jd - light_travel_time / 86400, x, y, z);
     }
 
-    // If the body was the Moon, don't forget the add the position of the Earth-Moon centre of mass
-    if (is_moon) {
-        *x += earth_pos_x;
-        *y += earth_pos_y;
-        *z += earth_pos_z;
+    // Look up the Earth-Moon centre of mass position, a short time in the future
+    // We use this to calculate the Earth's velocity vector, which is needed to correct for aberration
+    // (see eqn 7.119 of the Explanatory Supplement)
+    const double eb_dot_timestep = 1e-6; // days
+    const double eb_dot_timestep_sec = eb_dot_timestep * 86400;
+    jpl_computeXYZ(2, jd + eb_dot_timestep, &EMX_future, &EMY_future, &EMZ_future);
+    jpl_computeXYZ(9, jd + eb_dot_timestep, &moon_pos_x_future, &moon_pos_y_future, &moon_pos_z_future);
+    earth_pos_x_future = EMX_future - moon_earth_mass_ratio * moon_pos_x_future;
+    earth_pos_y_future = EMY_future - moon_earth_mass_ratio * moon_pos_y_future;
+    earth_pos_z_future = EMZ_future - moon_earth_mass_ratio * moon_pos_z_future;
+
+    // Equation (7.118) of the Explanatory Supplement - correct for aberration
+    if (!is_earth) {
+        const double u1[3] = {
+                *x - earth_pos_x,
+                *y - earth_pos_y,
+                *z - earth_pos_z
+        };
+        const double u1_mag = gsl_hypot3(u1[0], u1[1], u1[2]);
+        const double u[3] = {u1[0] / u1_mag, u1[1] / u1_mag, u1[2] / u1_mag};
+        const double eb_dot[3] = {
+                earth_pos_x_future - earth_pos_x,
+                earth_pos_y_future - earth_pos_y,
+                earth_pos_z_future - earth_pos_z
+        };
+
+        // Speed of light in AU per time step
+        const double c = GSL_CONST_MKSA_SPEED_OF_LIGHT / GSL_CONST_MKSA_ASTRONOMICAL_UNIT * eb_dot_timestep_sec;
+        const double V[3] = {eb_dot[0] / c, eb_dot[1] / c, eb_dot[2] / c};
+        const double V_mag = gsl_hypot3(V[0], V[1], V[2]);
+        const double beta = sqrt(1 - gsl_pow_2(V_mag));
+        const double f1 = u[0] * V[0] + u[1] * V[1] + u[2] * V[2];
+        const double f2 = 1 + f1 / (1 + beta);
+
+        // Correct for aberration
+        *x = earth_pos_x + (beta * u1[0] + f2 * u1_mag * V[0]) / (1 + f1);
+        *y = earth_pos_y + (beta * u1[1] + f2 * u1_mag * V[1]) / (1 + f1);
+        *z = earth_pos_z + (beta * u1[2] + f2 * u1_mag * V[2]) / (1 + f1);
     }
 
     // Populate other quantities, like the brightness, RA and Dec of the object, based on its XYZ position
-    magnitudeEstimate(bodyId, *x, *y, *z, earth_pos_x, earth_pos_y, earth_pos_z,
-                      sun_pos_x, sun_pos_y, sun_pos_z, ra, dec, mag, phase, angSize, phySize,
+    magnitudeEstimate(bodyId, *x, *y, *z, earth_pos_x, earth_pos_y, earth_pos_z, sun_pos_x, sun_pos_y, sun_pos_z, ra,
+                      dec, mag, phase, angSize, phySize,
                       albedo, sunDist, earthDist, sunAngDist, theta_eso, eclipticLongitude, eclipticLatitude,
                       eclipticDistance, i);
 }
diff --git a/src/ephemCalc/orbitalElements.h b/src/ephemCalc/orbitalElements.h
@@ -1,7 +1,7 @@
 // orbitalElements.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/listTools/ltDict.c b/src/listTools/ltDict.c
@@ -1,7 +1,7 @@
 // ltDict.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/listTools/ltDict.h b/src/listTools/ltDict.h
@@ -1,7 +1,7 @@
 // ltDict.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/listTools/ltList.c b/src/listTools/ltList.c
@@ -1,7 +1,7 @@
 // ltList.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/listTools/ltList.h b/src/listTools/ltList.h
@@ -1,7 +1,7 @@
 // ltList.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/listTools/ltMemory.c b/src/listTools/ltMemory.c
@@ -1,7 +1,7 @@
 // ltMemory.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/listTools/ltMemory.h b/src/listTools/ltMemory.h
@@ -1,7 +1,7 @@
 // ltMemory.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/listTools/ltStringProc.c b/src/listTools/ltStringProc.c
@@ -1,7 +1,7 @@
 // ltStringProc.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/listTools/ltStringProc.h b/src/listTools/ltStringProc.h
@@ -1,7 +1,7 @@
 // ltStringProc.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/main.c b/src/main.c
@@ -1,7 +1,7 @@
 // main.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/mathsTools/julianDate.c b/src/mathsTools/julianDate.c
@@ -1,7 +1,7 @@
 // julianDate.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/mathsTools/julianDate.h b/src/mathsTools/julianDate.h
@@ -1,7 +1,7 @@
 // julianDate.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/mathsTools/precess_equinoxes.c b/src/mathsTools/precess_equinoxes.c
@@ -1,7 +1,7 @@
 // precess_equinoxes.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/mathsTools/precess_equinoxes.h b/src/mathsTools/precess_equinoxes.h
@@ -1,7 +1,7 @@
 // precess_equinoxes.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/mathsTools/sphericalAst.c b/src/mathsTools/sphericalAst.c
@@ -1,7 +1,7 @@
 // sphericalAst.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/mathsTools/sphericalAst.h b/src/mathsTools/sphericalAst.h
@@ -1,7 +1,7 @@
 // sphericalAst.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/settings/settings.c b/src/settings/settings.c
@@ -1,7 +1,7 @@
 // settings.c
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //

diff --git a/src/settings/settings.h b/src/settings/settings.h
@@ -1,7 +1,7 @@
 // settings.h
 // 
 // -------------------------------------------------
-// Copyright 2015-2022 Dominic Ford
+// Copyright 2015-2024 Dominic Ford
 //
 // This file is part of EphemerisCompute.
 //