From 5ba4109f0021e7b2c9c66821461742a339e07355 Mon Sep 17 00:00:00 2001 From: pants Date: Wed, 7 Dec 2016 13:29:51 -0500 Subject: added support for hyperuniform lattices using existing hard-sphere jamming routine --- src/graph_genfunc.c | 53 +---------------------------------------------------- 1 file changed, 1 insertion(+), 52 deletions(-) (limited to 'src/graph_genfunc.c') 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; } -- cgit v1.2.3-70-g09d2