summaryrefslogtreecommitdiff
path: root/src/graph_genfunc.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/graph_genfunc.c')
-rw-r--r--src/graph_genfunc.c53
1 files changed, 1 insertions, 52 deletions
diff --git a/src/graph_genfunc.c b/src/graph_genfunc.c
index 54d9daa..8dcd0df 100644
--- a/src/graph_genfunc.c
+++ b/src/graph_genfunc.c
@@ -20,58 +20,7 @@ double g(double rho, double dist) {
double *genfunc_hyperuniform(unsigned int L, bound_t boundary, gsl_rng *r, unsigned int *num) {
*num = 2 * pow(L / 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);
- double pp = 1;
- for (unsigned int j = 0; j < i; j++) {
- double ds0, ds1, ds2, ds3, ds4, ds5, ds6, ds7, ds8;
- ds0 = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1],2);
- ds1 = pow(x-lattice[2*j] + 1,2)+pow(y-lattice[2*j+1],2);
- ds2 = pow(x-lattice[2*j] - 1,2)+pow(y-lattice[2*j+1],2);
- ds3 = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1] + 1,2);
- ds4 = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1] - 1,2);
- ds5 = pow(x-lattice[2*j] + 1,2)+pow(y-lattice[2*j+1] + 1,2);
- ds6 = pow(x-lattice[2*j] + 1,2)+pow(y-lattice[2*j+1] - 1,2);
- ds7 = pow(x-lattice[2*j] - 1,2)+pow(y-lattice[2*j+1] + 1,2);
- ds8 = pow(x-lattice[2*j] - 1,2)+pow(y-lattice[2*j+1] - 1,2);
- pp *= g(rho, ds0) * g(rho, ds1) * g(rho, ds2) * g(rho, ds3) * g(rho, ds4) * g(rho, ds5) * g(rho, ds6) * g(rho, ds7) * g(rho, ds8);
- }
- if (pp > gsl_ran_flat(r, 0, 1)) reject = false;
- }
- lattice[2*start + 2 * i] = x;
- lattice[2*start + 2 * i + 1] = y;
- }
+ double *lattice = spheres(*num, 0.8, 0.5, 0.9, 100, 100000);
return lattice;
}