From 813bbf1c8e84f55cb9cdff9d79b5302d6767aa12 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Fri, 12 Nov 2021 15:14:56 -0500 Subject: [PATCH] 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 --- generate/codegen.c | 35 ++++++++------- generate/gravsim/gravsim_test.c | 10 ++--- generate/gravsim/pluto_gravsim.h | 28 ++++++++++++ generate/template/astronomy.c | 65 +++++++++++----------------- source/c/astronomy.c | 74 ++++++++++++++------------------ 5 files changed, 108 insertions(+), 104 deletions(-) create mode 100644 generate/gravsim/pluto_gravsim.h diff --git a/generate/codegen.c b/generate/codegen.c index ebc13f13..c352f94f 100644 --- 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 */ diff --git a/generate/gravsim/gravsim_test.c b/generate/gravsim/gravsim_test.c index 5402e545..bac565a7 100644 --- 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; diff --git a/generate/gravsim/pluto_gravsim.h b/generate/gravsim/pluto_gravsim.h new file mode 100644 index 00000000..016646e8 --- /dev/null +++ b/generate/gravsim/pluto_gravsim.h @@ -0,0 +1,28 @@ +/* + pluto_gravsim.h - Don Cross + 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 diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 0e16551d..977c9894 100644 --- 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]; diff --git a/source/c/astronomy.c b/source/c/astronomy.c index eeaed6f0..221296b2 100644 --- 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];