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.c36
1 files changed, 18 insertions, 18 deletions
diff --git a/src/graph_genfunc.c b/src/graph_genfunc.c
index 991e127..39a8480 100644
--- a/src/graph_genfunc.c
+++ b/src/graph_genfunc.c
@@ -13,6 +13,10 @@ double *genfunc_uniform(unsigned int L, bound_t boundary, gsl_rng *r, unsigned i
return lattice;
}
+double g(double rho, double dist) {
+ return 1 - gsl_sf_exp(-M_PI * rho * dist);
+}
+
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);
@@ -49,25 +53,21 @@ double *genfunc_hyperuniform(unsigned int L, bound_t boundary, gsl_rng *r, unsig
while(reject) {
x = gsl_ran_flat(r, 0, 1);
y = gsl_ran_flat(r, 0, 1);
- reject = false;
+ double pp = 1;
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);
- double min_val = 100;
- for (unsigned int k = 0; k < 5; k++) {
- if (min_val > ds[k]) {
- min_val = ds[k];
- }
- }
- if (1-gsl_sf_exp(-M_PI * rho * min_val) < gsl_ran_flat(r, 0, 1)) {
- reject = true;
- break;
- }
- }
+ 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;