update find_peaks algorithm

Previously, the algorithm used to find peaks was to search for all peaks in the
Hough transform above some constant fraction of the highest peak. This
algorithm could have issues finding smaller peaks away from the highest peak.
The new algorithm instead finds the highest peak in the Hough transform and
then recomputes the Hough transform ignoring all PMT hits within the Cerenkov
cone of the first peak. The next peak is found from this transform and the
process is iteratively repeated until a certain number of peaks are found.

One disadvantage of this new system is that it will *always* find the same
number of peaks and this will usually be greater than the actual number of
rings in the event. This is not a problem though since when fitting the event
we loop over all possible peaks and do a quick fit to determine the starting
point and so false positives are OK because the real peaks will fit better
during this quick fit.

Another potential issue with this new method is that by rejecting all PMT hits
within the Cerenkov cone of the first peak we could miss a second peak very
close to the first peak. This is partially mitigated by the fact that when we
loop over all possible combinations of the particle ids and directions we allow
each peak to be used more than once. For example, when fitting for the
hypothesis that an event is caused by two electrons and one muon and given two
possible directions 1 and 2, we will fit for the following possible direction
combinations:

    1 1 1
    1 1 2
    1 2 1
    1 2 2
    2 2 1
    2 2 2

Therefore if there is a second ring close to the first it is possible to fit it
correctly since we will seed the quick fit with two particles pointing in the
same direction.

This commit also adds a few tests for new functions and changes the energy step
size during the quick fit to 10% of the starting energy value.

8 files changed, 568 insertions, 39 deletions

diff --git a/src/find_peaks.c b/src/find_peaks.c
index 05d8c15..e9030a0 100644
--- a/src/find_peaks.c
+++ b/src/find_peaks.c
@@ -6,6 +6,7 @@
 #include <stdlib.h> /* for malloc() */
 #include "optics.h"
 #include <math.h> /* for exp() */
+#include "misc.h"
 
 typedef struct peak {
     size_t i;
@@ -115,14 +116,37 @@ end:
     return;
 }
 
-void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, size_t m, double *result)
+void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, size_t m, double *result, double *last, size_t len)
 {
     /* Computes the "Hough transform" of the event `ev` and stores it in `result`. */
     size_t i, j, k;
-    double dir[3], pmt_dir[3], cos_theta;
+    double dir[3], pmt_dir[3], cos_theta2;
+    int skip;
+    double *sin_theta, *cos_theta, *sin_phi, *cos_phi;
+    double wavelength0, n_d2o;
 
-    double wavelength0 = 400.0;
-    double n_d2o = get_index_snoman_d2o(wavelength0);
+    sin_theta = malloc(sizeof(double)*n);
+    cos_theta = malloc(sizeof(double)*n);
+    sin_phi = malloc(sizeof(double)*n);
+    cos_phi = malloc(sizeof(double)*n);
+
+    wavelength0 = 400.0;
+    n_d2o = get_index_snoman_d2o(wavelength0);
+
+    /* Precompute sin(theta), cos(theta), sin(phi), and cos(phi) to speed
+     * things up. */
+    for (i = 0; i < n; i++) {
+        sin_theta[i] = sin(x[i]);
+        cos_theta[i] = cos(x[i]);
+    }
+
+    for (i = 0; i < m; i++) {
+        sin_phi[i] = sin(y[i]);
+        cos_phi[i] = cos(y[i]);
+    }
+
+    /* Zero out the result array. */
+    for (i = 0; i < n*m; i++) result[i] = 0.0;
 
     for (i = 0; i < MAX_PMTS; i++) {
         if (pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue;
@@ -131,55 +155,82 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n,
 
         normalize(pmt_dir);
 
+        /* Ignore PMT hits within the Cerenkov cone of previously found peaks. */
+        skip = 0;
+        for (j = 0; j < len; j++) {
+            if (DOT(pmt_dir,last+3*j) > 1/n_d2o) {
+                skip = 1;
+                break;
+            }
+        }
+
+        if (skip) continue;
+
         for (j = 0; j < n; j++) {
             for (k = 0; k < m; k++) {
-                double r = 1+x[j]*x[j]+y[k]*y[k];
-                dir[0] = 2*x[j]/r;
-                dir[1] = 2*y[k]/r;
-                dir[2] = (-1+x[j]*x[j]+y[k]*y[k])/r;
+                dir[0] = sin_theta[j]*cos_phi[k];
+                dir[1] = sin_theta[j]*sin_phi[k];
+                dir[2] = cos_theta[j];
 
-                cos_theta = DOT(pmt_dir,dir);
+                cos_theta2 = DOT(pmt_dir,dir);
 
-                result[j*m + k] += ev->pmt_hits[i].qhs*exp(-fabs(cos_theta-1.0/n_d2o)/0.1);
+                result[j*n + k] += ev->pmt_hits[i].qhs*exp(-fabs(cos_theta2-1.0/n_d2o)/0.1);
             }
         }
     }
+
+    free(sin_theta);
+    free(cos_theta);
+    free(sin_phi);
+    free(cos_phi);
 }
 
-void find_peaks(event *ev, double *pos, size_t n, size_t m, double *peak_theta, double *peak_phi, size_t *npeaks, size_t max_peaks, double threshold)
+void find_peaks(event *ev, double *pos, size_t n, size_t m, double *peak_theta, double *peak_phi, size_t *npeaks, size_t max_peaks)
 {
+    /* Finds rings in the event `ev` by looking for peaks in the Hough
+     * transform. The directions of potential particle rings are stored in the
+     * arrays `peak_theta` and `peak_phi`. The number of rings found are stored
+     * in the variable `npeaks`.
+     *
+     * This function first computes the Hough transform of the event `ev` and
+     * looks for the highest peak. The next peak is then found by computing the
+     * Hough transform ignoring any PMTs within the Cerenkov cone of the first
+     * peak. This process is repeated until `max_peaks` peaks are found. */
     size_t i;
+    double *x, *y, *result, *dir;
+    size_t *imax, *jmax;
+    size_t max;
 
-    double *x = calloc(n,sizeof(double));
-    double *y = calloc(m,sizeof(double));
-    double *result = calloc(n*m,sizeof(double));
+    x = calloc(n,sizeof(double));
+    y = calloc(m,sizeof(double));
+    result = malloc(n*m*sizeof(double));
+    dir = calloc(max_peaks*3,sizeof(double));
 
     for (i = 0; i < n; i++) {
-        x[i] = -10 + 20.0*i/(n-1);
+        x[i] = i*M_PI/(n-1);
     }
 
     for (i = 0; i < m; i++) {
-        y[i] = -10 + 20.0*i/(m-1);
+        y[i] = i*2*M_PI/(m-1);
     }
 
-    get_hough_transform(ev,pos,x,y,n,m,result);
+    imax = calloc(max_peaks,sizeof(size_t));
+    jmax = calloc(max_peaks,sizeof(size_t));
 
-    size_t *imax = calloc(max_peaks,sizeof(size_t));
-    size_t *jmax = calloc(max_peaks,sizeof(size_t));
-
-    find_peaks_array(result,n,m,imax,jmax,npeaks,max_peaks,threshold);
-
-    for (i = 0; i < *npeaks; i++) {
-        double R, theta;
-        R = sqrt(x[imax[i]]*x[imax[i]]+y[jmax[i]]*y[jmax[i]]);
-        theta = atan2(y[jmax[i]],x[imax[i]]);
-        peak_theta[i] = 2*atan(1/R);
-        peak_phi[i] = theta;
+    for (i = 0; i < max_peaks; i++) {
+        get_hough_transform(ev,pos,x,y,n,m,result,dir,i);
+        max = argmax(result,n*m);
+        peak_theta[i] = x[max/m];
+        peak_phi[i] = y[max % m];
+        get_dir(dir+i*3,peak_theta[i],peak_phi[i]);
     }
 
+    *npeaks = max_peaks;
+
     free(imax);
     free(jmax);
     free(x);
     free(y);
     free(result);
+    free(dir);
 }
diff --git a/src/find_peaks.h b/src/find_peaks.h
index 2abf6fe..7ca798b 100644
--- a/src/find_peaks.h
+++ b/src/find_peaks.h
@@ -5,7 +5,7 @@
 #include <stddef.h> /* for size_t */
 
 void find_peaks_array(double *x, size_t n, size_t m, size_t *imax, size_t *jmax, size_t *npeaks, size_t max_peaks, double threshold);
-void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, size_t m, double *result);
-void find_peaks(event *ev, double *pos, size_t n, size_t m, double *peak_theta, double *peak_phi, size_t *npeaks, size_t max_peaks, double threshold);
+void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, size_t m, double *result, double *last, size_t len);
+void find_peaks(event *ev, double *pos, size_t n, size_t m, double *peak_theta, double *peak_phi, size_t *npeaks, size_t max_peaks);
 
 #endif
diff --git a/src/fit.c b/src/fit.c
index cb5a11d..79d8c8e 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -5296,7 +5296,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double
     ub[3] = GTVALID;
 
     /* Find the peaks in the Hough transform of the event. */
-    find_peaks(ev, pos, 100, 100, peak_theta, peak_phi, &npeaks, LEN(peak_theta), 0.5);
+    find_peaks(ev, pos, 100, 100, peak_theta, peak_phi, &npeaks, LEN(peak_theta));
 
     result = malloc(sizeof(size_t)*n*ipow(npeaks,n));
 
@@ -5321,7 +5321,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double
     fpars.hit_only = 0;
     fpars.print = 0;
     fpars.charge_only = 0;
-    nlopt_set_ftol_abs(opt, 1e-2);
+    nlopt_set_ftol_abs(opt, 1e-5);
     nlopt_set_maxeval(opt, 10000);
 
     fpars.n = n;
@@ -5365,7 +5365,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double
             if (T0 < Tmin2) T0 = Tmin2;
 
             x[4+3*j] = T0;
-            ss[4+3*j] = x[4+3*j]*0.02;
+            ss[4+3*j] = x[4+3*j]*0.1;
             ss[5+3*j] = 0.01;
             ss[6+3*j] = 0.01;
             lb[4+3*j] = Tmin;
diff --git a/src/misc.c b/src/misc.c
index a610944..ce27274 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -742,3 +742,59 @@ void combinations_with_replacement(size_t n, size_t r, size_t *result, size_t *l
 
     free(tmp);
 }
+
+size_t argmax(double *a, size_t n)
+{
+    /* Returns the index of the maximum of the array `a`. */
+    size_t i, index;
+    double max;
+
+    if (n < 2) return 0;
+
+    index = 0;
+    max = a[0];
+    for (i = 1; i < n; i++) {
+        if (a[i] > max) {
+            max = a[i];
+            index = i;
+        }
+    }
+
+    return index;
+}
+
+size_t argmin(double *a, size_t n)
+{
+    /* Returns the index of the minimum of the array `a`. */
+    size_t i, index;
+    double min;
+
+    if (n < 2) return 0;
+
+    index = 0;
+    min = a[0];
+    for (i = 1; i < n; i++) {
+        if (a[i] < min) {
+            min = a[i];
+            index = i;
+        }
+    }
+
+    return index;
+}
+
+void get_dir(double *dir, double theta, double phi)
+{
+    /* Compute the 3-dimensional unit vector `dir` from the polar angle `theta`
+     * and azimuthal angle `phi`. */
+    double sin_theta, cos_theta, sin_phi, cos_phi;
+
+    cos_theta = cos(theta);
+    sin_theta = sin(theta);
+    cos_phi = cos(phi);
+    sin_phi = sin(phi);
+
+    dir[0] = sin_theta*cos_phi;
+    dir[1] = sin_theta*sin_phi;
+    dir[2] = cos_theta;
+}
diff --git a/src/misc.h b/src/misc.h
index 02df398..20197bb 100644
--- a/src/misc.h
+++ b/src/misc.h
@@ -31,5 +31,8 @@ size_t ipow(size_t base, size_t exp);
 void product(size_t n, size_t r, size_t *result);
 void unique_vertices(int *id, size_t n, size_t npeaks, size_t *result, size_t *nvertices);
 void combinations_with_replacement(size_t n, size_t r, size_t *result, size_t *len);
+size_t argmax(double *a, size_t n);
+size_t argmin(double *a, size_t n);
+void get_dir(double *dir, double theta, double phi);
 
 #endif
diff --git a/src/test-find-peaks.c b/src/test-find-peaks.c
index c7c2c1d..c82f9ec 100644
--- a/src/test-find-peaks.c
+++ b/src/test-find-peaks.c
@@ -9,12 +9,13 @@
 #include <math.h> /* for M_PI */
 #include <errno.h> /* for errno */
 #include <string.h> /* for strerror() */
+#include "vector.h"
 
 #define EV_RECORD      0x45562020
 #define MCTK_RECORD    0x4d43544b
 #define MCVX_RECORD    0x4d435658
 
-void plot_3d(double *x, double *y, double *z, size_t n, size_t m)
+void plot_hough_transform(double *x, double *y, double *z, size_t n, size_t m)
 {
     size_t i, j;
 
@@ -49,6 +50,51 @@ void plot_3d(double *x, double *y, double *z, size_t n, size_t m)
     }
 }
 
+void plot_find_peaks(event *ev, double *peak_theta, double *peak_phi, size_t n)
+{
+    size_t i;
+
+    FILE *pipe = popen("gnuplot -p", "w");
+
+    if (!pipe) {
+        fprintf(stderr, "error running gnuplot command: %s\n", strerror(errno));
+        exit(1);
+    }
+
+    fprintf(pipe, "set macros\n");
+    fprintf(pipe, "load 'viridis.pal'\n");
+    fprintf(pipe, "set title 'Hough Transform'\n");
+    /* Not entirely sure what these do, but following the instructions from
+     * http://lowrank.net/gnuplot/plot3d2-e.html. */
+    fprintf(pipe, "set xrange [0:6.28]\n");
+    fprintf(pipe, "set yrange [3.14:0]\n");
+    fprintf(pipe, "plot '-' u 1:2:3:4 with circles palette fillstyle solid, '-' u 1:2 with circles lc rgb \"red\" lw 2\n");
+
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (ev->pmt_hits[i].hit) {
+            double r = NORM(pmts[i].pos);
+            double theta = acos(pmts[i].pos[2]/r);
+            double phi = atan2(pmts[i].pos[1],pmts[i].pos[0]);
+            phi = (phi < 0) ? phi + 2*M_PI: phi;
+            fprintf(pipe, "%.10g %.10g %.10g %.10g\n", phi, theta, 0.01, ev->pmt_hits[i].qhs);
+        }
+    }
+    fprintf(pipe,"e\n");
+
+    for (i = 0; i < n; i++) {
+        fprintf(pipe, "%.10g %.10g\n", peak_phi[i], peak_theta[i]);
+    }
+    fprintf(pipe,"e\n");
+
+    /* Pause so that you can rotate the 3D graph. */
+    fprintf(pipe,"pause mouse keypress\n");
+
+    if (pclose(pipe)) {
+        fprintf(stderr, "error closing gnuplot command: %s\n", strerror(errno));
+        exit(1);
+    }
+}
+
 int get_event(zebraFile *f, event *ev, bank *b)
 {
     /* Read all the PMT banks from the zebra file and update `ev`.
@@ -193,16 +239,16 @@ int main(int argc, char **argv)
             double *result = calloc(n*m,sizeof(double));
 
             for (i = 0; i < n; i++) {
-                x[i] = -10 + 20.0*i/(n-1);
+                x[i] = i*M_PI/(n-1);
             }
 
             for (i = 0; i < m; i++) {
-                y[i] = -10 + 20.0*i/(m-1);
+                y[i] = i*2*M_PI/(m-1);
             }
 
-            get_hough_transform(&ev,pos,x,y,n,m,result);
+            get_hough_transform(&ev,pos,x,y,n,m,result,0,0);
 
-            find_peaks(&ev,pos,n,m,peak_theta,peak_phi,&npeaks,LEN(peak_theta),0.5);
+            find_peaks(&ev,pos,n,m,peak_theta,peak_phi,&npeaks,LEN(peak_theta));
 
             printf("gtid %i\n", ev.gtid);
             for (i = 0; i < npeaks; i++) {
@@ -210,7 +256,7 @@ int main(int argc, char **argv)
             }
 
             if (plot)
-                plot_3d(x,y,result,n,m);
+                plot_find_peaks(&ev,peak_theta,peak_phi,npeaks);
 
             /* Skip reading in the next bank since get_event() already read in
              * the next bank. */
diff --git a/src/test.c b/src/test.c
index e1c9e74..7053a94 100644
--- a/src/test.c
+++ b/src/test.c
@@ -25,6 +25,7 @@
 #include "sno.h"
 #include "quad.h"
 #include "find_peaks.h"
+#include <gsl/gsl_sort.h>
 
 typedef int testFunction(char *err);
 
@@ -1893,6 +1894,88 @@ int test_combinations_with_replacement(char *err)
     return 0;
 }
 
+int test_get_dir(char *err)
+{
+    /* Tests the get_dir() function. */
+    size_t i, j;
+    double dir[3], expected_dir[3];
+    double theta, phi;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        theta = genrand_real2()*M_PI;
+        phi = genrand_real2()*2*M_PI;
+
+        get_dir(dir,theta,phi);
+
+        expected_dir[0] = sin(theta)*cos(phi);
+        expected_dir[1] = sin(theta)*sin(phi);
+        expected_dir[2] = cos(theta);
+
+        for (j = 0; j < 3; j++) {
+            if (!isclose(dir[j],expected_dir[j],0,1e-9)) {
+                sprintf(err, "dir[%zu] returned %.2f but expected %.2f", j, dir[j], expected_dir[j]);
+                return 1;
+            }
+        }
+    }
+
+    return 0;
+}
+
+int test_argmax(char *err)
+{
+    /* Tests the argmax() function. */
+    size_t i, j;
+    double a[100];
+    size_t max, p;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        for (j = 0; j < 100; j++)
+            a[j] = genrand_real2();
+
+        max = argmax(a,LEN(a));
+
+        gsl_sort_largest_index(&p,1,a,1,LEN(a));
+
+        if (max != p) {
+            sprintf(err, "argmax() returned %zu but expected %zu", max, p);
+            return 1;
+        }
+    }
+
+    return 0;
+}
+
+int test_argmin(char *err)
+{
+    /* Tests the argmin() function. */
+    size_t i, j;
+    double a[100];
+    size_t min, p;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        for (j = 0; j < 100; j++)
+            a[j] = genrand_real2();
+
+        min = argmin(a,LEN(a));
+
+        gsl_sort_smallest_index(&p,1,a,1,LEN(a));
+
+        if (min != p) {
+            sprintf(err, "argmin() returned %zu but expected %zu", min, p);
+            return 1;
+        }
+    }
+
+    return 0;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -1940,6 +2023,9 @@ struct tests {
     {test_find_peaks_highest, "test_find_peaks_highest"},
     {test_find_peaks_sorted, "test_find_peaks_sorted"},
     {test_combinations_with_replacement, "test_combinations_with_replacement"},
+    {test_get_dir, "test_get_dir"},
+    {test_argmax, "test_argmax"},
+    {test_argmin, "test_argmin"},
 };
 
 int main(int argc, char **argv)
diff --git a/src/viridis.pal b/src/viridis.pal
new file mode 100644
index 0000000..ca3efc2
--- /dev/null
+++ b/src/viridis.pal
@@ -0,0 +1,287 @@
+# New matplotlib colormaps by Nathaniel J. Smith, Stefan van der Walt,
+# and (in the case of viridis) Eric Firing.
+#
+# This file and the colormaps in it are released under the CC0 license /
+# public domain dedication. We would appreciate credit if you use or
+# redistribute these colormaps, but do not impose any legal restrictions.
+#
+# To the extent possible under law, the persons who associated CC0 with
+# mpl-colormaps have waived all copyright and related or neighboring rights
+# to mpl-colormaps.
+#
+# You should have received a copy of the CC0 legalcode along with this
+# work.  If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.
+
+#https://github.com/BIDS/colormap/blob/master/colormaps.py
+
+
+# line styles
+set style line  1 lt 1 lc rgb '#440154' # dark purple
+set style line  2 lt 1 lc rgb '#472c7a' # purple
+set style line  3 lt 1 lc rgb '#3b518b' # blue
+set style line  4 lt 1 lc rgb '#2c718e' # blue
+set style line  5 lt 1 lc rgb '#21908d' # blue-green
+set style line  6 lt 1 lc rgb '#27ad81' # green
+set style line  7 lt 1 lc rgb '#5cc863' # green
+set style line  8 lt 1 lc rgb '#aadc32' # lime green
+set style line  9 lt 1 lc rgb '#fde725' # yellow
+
+
+# palette
+set palette defined (\
+0   0.267004 0.004874 0.329415,\
+1   0.268510 0.009605 0.335427,\
+2   0.269944 0.014625 0.341379,\
+3   0.271305 0.019942 0.347269,\
+4   0.272594 0.025563 0.353093,\
+5   0.273809 0.031497 0.358853,\
+6   0.274952 0.037752 0.364543,\
+7   0.276022 0.044167 0.370164,\
+8   0.277018 0.050344 0.375715,\
+9   0.277941 0.056324 0.381191,\
+10  0.278791 0.062145 0.386592,\
+11  0.279566 0.067836 0.391917,\
+12  0.280267 0.073417 0.397163,\
+13  0.280894 0.078907 0.402329,\
+14  0.281446 0.084320 0.407414,\
+15  0.281924 0.089666 0.412415,\
+16  0.282327 0.094955 0.417331,\
+17  0.282656 0.100196 0.422160,\
+18  0.282910 0.105393 0.426902,\
+19  0.283091 0.110553 0.431554,\
+20  0.283197 0.115680 0.436115,\
+21  0.283229 0.120777 0.440584,\
+22  0.283187 0.125848 0.444960,\
+23  0.283072 0.130895 0.449241,\
+24  0.282884 0.135920 0.453427,\
+25  0.282623 0.140926 0.457517,\
+26  0.282290 0.145912 0.461510,\
+27  0.281887 0.150881 0.465405,\
+28  0.281412 0.155834 0.469201,\
+29  0.280868 0.160771 0.472899,\
+30  0.280255 0.165693 0.476498,\
+31  0.279574 0.170599 0.479997,\
+32  0.278826 0.175490 0.483397,\
+33  0.278012 0.180367 0.486697,\
+34  0.277134 0.185228 0.489898,\
+35  0.276194 0.190074 0.493001,\
+36  0.275191 0.194905 0.496005,\
+37  0.274128 0.199721 0.498911,\
+38  0.273006 0.204520 0.501721,\
+39  0.271828 0.209303 0.504434,\
+40  0.270595 0.214069 0.507052,\
+41  0.269308 0.218818 0.509577,\
+42  0.267968 0.223549 0.512008,\
+43  0.266580 0.228262 0.514349,\
+44  0.265145 0.232956 0.516599,\
+45  0.263663 0.237631 0.518762,\
+46  0.262138 0.242286 0.520837,\
+47  0.260571 0.246922 0.522828,\
+48  0.258965 0.251537 0.524736,\
+49  0.257322 0.256130 0.526563,\
+50  0.255645 0.260703 0.528312,\
+51  0.253935 0.265254 0.529983,\
+52  0.252194 0.269783 0.531579,\
+53  0.250425 0.274290 0.533103,\
+54  0.248629 0.278775 0.534556,\
+55  0.246811 0.283237 0.535941,\
+56  0.244972 0.287675 0.537260,\
+57  0.243113 0.292092 0.538516,\
+58  0.241237 0.296485 0.539709,\
+59  0.239346 0.300855 0.540844,\
+60  0.237441 0.305202 0.541921,\
+61  0.235526 0.309527 0.542944,\
+62  0.233603 0.313828 0.543914,\
+63  0.231674 0.318106 0.544834,\
+64  0.229739 0.322361 0.545706,\
+65  0.227802 0.326594 0.546532,\
+66  0.225863 0.330805 0.547314,\
+67  0.223925 0.334994 0.548053,\
+68  0.221989 0.339161 0.548752,\
+69  0.220057 0.343307 0.549413,\
+70  0.218130 0.347432 0.550038,\
+71  0.216210 0.351535 0.550627,\
+72  0.214298 0.355619 0.551184,\
+73  0.212395 0.359683 0.551710,\
+74  0.210503 0.363727 0.552206,\
+75  0.208623 0.367752 0.552675,\
+76  0.206756 0.371758 0.553117,\
+77  0.204903 0.375746 0.553533,\
+78  0.203063 0.379716 0.553925,\
+79  0.201239 0.383670 0.554294,\
+80  0.199430 0.387607 0.554642,\
+81  0.197636 0.391528 0.554969,\
+82  0.195860 0.395433 0.555276,\
+83  0.194100 0.399323 0.555565,\
+84  0.192357 0.403199 0.555836,\
+85  0.190631 0.407061 0.556089,\
+86  0.188923 0.410910 0.556326,\
+87  0.187231 0.414746 0.556547,\
+88  0.185556 0.418570 0.556753,\
+89  0.183898 0.422383 0.556944,\
+90  0.182256 0.426184 0.557120,\
+91  0.180629 0.429975 0.557282,\
+92  0.179019 0.433756 0.557430,\
+93  0.177423 0.437527 0.557565,\
+94  0.175841 0.441290 0.557685,\
+95  0.174274 0.445044 0.557792,\
+96  0.172719 0.448791 0.557885,\
+97  0.171176 0.452530 0.557965,\
+98  0.169646 0.456262 0.558030,\
+99  0.168126 0.459988 0.558082,\
+100 0.166617 0.463708 0.558119,\
+101 0.165117 0.467423 0.558141,\
+102 0.163625 0.471133 0.558148,\
+103 0.162142 0.474838 0.558140,\
+104 0.160665 0.478540 0.558115,\
+105 0.159194 0.482237 0.558073,\
+106 0.157729 0.485932 0.558013,\
+107 0.156270 0.489624 0.557936,\
+108 0.154815 0.493313 0.557840,\
+109 0.153364 0.497000 0.557724,\
+110 0.151918 0.500685 0.557587,\
+111 0.150476 0.504369 0.557430,\
+112 0.149039 0.508051 0.557250,\
+113 0.147607 0.511733 0.557049,\
+114 0.146180 0.515413 0.556823,\
+115 0.144759 0.519093 0.556572,\
+116 0.143343 0.522773 0.556295,\
+117 0.141935 0.526453 0.555991,\
+118 0.140536 0.530132 0.555659,\
+119 0.139147 0.533812 0.555298,\
+120 0.137770 0.537492 0.554906,\
+121 0.136408 0.541173 0.554483,\
+122 0.135066 0.544853 0.554029,\
+123 0.133743 0.548535 0.553541,\
+124 0.132444 0.552216 0.553018,\
+125 0.131172 0.555899 0.552459,\
+126 0.129933 0.559582 0.551864,\
+127 0.128729 0.563265 0.551229,\
+128 0.127568 0.566949 0.550556,\
+129 0.126453 0.570633 0.549841,\
+130 0.125394 0.574318 0.549086,\
+131 0.124395 0.578002 0.548287,\
+132 0.123463 0.581687 0.547445,\
+133 0.122606 0.585371 0.546557,\
+134 0.121831 0.589055 0.545623,\
+135 0.121148 0.592739 0.544641,\
+136 0.120565 0.596422 0.543611,\
+137 0.120092 0.600104 0.542530,\
+138 0.119738 0.603785 0.541400,\
+139 0.119512 0.607464 0.540218,\
+140 0.119423 0.611141 0.538982,\
+141 0.119483 0.614817 0.537692,\
+142 0.119699 0.618490 0.536347,\
+143 0.120081 0.622161 0.534946,\
+144 0.120638 0.625828 0.533488,\
+145 0.121380 0.629492 0.531973,\
+146 0.122312 0.633153 0.530398,\
+147 0.123444 0.636809 0.528763,\
+148 0.124780 0.640461 0.527068,\
+149 0.126326 0.644107 0.525311,\
+150 0.128087 0.647749 0.523491,\
+151 0.130067 0.651384 0.521608,\
+152 0.132268 0.655014 0.519661,\
+153 0.134692 0.658636 0.517649,\
+154 0.137339 0.662252 0.515571,\
+155 0.140210 0.665859 0.513427,\
+156 0.143303 0.669459 0.511215,\
+157 0.146616 0.673050 0.508936,\
+158 0.150148 0.676631 0.506589,\
+159 0.153894 0.680203 0.504172,\
+160 0.157851 0.683765 0.501686,\
+161 0.162016 0.687316 0.499129,\
+162 0.166383 0.690856 0.496502,\
+163 0.170948 0.694384 0.493803,\
+164 0.175707 0.697900 0.491033,\
+165 0.180653 0.701402 0.488189,\
+166 0.185783 0.704891 0.485273,\
+167 0.191090 0.708366 0.482284,\
+168 0.196571 0.711827 0.479221,\
+169 0.202219 0.715272 0.476084,\
+170 0.208030 0.718701 0.472873,\
+171 0.214000 0.722114 0.469588,\
+172 0.220124 0.725509 0.466226,\
+173 0.226397 0.728888 0.462789,\
+174 0.232815 0.732247 0.459277,\
+175 0.239374 0.735588 0.455688,\
+176 0.246070 0.738910 0.452024,\
+177 0.252899 0.742211 0.448284,\
+178 0.259857 0.745492 0.444467,\
+179 0.266941 0.748751 0.440573,\
+180 0.274149 0.751988 0.436601,\
+181 0.281477 0.755203 0.432552,\
+182 0.288921 0.758394 0.428426,\
+183 0.296479 0.761561 0.424223,\
+184 0.304148 0.764704 0.419943,\
+185 0.311925 0.767822 0.415586,\
+186 0.319809 0.770914 0.411152,\
+187 0.327796 0.773980 0.406640,\
+188 0.335885 0.777018 0.402049,\
+189 0.344074 0.780029 0.397381,\
+190 0.352360 0.783011 0.392636,\
+191 0.360741 0.785964 0.387814,\
+192 0.369214 0.788888 0.382914,\
+193 0.377779 0.791781 0.377939,\
+194 0.386433 0.794644 0.372886,\
+195 0.395174 0.797475 0.367757,\
+196 0.404001 0.800275 0.362552,\
+197 0.412913 0.803041 0.357269,\
+198 0.421908 0.805774 0.351910,\
+199 0.430983 0.808473 0.346476,\
+200 0.440137 0.811138 0.340967,\
+201 0.449368 0.813768 0.335384,\
+202 0.458674 0.816363 0.329727,\
+203 0.468053 0.818921 0.323998,\
+204 0.477504 0.821444 0.318195,\
+205 0.487026 0.823929 0.312321,\
+206 0.496615 0.826376 0.306377,\
+207 0.506271 0.828786 0.300362,\
+208 0.515992 0.831158 0.294279,\
+209 0.525776 0.833491 0.288127,\
+210 0.535621 0.835785 0.281908,\
+211 0.545524 0.838039 0.275626,\
+212 0.555484 0.840254 0.269281,\
+213 0.565498 0.842430 0.262877,\
+214 0.575563 0.844566 0.256415,\
+215 0.585678 0.846661 0.249897,\
+216 0.595839 0.848717 0.243329,\
+217 0.606045 0.850733 0.236712,\
+218 0.616293 0.852709 0.230052,\
+219 0.626579 0.854645 0.223353,\
+220 0.636902 0.856542 0.216620,\
+221 0.647257 0.858400 0.209861,\
+222 0.657642 0.860219 0.203082,\
+223 0.668054 0.861999 0.196293,\
+224 0.678489 0.863742 0.189503,\
+225 0.688944 0.865448 0.182725,\
+226 0.699415 0.867117 0.175971,\
+227 0.709898 0.868751 0.169257,\
+228 0.720391 0.870350 0.162603,\
+229 0.730889 0.871916 0.156029,\
+230 0.741388 0.873449 0.149561,\
+231 0.751884 0.874951 0.143228,\
+232 0.762373 0.876424 0.137064,\
+233 0.772852 0.877868 0.131109,\
+234 0.783315 0.879285 0.125405,\
+235 0.793760 0.880678 0.120005,\
+236 0.804182 0.882046 0.114965,\
+237 0.814576 0.883393 0.110347,\
+238 0.824940 0.884720 0.106217,\
+239 0.835270 0.886029 0.102646,\
+240 0.845561 0.887322 0.099702,\
+241 0.855810 0.888601 0.097452,\
+242 0.866013 0.889868 0.095953,\
+243 0.876168 0.891125 0.095250,\
+244 0.886271 0.892374 0.095374,\
+245 0.896320 0.893616 0.096335,\
+246 0.906311 0.894855 0.098125,\
+247 0.916242 0.896091 0.100717,\
+248 0.926106 0.897330 0.104071,\
+249 0.935904 0.898570 0.108131,\
+250 0.945636 0.899815 0.112838,\
+251 0.955300 0.901065 0.118128,\
+252 0.964894 0.902323 0.123941,\
+253 0.974417 0.903590 0.130215,\
+254 0.983868 0.904867 0.136897,\
+255 0.993248 0.906157 0.143936)