diff options
Diffstat (limited to 'src/graph_genfunc.c')
| -rw-r--r-- | src/graph_genfunc.c | 36 |
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; |