From 0ad947f800bcbe2c488d2d5cbcdb16c46e6d3857 Mon Sep 17 00:00:00 2001
From: Jaron <jaron@kent-dobias.com>
Date: Tue, 8 Nov 2016 08:17:26 -0500
Subject: various changes, including adding central moments and changing the
 fuse thresholds to long doubles

---
 src/cen_anal_process.c | 154 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 154 insertions(+)
 create mode 100644 src/cen_anal_process.c

(limited to 'src/cen_anal_process.c')

diff --git a/src/cen_anal_process.c b/src/cen_anal_process.c
new file mode 100644
index 0000000..98f68ac
--- /dev/null
+++ b/src/cen_anal_process.c
@@ -0,0 +1,154 @@
+
+#include "fracture.h"
+#include <gsl/gsl_sf_erf.h>
+#include <gsl/gsl_sf_laguerre.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+
+void get_mean(uint_t len, uint32_t *data, double result[2]) {
+	uint_t total = 0;
+	double mean = 0;
+	double mean_err = 0;
+
+	for (uint_t i = 0; i < len; i++) {
+		uint32_t datai = data[i];
+
+		total += datai;
+		mean += i * datai;
+		mean_err += gsl_pow_2(i) * datai; 
+	}
+
+	double meanf = mean / total;
+	double meanf_err = meanf * sqrt(mean_err / gsl_pow_2(mean) + 1 / total);
+
+	result[0] = meanf;
+	result[1] = meanf_err;
+}
+
+void get_mom(uint_t len, uint_t n, uint32_t *data, double mean[2], double result[2]) {
+	uint_t total = 0;
+	double mom = 0;
+	double mom_err = 0;
+
+	for (uint_t i = 0; i < len; i++) {
+		uint32_t datai = data[i];
+		double in = pow(fabs(((double)i) - mean[0]), n);
+
+		total += datai;
+		mom += in * datai;
+		mom_err += gsl_pow_2(in) * datai; 
+	}
+
+	double momf = mom / total;
+	double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total);
+
+	result[0] = momf;
+	result[1] = momf_err;
+}
+
+int main(int argc, char *argv[]) {
+	uint_t nc = argc - 1;
+	uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
+	double *betas_c = (double *)malloc(nc * sizeof(double));
+	double *vals_c1 = (double *)malloc(nc * sizeof(double));
+	double *errors_c1 = (double *)malloc(nc * sizeof(double));
+	double *vals_c2 = (double *)malloc(nc * sizeof(double));
+	double *errors_c2 = (double *)malloc(nc * sizeof(double));
+	double *vals_c3 = (double *)malloc(nc * sizeof(double));
+	double *errors_c3 = (double *)malloc(nc * sizeof(double));
+	double *vals_c4 = (double *)malloc(nc * sizeof(double));
+	double *errors_c4 = (double *)malloc(nc * sizeof(double));
+
+	char *out_filename = (char *)malloc(100 * sizeof(char));
+	uint_t out_filename_len = 0;
+
+	for (uint_t i = 0; i < nc; i++) {
+		char *fn = argv[1 + i];
+		char *buff = (char *)malloc(20 * sizeof(char));
+		uint_t pos = 0; uint_t c = 0;
+		while (c < 5) {
+			if (fn[pos] == '_') {
+				c++;
+			}
+			if (i == 0) {
+				out_filename[pos] = fn[pos];
+				out_filename_len++;
+			}
+			pos++;
+		}
+		c = 0;
+		while (fn[pos] != '_') {
+			buff[c] = fn[pos];
+			pos++;
+			c++;
+		}
+		buff[c] = '\0';
+		Ls_c[i] = atoi(buff);
+		c = 0;
+		pos++;
+		while (fn[pos] != '_') {
+			buff[c] = fn[pos];
+			pos++;
+			c++;
+		}
+		buff[c] = '\0';
+		betas_c[i] = atof(buff);
+		free(buff);
+
+		uint_t dist_len = 4 * pow(Ls_c[i], 2);
+		uint32_t *dist = malloc(dist_len * sizeof(uint32_t));
+		FILE *dist_file = fopen(fn, "rb");
+		fread(dist, sizeof(uint32_t), dist_len, dist_file);
+		fclose(dist_file);
+		{
+			double mom1[2];
+			get_mean(dist_len, dist, mom1);
+			vals_c1[i] = mom1[0];
+			errors_c1[i] = mom1[1];
+
+			double mom2[2];
+			get_mom(dist_len, 2, dist, mom1, mom2);
+			vals_c2[i] = mom2[0];
+			errors_c2[i] = mom2[1];
+
+			double mom3[2];
+			get_mom(dist_len, 3, dist, mom1, mom3);
+			vals_c3[i] = mom3[0];
+			errors_c3[i] = mom3[1];
+
+			double mom4[2];
+			get_mom(dist_len, 4, dist, mom1, mom4);
+			vals_c4[i] = mom4[0];
+			errors_c4[i] = mom4[1];
+		}
+		free(dist);
+	}
+
+	out_filename[out_filename_len-1] = '.';
+	out_filename[out_filename_len] = 't';
+	out_filename[out_filename_len+1] = 'x';
+	out_filename[out_filename_len+2] = 't';
+	out_filename[out_filename_len+3] = '\0';
+
+	FILE *out_file = fopen(out_filename, "w");
+
+	for (uint_t i = 0; i < nc; i++) {
+		fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]);
+	}
+
+	fclose(out_file);
+
+
+	free(Ls_c);
+	free(betas_c);
+	free(vals_c1);
+	free(errors_c1);
+	free(vals_c2);
+	free(errors_c2);
+	free(vals_c3);
+	free(errors_c3);
+
+	return 0;
+}
+
-- 
cgit v1.2.3-70-g09d2