Include the gravsim check in the unit tests.

I have gravsim_test.c to the point where it calculates a standard deviation of error between TOP2013 and Astronomy Engine for calculating the position of Pluto over 10 worst-case samples. My baseline is now 0.205303 arcminutes of heliocentric position error. For Runge-Kutta (or some other method) to be an improvement, it has to beat that score without incurring significant extra work or larger memory consumption.
2026-05-19 14:27:52 -04:00 · 2021-11-10 21:00:38 -05:00
parent 1bebad6602
commit cc6a32bb98
4 changed files with 168 additions and 31 deletions
--- a/generate/gravsim/gravsim_test.c
+++ b/generate/gravsim/gravsim_test.c
@@ -6,19 +6,26 @@
 */

 #include <stdio.h>
+#include <string.h>
 #include <math.h>
 #include "astronomy.h"
 #include "top2013.h"
+#include "codegen.h"        // for CHECK macro

-int TestPluto(double tt, const top_model_t *model)
+static int DebugMode;
+
+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;
-    double dx, dy, dz, arcmin;
+    double dx, dy, dz;
+    int error;

-    printf("tt = %0.1lf\n", tt);
+    *arcmin = 999.0;    /* bogus value in case of failure */
+
+    if (DebugMode) printf("tt = %0.1lf\n", tt);

    time = Astronomy_TerrestrialTime(tt);
    pos = Astronomy_HelioVector(BODY_PLUTO, time);
@@ -28,47 +35,70 @@ int TestPluto(double tt, const top_model_t *model)
        return 1;
    }

-    printf("HelioVector : pos=(%23.16lf, %23.16lf, %23.16lf)\n", pos.x, pos.y, pos.z);
+    if (DebugMode) printf("HelioVector : pos=(%23.16lf, %23.16lf, %23.16lf)\n", pos.x, pos.y, pos.z);

    /* Compare with untruncated TOP2013 model. */
-    if (TopCalcElliptical(model, tt, &ellip)) return 1;
-    if (TopEcliptic(model->planet, &ellip, &ecl)) return 1;
-    if (TopEquatorial(&ecl, &equ)) return 1;
-    printf("TOP2013     : pos=(%23.16lf, %23.16lf, %23.16lf)\n", equ.x, equ.y, equ.z);
+    CHECK(TopCalcElliptical(model, tt, &ellip));
+    CHECK(TopEcliptic(model->planet, &ellip, &ecl));
+    CHECK(TopEquatorial(&ecl, &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). */
    dx = equ.x - pos.x;
    dy = equ.y - pos.y;
    dz = equ.z - pos.z;
-    arcmin = (RAD2DEG * 60.0) * sqrt((dx*dx + dy*dy + dz*dz) / (equ.x*equ.x + equ.y*equ.y + equ.z*equ.z));
-    printf("arcmin_error = %0.6lf\n", arcmin);
-    printf("\n");
+    *arcmin = (RAD2DEG * 60.0) * sqrt((dx*dx + dy*dy + dz*dz) / (equ.x*equ.x + equ.y*equ.y + equ.z*equ.z));
+    if (DebugMode) printf("arcmin_error = %0.6lf\n", *arcmin);
+    if (DebugMode) printf("\n");

-    return 0;
+    error = 0;
+fail:
+    return error;
 }

-int main(void)
+int main(int argc, const char *argv[])
 {
    top_model_t model;
-    const double delta_tt = 18250.0;
-    int error;
+    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;
+
+    DebugMode = (argc > 1) && !strcmp(argv[1], "-v");

    TopInitModel(&model);
-    error = TopLoadModel(&model, "../TOP2013.dat", 9);
-    if (error) goto fail;
+    CHECK(TopLoadModel(&model, "../TOP2013.dat", 9));

-    /*
-        The PlutoStateTable has exact values at tt = { ... , -36500, 0, +36500, ... }.
-        I want to exercise the calculator at the midpoints tt = { ..., -18250, +18250, ... }.
-        We also exercise at the exact points where we coincide with TOP2013, to verify we get
-        exact match.
-    */
-    for (int n = -5; n <= 5; ++n)
+    /* Verify that we are in sync at exact checkpoints. */
+    if (DebugMode) printf("(Expect exact matches...)\n");
+    for (n = -5; n < 5; ++n)
    {
-        double tt = n * delta_tt;
-        if (TestPluto(tt, &model)) goto fail;
+        tt = n * PLUTO_TIME_STEP;
+        CHECK(TestPluto(tt, &model, &arcmin));
+        if (arcmin > 1.0e-6)
+            FAIL("EXCESSIVE ERROR!\n");
    }

+    /* Now check worst-case at near-midpoints. */
+    if (DebugMode) printf("(Expect simulation differences...)\n");
+    count = 0;
+    dev = 0.0;
+    for (n = -5; n < 5; ++n)
+    {
+        /* Calculate a worst case value: the midpoint of a step, near the midpoint of a segment. */
+        /* Each segment is 36500 days, and each step is 250 days. */
+        tt = ((n + 0.5) * PLUTO_TIME_STEP) + (PLUTO_DT / 2.0);
+        CHECK(TestPluto(tt, &model, &arcmin));
+        dev += arcmin * arcmin;
+        ++count;
+    }
+
+    dev = sqrt(dev / count);
+    printf("gravsim_test.c: Pluto score = %0.6lf arcmin, over %d data points.\n", dev, count);
+    if (dev > 0.205303)
+        FAIL("gravsim_test.c: EXCESSIVE ERROR\n");
+
+    error = 0;
 fail:
    TopFreeModel(&model);
    return error;
--- a/generate/gravsim/original.txt
+++ b/generate/gravsim/original.txt
@@ -0,0 +1,103 @@
+(Expect exact matches...)
+tt = -182500.0
+HelioVector : pos=(   -14.1723844533378998,    -26.0054690135835997,     -3.8387026446525998)
+TOP2013     : pos=(   -14.1723844533378838,    -26.0054690135836388,     -3.8387026446526260)
+arcmin_error = 0.000000
+
+tt = -146000.0
+HelioVector : pos=(    40.9468572586402999,     25.9049735920209017,     -4.2563362404987997)
+TOP2013     : pos=(    40.9468572586402928,     25.9049735920209017,     -4.2563362404987810)
+arcmin_error = 0.000000
+
+tt = -109500.0
+HelioVector : pos=(   -25.5839689598008988,     22.0699164999425008,     14.5902026036779997)
+TOP2013     : pos=(   -25.5839689598009237,     22.0699164999424546,     14.5902026036780406)
+arcmin_error = 0.000000
+
+tt = -73000.0
+HelioVector : pos=(    36.4035708396756021,    -11.7473067389593009,    -14.6304139635222992)
+TOP2013     : pos=(    36.4035708396755808,    -11.7473067389593133,    -14.6304139635222619)
+arcmin_error = 0.000000
+
+tt = -36500.0
+HelioVector : pos=(    10.2436041239516999,     44.5280986402284995,     10.8048664487065995)
+TOP2013     : pos=(    10.2436041239517408,     44.5280986402285350,     10.8048664487065693)
+arcmin_error = 0.000000
+
+tt = 0.0
+HelioVector : pos=(    -9.8753695807738993,    -27.9789262247366999,     -5.7537118247043004)
+TOP2013     : pos=(    -9.8753695807739028,    -27.9789262247367461,     -5.7537118247042978)
+arcmin_error = 0.000000
+
+tt = 36500.0
+HelioVector : pos=(    39.7009143866163967,     28.4327664903825017,     -3.0906026170880998)
+TOP2013     : pos=(    39.7009143866164109,     28.4327664903825230,     -3.0906026170880683)
+arcmin_error = 0.000000
+
+tt = 73000.0
+HelioVector : pos=(   -27.3620419812794999,     18.4265651225706009,     13.9975343005914006)
+TOP2013     : pos=(   -27.3620419812795213,     18.4265651225705902,     13.9975343005913544)
+arcmin_error = 0.000000
+
+tt = 109500.0
+HelioVector : pos=(    38.3556091850032033,     -8.7643800131841996,    -14.2951819118807002)
+TOP2013     : pos=(    38.3556091850032033,     -8.7643800131842013,    -14.2951819118807286)
+arcmin_error = 0.000000
+
+tt = 146000.0
+HelioVector : pos=(     7.3929490279056003,     44.3826789515344018,     11.6295002148542999)
+TOP2013     : pos=(     7.3929490279056296,     44.3826789515344231,     11.6295002148543460)
+arcmin_error = 0.000000
+
+(Expect simulation differences...)
+tt = -164125.0
+HelioVector : pos=(    35.4241202565997995,    -13.0838961196795953,    -14.7486042553514931)
+TOP2013     : pos=(    35.4240537144631702,    -13.0815716642804176,    -14.7477155126910375)
+arcmin_error = 0.211103
+
+tt = -127625.0
+HelioVector : pos=(    11.5241325901806722,     44.5182113626220612,     10.4111066092492237)
+TOP2013     : pos=(    11.5231470636593958,     44.5175128022429547,     10.4109656507676931)
+arcmin_error = 0.088676
+
+tt = -91125.0
+HelioVector : pos=(   -11.6243877502725041,    -27.2611743158390816,     -5.0012598618039528)
+TOP2013     : pos=(   -11.6242275682680951,    -27.2645343523108856,     -5.0023119538620842)
+arcmin_error = 0.403100
+
+tt = -54625.0
+HelioVector : pos=(    40.2183556155567814,     27.4593327837268362,     -3.5484622626858044)
+TOP2013     : pos=(    40.2172449951831297,     27.4597185252380598,     -3.5480786789161254)
+arcmin_error = 0.087072
+
+tt = -18125.0
+HelioVector : pos=(   -26.7100827824359790,     19.8981115521309739,     14.2551462887341351)
+TOP2013     : pos=(   -26.7092765882112495,     19.8984453657328935,     14.2550034970384640)
+arcmin_error = 0.083899
+
+tt = 18375.0
+HelioVector : pos=(    37.6125960087146254,     -9.9763217209973760,    -14.4459503349835696)
+TOP2013     : pos=(    37.6112540047620385,     -9.9737364446586518,    -14.4446487322043033)
+arcmin_error = 0.264249
+
+tt = 54875.0
+HelioVector : pos=(     8.5166532438548153,     44.4557557474244902,     11.3080656988418848)
+TOP2013     : pos=(     8.5162250004929234,     44.4545071614106178,     11.3076670772183938)
+arcmin_error = 0.101603
+
+tt = 91375.0
+HelioVector : pos=(    -7.4175021427046914,    -28.7885128261721235,     -6.7522720720512748)
+TOP2013     : pos=(    -7.4167093935591151,    -28.7898441810303503,     -6.7527408441642400)
+arcmin_error = 0.182544
+
+tt = 127875.0
+HelioVector : pos=(    38.9369560746144501,     29.7760202114865038,     -2.4364778277987855)
+TOP2013     : pos=(    38.9359077127982118,     29.7755990080719428,     -2.4364552274540632)
+arcmin_error = 0.079157
+
+tt = 164375.0
+HelioVector : pos=(   -28.2074298286596807,     16.3566948758541137,     13.6092594206036281)
+TOP2013     : pos=(   -28.2049513580399562,     16.3578962339715197,     13.6090548255666199)
+arcmin_error = 0.268730
+
+gravsim_test.c: Pluto score = 0.205303 arcmin, over 10 data points.
--- a/generate/gravsim/run
+++ b/generate/gravsim/run
@@ -1,8 +1,10 @@
 #!/bin/bash

-pushd .. > /dev/null
-./generate source || exit 1
-popd > /dev/null
+if [[ "$1" != "skipgen" ]]; then
+    pushd .. > /dev/null
+    ./generate source || exit $?
+    popd > /dev/null
+fi

 gcc -Wall -Werror -o gravsim_test \
    -D ASTRONOMY_ENGINE_GRAVSIM_LOG \
@@ -12,7 +14,8 @@ gcc -Wall -Werror -o gravsim_test \
    -I ../vsop/ \
    ../../source/c/astronomy.c \
    ../top2013/top2013.c \
-    gravsim_test.c -lm || exit 1
+    gravsim_test.c -lm || exit $?

 rm -f gravsim.log
-./gravsim_test
+./gravsim_test || exit $?
+exit 0
--- a/generate/run
+++ b/generate/run
@@ -131,6 +131,7 @@ echo "Generating EQJ/GAL conversion test data."
 echo ""

 ./makedoc || exit $?
+( cd gravsim && ./run skipgen && cd .. ) || exit $?
 ./unit_test_csharp $1 || exit $?
 ./unit_test_js $1 || exit $?
 ./unit_test_c $1 || exit $?