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
|
#include "fracture.h"
double *genfunc_uniform(unsigned int L, bound_t boundary, gsl_rng *r, unsigned int *num) {
*num = pow(L / 2 + 1, 2) + pow((L + 1) / 2, 2);
double *lattice = (double *)malloc(2 * (*num) * sizeof(double));
for (unsigned int i = 0; i < (*num); i++) {
lattice[2*i] = gsl_ran_flat(r, 0, 1);
lattice[2*i+1] = gsl_ran_flat(r, 0, 1);
}
return lattice;
}
double *genfunc_hyperuniform(unsigned int L, bound_t boundary, gsl_rng *r, unsigned int *num) {
*num = pow(L / 2 + 1, 2) + pow((L + 1) / 2, 2);
// necessary to prevent crashing when underflow occurs
gsl_set_error_handler_off();
double *lattice = (double *)malloc(2 * (*num) * sizeof(double));
double rho = *num;
unsigned int to_gen = *num;
unsigned int start = 0;
if (boundary == EMBEDDED_BOUND) {
for (unsigned int i = 0; i < L / 2; i++) {
lattice[2 * i + 1] = 0;
lattice[2 * i] = (2. * i + 1.) / L;
lattice[L + 2 * i + 1] = 1;
lattice[L + 2 * i] = (2. * i + 1.) / L;
lattice[2 * L + 2 * i + 1] = (2. * i + 1.) / L;
lattice[2 * L + 2 * i] = 0;
lattice[3 * L + 2 * i + 1] = (2. * i + 1.) / L;
lattice[3 * L + 2 * i] = 1;
}
to_gen -= 2 * L;
start = 2 * L;
}
for (unsigned int i = 0; i < to_gen; i++) {
bool reject = true;
double x, y;
while(reject) {
x = gsl_ran_flat(r, 0, 1);
y = gsl_ran_flat(r, 0, 1);
reject = false;
for (unsigned int j = 0; j < i; j++) {
double *ds = (double *)malloc(5 * sizeof(double));
ds[0] = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1],2);
ds[1] = pow(x-lattice[2*j] + 1,2)+pow(y-lattice[2*j+1],2);
ds[2] = pow(x-lattice[2*j] - 1,2)+pow(y-lattice[2*j+1],2);
ds[3] = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1] + 1,2);
ds[4] = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1] - 1,2);
unsigned int min_pos = 0; double min_val = 100;
for (unsigned int k = 0; k < 5; k++) {
if (min_val > ds[k]) {
min_pos = k; min_val = ds[k];
}
}
if (1-gsl_sf_exp(-M_PI * rho * min_val) < gsl_ran_flat(r, 0, 1)) {
reject = true;
break;
}
}
}
lattice[2*start + 2 * i] = x;
lattice[2*start + 2 * i + 1] = y;
}
return lattice;
}
void randfunc_flat(gsl_rng *r, double *x, double *y) {
*x = gsl_ran_flat(r, 0, 1);
*y = gsl_ran_flat(r, 0, 1);
}
void randfunc_gaus(gsl_rng *r, double *x, double *y) {
*x = 100;
*y = 100;
double sigma = 0.25;
while (fabs(*x) > 0.5 || fabs(*y) > 0.5) {
gsl_ran_bivariate_gaussian(r, sigma, sigma, 0, x, y);
}
*x += 0.5;
*y += 0.5;
}
|
