#include "fracture.h"

double *genfunc_uniform(unsigned int L, bound_t boundary, 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, bound_t boundary, 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;
	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);
			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);
				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;
				}
			}
		}
		lattice[2*start + 2 * i] = x;
		lattice[2*start + 2 * i + 1] = y;
	}

	return lattice;
}