summaryrefslogtreecommitdiff
path: root/src/anal_process.c
blob: c66aea085cf7b48364e0c15ff86ce40a6c168a25 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127

#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 mom(uint_t len, uint_t n, uint32_t *data, 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(i, 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));

	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];
			mom(dist_len, 1, dist, mom1);
			vals_c1[i] = mom1[0];
			errors_c1[i] = mom1[1];

			double mom2[2];
			mom(dist_len, 2, dist, mom2);
			vals_c2[i] = mom2[0];
			errors_c2[i] = mom2[1];

			double mom3[2];
			mom(dist_len, 3, dist, mom3);
			vals_c3[i] = mom3[0];
			errors_c3[i] = mom3[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\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]);
	}

	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;
}