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
|
#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_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 < 4) {
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 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\n", Ls_c[i], betas_c[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i]);
}
fclose(out_file);
free(Ls_c);
free(betas_c);
free(vals_c2);
free(errors_c2);
free(vals_c3);
free(errors_c3);
return 0;
}
|
