aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/find_peaks.c107
-rw-r--r--src/find_peaks.h4
-rw-r--r--src/fit.c6
-rw-r--r--src/misc.c56
-rw-r--r--src/misc.h3
-rw-r--r--src/test-find-peaks.c58
-rw-r--r--src/test.c86
-rw-r--r--src/viridis.pal287
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)