From 2bb0740b68fdb62d45adc00204b3990ca42ade77 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 22 Aug 2016 10:11:14 -0400 Subject: started repo again without all the data files gunking the works --- src/randfuncs.c | 71 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) create mode 100644 src/randfuncs.c (limited to 'src/randfuncs.c') diff --git a/src/randfuncs.c b/src/randfuncs.c new file mode 100644 index 0000000..260c047 --- /dev/null +++ b/src/randfuncs.c @@ -0,0 +1,71 @@ + +#include "fracture.h" + +double *genfunc_uniform(unsigned int L, 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, 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; + for (unsigned int i = 0; i < (*num); 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 * i] = x; + lattice[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; +} -- cgit v1.2.3-70-g09d2