Pluto gravity sim: refactor constants for sharing.

Reworked the Pluto gravity sim constants so they are defined in one place: a new header file gravsim/pluto_gravsim.h. Then the code generator writes the #defines to the C code, instead of having two independent versions of the same constants. I will continue down the road of having a single-source-of-truth for these constants across all 4 supported languages. Also, confusingly, I had one constant called PLUTO_DT in codegen.c that was called PLUTO_TIME_STEP in astronomy.c. Also, astronomy.c had a different constant PLUTO_DT that didn't mean the same thing. I reworked the naming to be consistent in all places. I already had a TopPosition() function that knows how to calculate exact equatorial coordinates, so I eliminated the redundant logic from gravsim_test.c
2026-05-19 14:27:52 -04:00 · 2021-11-12 15:14:56 -05:00
parent dec8fd6f24
commit 813bbf1c8e
5 changed files with 108 additions and 104 deletions
--- a/generate/codegen.c
+++ b/generate/codegen.c
@@ -30,6 +30,7 @@
 #include "codegen.h"
 #include "top2013.h"
 #include "ephfile.h"
+#include "gravsim/pluto_gravsim.h"

 #define CG_MAX_LINE_LENGTH  200
 #define MAX_DATA_PER_LINE    20
@@ -1207,10 +1208,16 @@ fail:
    if (infile) fclose(infile);
    return error;
 }
-#define PLUTO_NUM_STATES    41
-#define PLUTO_TT1           (-730000.0)     /* 0001-04-30T12:00:00.000Z */
-#define PLUTO_TT2           (+730000.0)     /* 3998-09-03T12:00:00.000Z */
-#define PLUTO_DT            ((PLUTO_TT2 - PLUTO_TT1) / (PLUTO_NUM_STATES - 1))
+
+
+static int PlutoConstants_C(cg_context_t *context)
+{
+    fprintf(context->outfile, "#define PLUTO_NUM_STATES  %d\n", PLUTO_NUM_STATES);
+    fprintf(context->outfile, "#define PLUTO_TIME_STEP   %d\n", PLUTO_TIME_STEP);
+    fprintf(context->outfile, "#define PLUTO_DT          %d\n\n", PLUTO_DT);
+    fprintf(context->outfile, "#define PLUTO_NSTEPS      %d\n\n", PLUTO_NSTEPS);
+    return 0;
+}


 static int PlutoStateTable_C(cg_context_t *context, const top_model_t *model)
@@ -1220,14 +1227,12 @@ static int PlutoStateTable_C(cg_context_t *context, const top_model_t *model)
    top_rectangular_t equ;
    int i;

-    fprintf(context->outfile, "#define PLUTO_NUM_STATES  %d\n", PLUTO_NUM_STATES);
-    fprintf(context->outfile, "#define PLUTO_TIME_STEP   %0.0lf\n\n", PLUTO_DT);
    fprintf(context->outfile, "static const body_state_t PlutoStateTable[] =\n");
    fprintf(context->outfile, "{\n");

    for (i=0; i < PLUTO_NUM_STATES; ++i)
    {
-        tt = i*PLUTO_DT + PLUTO_TT1;
+        tt = i*PLUTO_TIME_STEP + PLUTO_TT1;
        CHECK(TopPosition(model, tt, &equ));

        fprintf(context->outfile,
@@ -1252,13 +1257,13 @@ static int PlutoStateTable_CSharp(cg_context_t *context, const top_model_t *mode
    int i;

    fprintf(context->outfile, "        private const int PLUTO_NUM_STATES = %d;\n", PLUTO_NUM_STATES);
-    fprintf(context->outfile, "        private const int PLUTO_TIME_STEP  = %0.0lf;\n\n", PLUTO_DT);
+    fprintf(context->outfile, "        private const int PLUTO_TIME_STEP  = %d;\n\n", PLUTO_TIME_STEP);
    fprintf(context->outfile, "        private static readonly body_state_t[] PlutoStateTable = new body_state_t[]\n");
    fprintf(context->outfile, "        {\n");

    for (i=0; i < PLUTO_NUM_STATES; ++i)
    {
-        tt = i*PLUTO_DT + PLUTO_TT1;
+        tt = i*PLUTO_TIME_STEP + PLUTO_TT1;
        CHECK(TopPosition(model, tt, &equ));

        fprintf(context->outfile,
@@ -1283,12 +1288,12 @@ static int PlutoStateTable_JS(cg_context_t *context, const top_model_t *model)
    int i;

    fprintf(context->outfile, "const PLUTO_NUM_STATES = %d;\n", PLUTO_NUM_STATES);
-    fprintf(context->outfile, "const PLUTO_TIME_STEP  = %0.0lf;\n\n", PLUTO_DT);
+    fprintf(context->outfile, "const PLUTO_TIME_STEP  = %d;\n\n", PLUTO_TIME_STEP);
    fprintf(context->outfile, "const PlutoStateTable: BodyStateTableEntry[] = [\n");

    for (i=0; i < PLUTO_NUM_STATES; ++i)
    {
-        tt = i*PLUTO_DT + PLUTO_TT1;
+        tt = i*PLUTO_TIME_STEP + PLUTO_TT1;
        CHECK(TopPosition(model, tt, &equ));

        fprintf(context->outfile,
@@ -1313,12 +1318,12 @@ static int PlutoStateTable_Python(cg_context_t *context, const top_model_t *mode
    int i;

    fprintf(context->outfile, "_PLUTO_NUM_STATES = %d\n", PLUTO_NUM_STATES);
-    fprintf(context->outfile, "_PLUTO_TIME_STEP  = %0.0lf\n\n", PLUTO_DT);
+    fprintf(context->outfile, "_PLUTO_TIME_STEP  = %d\n\n", PLUTO_TIME_STEP);
    fprintf(context->outfile, "_PlutoStateTable = [\n");

    for (i=0; i < PLUTO_NUM_STATES; ++i)
    {
-        tt = i*PLUTO_DT + PLUTO_TT1;
+        tt = i*PLUTO_TIME_STEP + PLUTO_TT1;
        CHECK(TopPosition(model, tt, &equ));

        fprintf(context->outfile,
@@ -1341,9 +1346,6 @@ static int PlutoStateTable(cg_context_t *context)
    top_model_t model;
    TopInitModel(&model);

-    if (PLUTO_DT != round(PLUTO_DT))
-        CHECK(LogError(context, "PlutoStateTable: PLUTO_DT = %lf is not an integer.\n", PLUTO_DT));
-
    CHECK(TopLoadModel(&model, "TOP2013.dat", 9));

    switch (context->language)
@@ -1858,6 +1860,7 @@ static const cg_directive_entry DirectiveTable[] =
    { "IAU_DATA",           OptIauData          },
    { "ADDSOL",             OptAddSol           },
    { "CONSTEL",            ConstellationData   },
+    { "C_PLUTO_CONST",      PlutoConstants_C    },
    { "PLUTO_TABLE",        PlutoStateTable     },
    { "JUPITER_MOONS",      JupiterMoons        },
    { NULL, NULL }  /* Marks end of list */
--- a/generate/gravsim/gravsim_test.c
+++ b/generate/gravsim/gravsim_test.c
@@ -11,6 +11,7 @@
 #include "astronomy.h"
 #include "top2013.h"
 #include "codegen.h"        // for CHECK macro
+#include "pluto_gravsim.h"

 static int DebugMode;

@@ -18,8 +19,7 @@ int TestPluto(double tt, const top_model_t *model, double *arcmin)
 {
    astro_time_t time;
    astro_vector_t pos;
-    top_elliptical_t ellip;
-    top_rectangular_t ecl, equ;
+    top_rectangular_t equ;
    double dx, dy, dz;
    int error;

@@ -38,9 +38,7 @@ int TestPluto(double tt, const top_model_t *model, double *arcmin)
    if (DebugMode) printf("HelioVector : pos=(%23.16lf, %23.16lf, %23.16lf)\n", pos.x, pos.y, pos.z);

    /* Compare with untruncated TOP2013 model. */
-    CHECK(TopCalcElliptical(model, tt, &ellip));
-    CHECK(TopEcliptic(model->planet, &ellip, &ecl));
-    CHECK(TopEquatorial(&ecl, &equ));
+    CHECK(TopPosition(model, tt, &equ));
    if (DebugMode) printf("TOP2013     : pos=(%23.16lf, %23.16lf, %23.16lf)\n", equ.x, equ.y, equ.z);

    /* Calculate relative error and scale in arcminute units (as seen from the Sun). */
@@ -59,8 +57,6 @@ fail:
 int main(int argc, const char *argv[])
 {
    top_model_t model;
-    const double PLUTO_TIME_STEP = 36500.0;     /* must keep in sync with same-name symbol in astronomy.c */
-    const double PLUTO_DT = 250.0;              /* must keep in sync with same-name symbol in astronomy.c */
    double arcmin, tt, dev;
    int error, n, count;

--- a/generate/gravsim/pluto_gravsim.h
+++ b/generate/gravsim/pluto_gravsim.h
@@ -0,0 +1,28 @@
+/*
+    pluto_gravsim.h  -  Don Cross <cosinekitty@gmail.com>
+    https://github.com/cosinekitty/astronomy
+
+    Constants that tune the gravitational simulation of Pluto's orbit.
+*/
+
+#ifndef __ASTRONOMY_PLUTO_GRAVSIM_H
+#define __ASTRONOMY_PLUTO_GRAVSIM_H
+
+#define PLUTO_NUM_STATES    41
+#define PLUTO_TT1           (-730000)     /* 0001-04-30T12:00:00.000Z */
+#define PLUTO_TT2           (+730000)     /* 3998-09-03T12:00:00.000Z */
+
+#if ((PLUTO_TT2 - PLUTO_TT1) % (PLUTO_NUM_STATES - 1)) != 0
+    #error PLUTO_TIME_STEP ratio must be an integer.
+#endif
+
+#define PLUTO_TIME_STEP     ((PLUTO_TT2 - PLUTO_TT1) / (PLUTO_NUM_STATES - 1))
+#define PLUTO_DT            250
+
+#if PLUTO_TIME_STEP % PLUTO_DT != 0
+    #error Invalid combination of Pluto time step, time increment.
+#endif
+
+#define PLUTO_NSTEPS    ((PLUTO_TIME_STEP / PLUTO_DT) + 1)
+
+#endif
--- a/generate/template/astronomy.c
+++ b/generate/template/astronomy.c
@@ -39,6 +39,8 @@ extern "C" {
 /** @cond DOXYGEN_SKIP */
 #define PI      3.14159265358979323846

+$ASTRO_C_PLUTO_CONST()
+
 typedef enum
 {
    FROM_2000,
@@ -54,6 +56,30 @@ typedef struct
 }
 terse_vector_t;

+typedef struct
+{
+    double          tt;   /* J2000 terrestrial time [days] */
+    terse_vector_t  r;    /* position [au] */
+    terse_vector_t  v;    /* velocity [au/day] */
+    terse_vector_t  a;    /* acceleration [au/day^2] */
+}
+body_grav_calc_t;
+
+typedef struct
+{
+    body_grav_calc_t   step[PLUTO_NSTEPS];
+}
+body_segment_t;
+
+typedef struct
+{
+    double          tt;  /* Terrestrial Time in J2000 days */
+    terse_vector_t  r;   /* position [au] */
+    terse_vector_t  v;   /* velocity [au/day] */
+}
+body_state_t;
+/** @endcond */
+
 static const terse_vector_t VecZero = { 0.0, 0.0, 0.0 };

 static terse_vector_t VecAdd(terse_vector_t a, terse_vector_t b)
@@ -142,15 +168,6 @@ static astro_vector_t PublicVec(astro_time_t time, terse_vector_t terse)
    return vector;
 }

-typedef struct
-{
-    double          tt;  /* Terrestrial Time in J2000 days */
-    terse_vector_t  r;   /* position [au] */
-    terse_vector_t  v;   /* velocity [au/day] */
-}
-body_state_t;
-/** @endcond */
-
 static const double DAYS_PER_TROPICAL_YEAR = 365.24217;
 static const double ASEC360 = 1296000.0;
 static const double ASEC2RAD = 4.848136811095359935899141e-6;
@@ -2308,16 +2325,6 @@ static astro_vector_t CalcSolarSystemBarycenter(astro_time_t time)

 /*------------------ begin Pluto integrator ------------------*/

-/** @cond DOXYGEN_SKIP */
-typedef struct
-{
-    double          tt;   /* J2000 terrestrial time [days] */
-    terse_vector_t  r;    /* position [au] */
-    terse_vector_t  v;    /* velocity [au/day] */
-    terse_vector_t  a;    /* acceleration [au/day^2] */
-} body_grav_calc_t;
-/** @endcond */
-
 $ASTRO_PLUTO_TABLE();


@@ -2498,26 +2505,6 @@ body_grav_calc_t GravSim(           /* out: [pos, vel, acc] of the simulated bod
    return calc2;
 }

-/* Allow customization of the time step inside a gravsim segment. */
-#ifndef PLUTO_DT
-#define PLUTO_DT 250
-#endif
-
-#if PLUTO_TIME_STEP % PLUTO_DT != 0
-    #error Invalid combination of Pluto time step, time increment.
-#endif
-
-#define PLUTO_NSTEPS    ((PLUTO_TIME_STEP / PLUTO_DT) + 1)
-
-/** @cond DOXYGEN_SKIP */
-typedef struct
-{
-    body_grav_calc_t   step[PLUTO_NSTEPS];
-}
-body_segment_t;
-/** @endcond */
-
-
 /* FIXFIXFIX - Using a global is not thread-safe. Either add thread-locks or change API to accept a cache pointer. */
 static body_segment_t *pluto_cache[PLUTO_NUM_STATES-1];

--- a/source/c/astronomy.c
+++ b/source/c/astronomy.c
@@ -39,6 +39,14 @@ extern "C" {
 /** @cond DOXYGEN_SKIP */
 #define PI      3.14159265358979323846

+#define PLUTO_NUM_STATES  41
+#define PLUTO_TIME_STEP   36500
+#define PLUTO_DT          250
+
+#define PLUTO_NSTEPS      147
+
+
+
 typedef enum
 {
    FROM_2000,
@@ -54,6 +62,30 @@ typedef struct
 }
 terse_vector_t;

+typedef struct
+{
+    double          tt;   /* J2000 terrestrial time [days] */
+    terse_vector_t  r;    /* position [au] */
+    terse_vector_t  v;    /* velocity [au/day] */
+    terse_vector_t  a;    /* acceleration [au/day^2] */
+}
+body_grav_calc_t;
+
+typedef struct
+{
+    body_grav_calc_t   step[PLUTO_NSTEPS];
+}
+body_segment_t;
+
+typedef struct
+{
+    double          tt;  /* Terrestrial Time in J2000 days */
+    terse_vector_t  r;   /* position [au] */
+    terse_vector_t  v;   /* velocity [au/day] */
+}
+body_state_t;
+/** @endcond */
+
 static const terse_vector_t VecZero = { 0.0, 0.0, 0.0 };

 static terse_vector_t VecAdd(terse_vector_t a, terse_vector_t b)
@@ -142,15 +174,6 @@ static astro_vector_t PublicVec(astro_time_t time, terse_vector_t terse)
    return vector;
 }

-typedef struct
-{
-    double          tt;  /* Terrestrial Time in J2000 days */
-    terse_vector_t  r;   /* position [au] */
-    terse_vector_t  v;   /* velocity [au/day] */
-}
-body_state_t;
-/** @endcond */
-
 static const double DAYS_PER_TROPICAL_YEAR = 365.24217;
 static const double ASEC360 = 1296000.0;
 static const double ASEC2RAD = 4.848136811095359935899141e-6;
@@ -3314,19 +3337,6 @@ static astro_vector_t CalcSolarSystemBarycenter(astro_time_t time)

 /*------------------ begin Pluto integrator ------------------*/

-/** @cond DOXYGEN_SKIP */
-typedef struct
-{
-    double          tt;   /* J2000 terrestrial time [days] */
-    terse_vector_t  r;    /* position [au] */
-    terse_vector_t  v;    /* velocity [au/day] */
-    terse_vector_t  a;    /* acceleration [au/day^2] */
-} body_grav_calc_t;
-/** @endcond */
-
-#define PLUTO_NUM_STATES  41
-#define PLUTO_TIME_STEP   36500
-
 static const body_state_t PlutoStateTable[] =
 {
    {  -730000.0, {-26.1182072321076, -14.3761681778250,   3.3844025152995}, { 1.6339372163656e-03, -2.7861699588508e-03, -1.3585880229445e-03} }
@@ -3550,26 +3560,6 @@ body_grav_calc_t GravSim(           /* out: [pos, vel, acc] of the simulated bod
    return calc2;
 }

-/* Allow customization of the time step inside a gravsim segment. */
-#ifndef PLUTO_DT
-#define PLUTO_DT 250
-#endif
-
-#if PLUTO_TIME_STEP % PLUTO_DT != 0
-    #error Invalid combination of Pluto time step, time increment.
-#endif
-
-#define PLUTO_NSTEPS    ((PLUTO_TIME_STEP / PLUTO_DT) + 1)
-
-/** @cond DOXYGEN_SKIP */
-typedef struct
-{
-    body_grav_calc_t   step[PLUTO_NSTEPS];
-}
-body_segment_t;
-/** @endcond */
-
-
 /* FIXFIXFIX - Using a global is not thread-safe. Either add thread-locks or change API to accept a cache pointer. */
 static body_segment_t *pluto_cache[PLUTO_NUM_STATES-1];