From 28749fd671aae37fa123196d0a0e95e2e4e11184 Mon Sep 17 00:00:00 2001
From: Don Cross <cosinekitty@gmail.com>
Date: Mon, 6 Jul 2020 14:47:07 -0400
Subject: [PATCH] Wrote scripts to plot Pluto distance, comparing NOVAS vs
 TOP2013.

The data plot confirms there is some kind of high-frequency error
that causes excessive local minima/maxima that fools the apside
finder.
---
 generate/.gitignore  |  1 +
 generate/ctest.c     | 58 ++++++++++++++++++++++++++++++++++++++++++++
 generate/generate.c  |  2 +-
 generate/plot_pluto  |  7 ++++++
 generate/plotdist.py | 39 +++++++++++++++++++++++++++++
 5 files changed, 106 insertions(+), 1 deletion(-)
 create mode 100755 generate/plot_pluto
 create mode 100755 generate/plotdist.py

diff --git a/generate/.gitignore b/generate/.gitignore
index 934615ff..87cf87ca 100644
--- a/generate/.gitignore
+++ b/generate/.gitignore
@@ -2,6 +2,7 @@ generate_c_docs
 generate
 lnxp1600p2200.405
 hygdata_v3.csv
+pluto_*.csv
 TOP2013.dat
 temp/
 html/
diff --git a/generate/ctest.c b/generate/ctest.c
index d452a611..2ce139d3 100644
--- a/generate/ctest.c
+++ b/generate/ctest.c
@@ -90,6 +90,7 @@ static int LocalSolarEclipseTest(void);
 static int LocalSolarEclipseTest1(void);
 static int LocalSolarEclipseTest2(void);
 static int Transit(void);
+static int DistancePlot(astro_body_t body, double ut1, double ut2, const char *filename);
 
 typedef int (* unit_test_func_t) (void);
 
@@ -179,6 +180,29 @@ int main(int argc, const char *argv[])
                 goto success;
             }
         }
+
+        if (argc == 6)
+        {
+            if (!strcmp(verb, "distplot"))
+            {
+                /* ctest distplot name tt1 tt2 */
+                astro_body_t body;
+                double ut1, ut2;
+
+                body = Astronomy_BodyCode(argv[2]);
+                if (body == BODY_INVALID)
+                    FAIL("Invalid body name '%s'\n", argv[2]);
+
+                if (1 != sscanf(argv[3], "%lf", &ut1))
+                    FAIL("Invalid tt1 value '%s'\n", argv[3]);
+
+                if (1 != sscanf(argv[4], "%lf", &ut2))
+                    FAIL("Invalid tt2 value '%s'\n", argv[4]);
+
+                CHECK(DistancePlot(body, ut1, ut2, argv[5]));
+                goto success;
+            }
+        }
     }
 
     FAIL("ctest: Invalid command line arguments.\n");
@@ -3210,3 +3234,37 @@ fail:
 }
 
 /*-----------------------------------------------------------------------------------------------------------*/
+
+static int DistancePlot(astro_body_t body, double ut1, double ut2, const char *filename)
+{
+    const int npoints = 100000;
+    int error = 1;
+    int i;
+    double ut;
+    astro_time_t time, j2000;
+    astro_func_result_t dist;
+    FILE *outfile = NULL;
+
+    outfile = fopen(filename, "wt");
+    if (outfile == NULL)
+        FAIL("DistancePlot: Cannot open output file: %s\n", filename);
+
+    j2000 = Astronomy_MakeTime(2000, 1, 1, 12, 0, 0.0);
+
+    fprintf(outfile, "\"tt\",\"distance\"\n");
+    for (i=0; i < npoints; ++i)
+    {
+        ut = ut1 + (((double)i)/((double)(npoints-1)) * (ut2 - ut1));
+        time = Astronomy_AddDays(j2000, ut);
+        dist = Astronomy_HelioDistance(body, time);
+        CHECK_STATUS(dist);
+        fprintf(outfile, "%0.16lf,%0.16lg\n", time.tt, dist.value);
+    }
+
+    error = 0;
+fail:
+    if (outfile != NULL) fclose(outfile);
+    return error;
+}
+
+/*-----------------------------------------------------------------------------------------------------------*/
diff --git a/generate/generate.c b/generate/generate.c
index ff2fb701..8a8f7168 100644
--- a/generate/generate.c
+++ b/generate/generate.c
@@ -1809,7 +1809,7 @@ static int DistancePlot(const char *name, double tt1, double tt2)
     double pos[3];
     double tt, dist;
     int i;
-    const int npoints = 1000;
+    const int npoints = 100000;
 
     CHECK(OpenEphem());
 
diff --git a/generate/plot_pluto b/generate/plot_pluto
new file mode 100755
index 00000000..3e8c0159
--- /dev/null
+++ b/generate/plot_pluto
@@ -0,0 +1,7 @@
+#!/bin/bash
+rm -f pluto_{novas,top}.csv
+./build || exit $?
+./generate distplot Pluto -50000 +50000 > pluto_novas.csv || exit $?
+./ctbuild || exit $?
+./ctest distplot Pluto -200000 +200000 pluto_top.csv || exit $?
+./plotdist.py pluto_novas.csv pluto_top.csv || exit $?
diff --git a/generate/plotdist.py b/generate/plotdist.py
new file mode 100755
index 00000000..ae4491bc
--- /dev/null
+++ b/generate/plotdist.py
@@ -0,0 +1,39 @@
+#!/usr/bin/env python3
+
+import sys
+import re
+import matplotlib.pyplot as plt
+
+
+def LoadCsv(inFileName):
+    xlist = []
+    ylist = []
+    with open(inFileName, 'rt') as infile:
+        lnum = 0
+        for line in infile:
+            line = line.strip()
+            lnum += 1
+            if lnum > 1:
+                token = line.split(',')
+                if len(token) < 2:
+                    raise Exception('Invalid CSV line {} in file {}'.format(lnum, inFileName))
+                xlist.append(float(token[0]))
+                ylist.append(float(token[1]))
+    return xlist, ylist
+
+
+def PlotData(novasFileName, topFileName):
+    novas_tlist, novas_rlist = LoadCsv(novasFileName)
+    top_tlist, top_rlist = LoadCsv(topFileName)
+    plt.plot(novas_tlist, novas_rlist, 'b.')
+    plt.plot(top_tlist, top_rlist, 'r.')
+    plt.show()
+    plt.close('all')
+    return 0
+
+
+if __name__ == '__main__':
+    if len(sys.argv) != 3:
+        print('USAGE: plotdist.py novas_dist.csv top_dist.csv')
+        sys.exit(1)
+    sys.exit(PlotData(sys.argv[1], sys.argv[2]))