started repo again without all the data files gunking the works

1 files changed, 71 insertions, 0 deletions

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;
+}