From cc6a32bb98b65efb2743cd5d9792b49852f4b5bc Mon Sep 17 00:00:00 2001 From: Don Cross Date: Wed, 10 Nov 2021 21:00:38 -0500 Subject: [PATCH] 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. --- generate/gravsim/gravsim_test.c | 82 +++++++++++++++++-------- generate/gravsim/original.txt | 103 ++++++++++++++++++++++++++++++++ generate/gravsim/run | 13 ++-- generate/run | 1 + 4 files changed, 168 insertions(+), 31 deletions(-) create mode 100644 generate/gravsim/original.txt diff --git a/generate/gravsim/gravsim_test.c b/generate/gravsim/gravsim_test.c index cf90c04a..d5fba8c0 100644 --- a/generate/gravsim/gravsim_test.c +++ b/generate/gravsim/gravsim_test.c @@ -6,19 +6,26 @@ */ #include +#include #include #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; diff --git a/generate/gravsim/original.txt b/generate/gravsim/original.txt new file mode 100644 index 00000000..38d44e57 --- /dev/null +++ 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. diff --git a/generate/gravsim/run b/generate/gravsim/run index 8676ffef..dd7fa28b 100755 --- 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 diff --git a/generate/run b/generate/run index ae87a782..0003e964 100755 --- 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 $?