/* See copyright information at the end of the file */

/*! @file simulation.c
 * @brief Generate simulation data.
 *
 * @code
 * Usage: simulation [r(otation)|t(ranslation)] [OPTIONS] dir
 * Options:
 * 	-b [bond_coeff]
 * 	-c [linear|exponential|geometric|lundy-mees]
 * 	-d [data file]
 * 	-l [log level] where log level is one of [0|NONE] or [1|LOG] or
 * [2|WARN] or [3|VERBOSE]
 *	-m [cg|fire|quickmin|sd]
 *	-p [apot]
 *	-v [v0] (in fJ)
 * 	-A [number of annealing steps]
 *	-B [beta for langevin schedules]
 * 	-C [acoll]
 *      -M [minimization step tolerance] (in fJ)
 *	-R [runtime (in LAMMPS units)]
 *	-S [langevin seed]
 *	-T [number of consecutive minimizations]
 * 	-x [translation vector, x component (units of apot)]
 * 	-y [translation vector, y component (units of apot)]
 *
 * Rotation mode options:
 * 	-f [first angle]
 * 	-s [angular step]
 * 	-t [last angle]
 *
 * Translation mode options:
 * 	-a [initial angle for translational simulation]
 * @endcode
 */

#include <assert.h>
#include <ctype.h>
#include <dirent.h>
#include <errno.h>
#include <fcntl.h>
#include <math.h>
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <unistd.h>

#include <sys/stat.h>

#define STD_IMPLEMENTATION
#include "std.h"
#define COMMON_IMPLEMENTATION
#include "common.h"
#include "inputs.h"

#define DATA "data/test.lmpdat"
#define ACOL 5.8
#define APOT 5.2
#define MIN_STYLE "fire"
#define TIMESTEP 500
#define PKB 1.38e-8
#define PV0 1e-8
#define MAIL "marioforzanini@lcm.mi.infn.it"

/*! Initial value for non dynamic metadata entries, i.e. those that
 * may be set by the user and are not automatically determined by
 * dir_metadata */
static Metadata default_metadata = {
	/*! One of the cooling schedules found in
	 *  https://doi.org/10.4236/am.2017.88090
	 */
	.cooling_fn = cooling_linear,
	/*! LAMMPS input data, generated with lattice_scheme_generator */
	.data = STRING_STATIC("data/test.lmpdat"),
	.min_style = STRING_STATIC("fire"),
	.acol = ACOL,
	.angle = 0.0,
	.apot = APOT,
	.beta = -1.0,
	.bond_coeff = 1e-7,
	.v0 = PV0,
	/*! Default temperature is chosen such that \f$kb * T = 10\% v0\f$ (in
	   K) */
	.temperature = PV0 / PKB / 10,
	.annealing = 0,
	.damping = 10 * TIMESTEP,
	.ntry = 0,
	.seed = 3548735,
	.timestep = TIMESTEP,
	.mode = MODE_UNKNOWN,
	.step = STEP,
	.from = 0,
	.to = 60,
	.translation = { 0, 0 },
	.tolerance = 1e-15,
};

/*! Number of cpus to use per node when run under SLURM */
static int cpus = 4;

/*! Take a list of lines and specialize LAMMPS commands
 * substituting the chosen parameters, writing the result to out.
 *
 * @see inputs.h
 * @todo It may be wiser to make this into a table to ease
 * mantainability, e.g.:
 * @code
 * typedef void (SubstituteFn*)(void);
 * struct Substitution {
 *     String pattern;
 *     SubstituteFn fn;
 * } substitutions[] = {
 *     {STRING_STATIC("read_data"), read_data_fn},
 *     ...
 * };
 * @endcode
 * But function pointers introduce indirection, probably hurting
 * performance, either way if the number of recognised patterns
 * becomes too large the tradeoff may be worthwhile. It would also be
 * wiser to put this into a hash table to avoid branching and linearly
 * searching the table.
 *
 * @param[in] out File to write the results to
 * @param[in] template Array of lines to be used as a template.
 * @param template_len Number of lines in the template.
 * @param curr_step Current step for the simulation, its meaning depends on
 * the simulation mode: it may be an angular step or a translation step
 * @param[in] m Metadata for the simulation.
 */
void
file_from_template(FILE out[static 1], String template[static 1],
	size_t template_len, int curr_step, Metadata *m)
{
	for (int i = 0; i < template_len; i++) {
		String line = template[i];

		if (string_starts_with(
			    STRING_LITERAL("<BEGIN_TRIES>"), line)) {
			if (m->ntry > 1) {
				fprintf(out,
					"variable Ntry_plus_one equal %d\n"
					"variable try loop 1 "
					"${Ntry_plus_one}\n"
					"label try_loop\n"
					"reset_timestep 0\n"
					"log logfire_r%d_t${try}.lammps\n",
					m->ntry + 1, curr_step);
			}
		} else if (string_starts_with(
				   STRING_LITERAL("<END_TRIES>"), line)) {
			if (m->ntry > 1) {
				fprintf(out,
					"write_data "
					"conf_finale_r%d_t${try}.data\n"
					"next try\n"
					"velocity all set 0.0 0.0 0.0\n"
					"jump SELF try_loop\n",
					curr_step);
			} else {
				fprintf(out,
					"write_data "
					"conf_finale_r%d_t0.data\n",
					curr_step);
			}
		} else if (string_starts_with(
				   STRING_LITERAL("write_data"), line)) {
			fprintf(out, "write_data conf_finale_r%d.data\n",
				curr_step);
		} else if (string_starts_with(
				   STRING_LITERAL("bond_coeff"), line)) {
			fprintf(out, "bond_coeff 1 %.6e %.2f\n",
				m->bond_coeff, m->acol);
		} else if (string_starts_with(
				   STRING_LITERAL("read_data"), line)) {
			fprintf(out, "read_data %s\n", basename(m->data.str));
		} else if (string_starts_with(
				   STRING_LITERAL(
					   "displace_atoms all rotate"),
				   line)) {
			fprintf(out,
				"displace_atoms all rotate 0.0 0.0 "
				"0.0 0.0 "
				"0.0 1.0 %.2f\n",
				m->angle);
		} else if (string_starts_with(
				   STRING_LITERAL("displace_atoms all move"),
				   line)) {
			fprintf(out,
				"displace_atoms all move %.3f %.3f 0.0\n",
				m->translation[0], m->translation[1]);
		} else if (string_starts_with(
				   STRING_LITERAL("<COOLING>"), line)) {
			if (m->annealing > 0) {
				m->cooling_fn(out, m->ntry, m->damping,
					m->seed, m->runtime, m->temperature,
					m->annealing, m->beta);
			}
		} else if (string_starts_with(
				   STRING_LITERAL("variable V0"), line)) {
			fprintf(out, "variable V0 equal %.6e\n", m->v0);
		} else if (string_starts_with(
				   STRING_LITERAL("variable Nstep"), line)) {
			fprintf(out, "variable Nstep equal %d\n", m->nstep);
		} else if (string_starts_with(
				   STRING_LITERAL("min_style"), line)) {
			fprintf(out, "min_style %s\n", m->min_style.str);
		} else if (string_starts_with(
				   STRING_LITERAL("variable c"), line)) {
			fprintf(out, "variable c equal 2*PI/%.2f\n", m->apot);
		} else if (string_starts_with(STRING_LITERAL("log "), line)) {
			fprintf(out, "log logfire_r%d_t0.lammps\n",
				curr_step);
		} else if (string_starts_with(
				   STRING_LITERAL("timestep "), line)) {
			fprintf(out, "timestep %zu\n", m->timestep);
		} else if (string_starts_with(
				   STRING_LITERAL("variable V0final equal"),
				   line)
			&& !string_ends_with(STRING_LITERAL("0"), line)) {
			fprintf(out, "variable V0final equal %.6e\n", m->v0);
		} else if (string_starts_with(
				   STRING_LITERAL("<UP OR DOWN LOG>"),
				   line)) {
			fprintf(out,
				"if \"V0final == 0\" then \"log "
				"logfire_r%d_v${i}_down.lammps\" else \"log "
				"logfire_%d_v${i}_up.lammps\"\n",
				curr_step, curr_step);
		} else if (string_starts_with(
				   STRING_LITERAL("variable forceTolerance"),
				   line)) {
			fprintf(out, "variable forceTolerance equal %e\n",
				m->tolerance);
		} else {
			fprintf(out, String_Fmt "\n", String_Arg(line));
		}
	}
}

/*! Create an executable file or open it for writing.
 *
 * @param[in] run_path Filename.
 * @param[out] run File pointer of the open file.
 * @return errno(3).
 */
int
fopen_run_script(const char run_path[static 1], FILE *run[static 1])
{
	int runfd;

	if ((runfd = open(run_path, O_WRONLY | O_CREAT,
		     S_IRWXU | S_IRWXG | S_IROTH))
		< 0) {
		if (errno == EEXIST) {
			if ((runfd = open(run_path, O_WRONLY,
				     S_IRWXG | S_IRWXU | S_IROTH))
				< 0) {
				return errno;
			}
		} else {
			DBG("Error opening %s: %s.\n Continuing anyway.\n",
				run_path, strerror(errno));
			return errno;
		}
	}

	if (!(*run = fdopen(runfd, "w"))) {
		return errno;
	}

	return 0;
}

/*! Angles simulation.
 *
 * Create nsdirname and, using angle_template, generate input files
 * for a simulation involving each angle in the interval [from, to].
 *
 * Also create a script orchestrating the simulation, if run under the
 * SLURM scheduler.
 *
 * @param[in] r Global Region.
 * @param[in] nsdirname Simulation directory (without trailing slash).
 * @param[in] m Metadata for the simulation
 * @return true if no errors occurred, false otherwise.
 */
bool
angles(Region *r, const char nsdirname[static 1], Metadata m[static 1])
{
	char *path, *run_path, *data_path, *basedata;
	size_t dirlen, path_len, run_len, data_len;
	int errnum;
	bool error_occurred = false;
	FILE *run = NULL;
	DIR *d;

	assert(m->to >= m->from);
	basedata = basename(m->data.str);
	dirlen = strlen(nsdirname);

	if (!(d = opendir(nsdirname))) {
		if (errno == ENOENT) {
			if (mkdir(nsdirname, S_IRWXU | S_IRWXG | S_IROTH)
				!= 0) {
				DBG("Cannot create directory %s: %s\n",
					nsdirname, strerror(errno));
				goto err;
			}
		} else {
			DBG("Cannot create directory %s: %s\n", nsdirname,
				strerror(errno));
			goto err;
		}
	}

	data_len = strlen(basedata) + dirlen + 2;
	data_path = region_malloc(r, data_len);
	if (snprintf(data_path, data_len, "%s/%s", nsdirname, basedata) < 0) {
		DBG("Encoding error.%s", "\n");
		goto err;
	}
	if ((errnum = cp(m->data.str, data_path)) != 0) {
		DBG("Error copying " String_Fmt " to %s: %s. Continuing "
		    "anyway.\n",
			String_Arg(m->data), strerror(errnum), data_path);
		error_occurred = true;
	}

	if (!write_metadata(r, nsdirname, m->data.str, m)) {
		DBG("Error writing metadata, continuing anyway.%s", "\n");
		error_occurred = true;
		goto err;
	}

	path_len = dirlen + LATTICE_FILENAME_LEN + 2;
	path = region_malloc(r, path_len);

	run_len = dirlen + RUN_FILENAME_LEN + 2;
	run_path = region_malloc(r, run_len);
	if (snprintf(run_path, run_len, "%s/" RUN_FILENAME, nsdirname) < 0) {
		DBG("Encoding error.%s", "\n");
		goto err;
	}

	if ((errnum = fopen_run_script(run_path, &run)) != 0) {
		DBG("Error opening %s: %s.\n Continuing anyway.\n", run_path,
			strerror(errnum));
		error_occurred = true;
	} else {
		fprintf(run,
			"#!/bin/sh\n"
			"#SBATCH --cpus-per-task=%d\n"
			"#SBATCH --nodes=1\n"
			"#SBATCH --mem-per-cpu=1G\n"
			"#SBATCH --mail-type=END,FAIL\n"
			"#SBATCH --mail-user=" MAIL "\n"
			"#SBATCH --array=%d-%d\n",
			cpus, m->from, m->to);
		fprintf(run, "cp %s \"$SLURM_TMPDIR\"\n", basedata);
		fprintf(run,
			"cp lattice_r${SLURM_ARRAY_TASK_ID}.in "
			"\"$SLURM_TMPDIR/\"\n"
			"(\n"
			"\tcd \"$SLURM_TMPDIR/\"\n"
			"\tmpirun lmp_mpi -i "
			"lattice_r${SLURM_ARRAY_TASK_ID}.in\n)\n"
			"cp $SLURM_TMPDIR/"
			"logfire_r${SLURM_ARRAY_TASK_ID}*.lammps .\n"
			"cp $SLURM_TMPDIR/"
			"conf_finale_r${SLURM_ARRAY_TASK_ID}*.data .\n");
		fclose(run);
	}

	for (int angle = m->from; angle <= m->to; angle++) {
		FILE *fp;
		assert(angle < 1000);
		if (snprintf(path, path_len, "%s/lattice_r%d.in", nsdirname,
			    angle)
			< 0) {
			DBG("Encoding error.%s", "\n");
			goto err;
		}
		if (!(fp = fopen(path, "w"))) {
			DBG("Error opening %s for writing: %s\n", path,
				strerror(errno));
			goto err;
		}
		m->angle = angle * m->step;
		file_from_template(
			fp, angle_template, angle_template_len, angle, m);
		fclose(fp);
	}

	if (d)
		closedir(d);
	return !error_occurred;
err:
	if (d)
		closedir(d);
	return false;
}

/*! Generate input files for translation simulations.
 *
 * Rigidly translate the system along a chosen direction, then relax
 * it. This should make it possible to approximate the curvature of
 * the minimum at (0,0).
 *
 * @param[in] r Global Region.
 * @param[in] nsdirname Simulation directory (without trailing slash).
 * @param[in] m Metadata for the simulation.
 * @return true if no errors occurred, false otherwise.
 */
bool
translate(Region *r, const char nsdirname[static 1],
	Metadata m[restrict static 1])
{
	DIR *d = NULL;
	char *path, *run_path, *data_path, *basedata;
	size_t path_len, dirlen, run_len, data_len;
	int errnum;
	bool error_occurred = false;
	double norm = 0,
	       unit_direction[2] = { m->translation[0], m->translation[1] };
	FILE *run = NULL;

	dirlen = strlen(nsdirname);
	assert(m->to >= m->from);
	basedata = basename(m->data.str);

	if (!(d = opendir(nsdirname))) {
		if (errno == ENOENT) {
			if (mkdir(nsdirname, S_IRWXU | S_IRWXG | S_IROTH)
				!= 0) {
				DBG("Cannot create directory %s: %s\n",
					nsdirname, strerror(errno));
				goto err;
			}
		} else {
			DBG("Cannot create directory %s: %s\n", nsdirname,
				strerror(errno));
			goto err;
		}
	}

	data_len = strlen(basedata) + dirlen + 2;
	data_path = region_malloc(r, data_len);
	if (snprintf(data_path, data_len, "%s/%s", nsdirname, basedata) < 0) {
		DBG("Encoding error.%s", "\n");
		goto err;
	}
	if ((errnum = cp(m->data.str, data_path)) != 0) {
		DBG("Error copying " String_Fmt " to %s: %s. Continuing "
		    "anyway.\n",
			String_Arg(m->data), strerror(errnum), data_path);
		error_occurred = true;
	}

	if (!write_metadata(r, nsdirname, m->data.str, m)) {
		DBG("Error writing %s/" METADATA_FILENAME
		    ": %s. Continuing anyway.\n",
			nsdirname, strerror(errno));
		error_occurred = true;
	}

	path_len = dirlen + LATTICE_FILENAME_LEN + 2;
	path = region_malloc(r, path_len);

	run_len = dirlen + RUN_FILENAME_LEN + 2;
	run_path = region_malloc(r, run_len);
	snprintf(run_path, run_len, "%s/" RUN_FILENAME, nsdirname);

	if ((errnum = fopen_run_script(run_path, &run)) != 0) {
		DBG("Error opening %s: %s. Continuing anyway.\n", run_path,
			strerror(errno));
		error_occurred = true;
	} else {
		fprintf(run,
			"#!/bin/sh\n"
			"#SBATCH --cpus-per-task=%d\n"
			"#SBATCH --nodes=1\n"
			"#SBATCH --mem-per-cpu=1G\n"
			"#SBATCH --mail-type=END,FAIL\n"
			"#SBATCH --mail-user=" MAIL "\n"
			"#SBATCH --array=0-%d\n",
			cpus, m->nstep);
		fprintf(run, "cp %s \"$SLURM_TMPDIR\"\n", basedata);
		fprintf(run,
			"cp lattice_r${SLURM_ARRAY_TASK_ID}.in "
			"\"$SLURM_TMPDIR/\"\n"
			"(\n"
			"\tcd \"$SLURM_TMPDIR/\"\n"
			"\tmpirun lmp_mpi -i "
			"lattice_r${SLURM_ARRAY_TASK_ID}.in\n)\n"
			"cp $SLURM_TMPDIR/"
			"logfire_r${SLURM_ARRAY_TASK_ID}*.lammps .\n"
			"cp $SLURM_TMPDIR/"
			"conf_finale_r${SLURM_ARRAY_TASK_ID}*.data .\n");
		fclose(run);
	}

	assert(m->nstep != 0);
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wfloat-equal"
	assert(m->translation[0] != 0 && m->translation[1] != 0);
#pragma GCC diagnostic pop
	norm = sqrt(m->translation[0] * m->translation[0]
		+ m->translation[1] * m->translation[1]);
	unit_direction[0] /= norm;
	unit_direction[1] /= norm;
	for (int step = 0; step <= m->nstep; step++) {
		FILE *fp;
		if (snprintf(path, path_len, "%s/lattice_r%d.in", nsdirname,
			    step)
			< 0) {
			DBG("Encoding error.%s", "\n");
			goto err;
		}
		if (!(fp = fopen(path, "w"))) {
			DBG("Error opening %s for writing: %s\n", path,
				strerror(errno));
			goto err;
		}
		m->mode |= MODE_ANGLES;
		m->translation[0] = unit_direction[0] * step * norm * m->apot
			/ m->nstep;
		m->translation[1] = unit_direction[1] * step * norm * m->apot
			/ m->nstep;
		file_from_template(
			fp, angle_template, angle_template_len, step, m);
		fclose(fp);
	}

	if (d)
		closedir(d);
	return !error_occurred;
err:
	if (d)
		closedir(d);
	return false;
}

static void
usage(const char program[static 1])
{
	fprintf(stderr,
		"Usage: %s [r(otation)|t(ranslation)] [OPTIONS] dir\n"
		"Options:\n"
		"\t-b [bond_coeff]\n"
		"\t-c [" LINEAR_NAME "|" EXPONENTIAL_NAME "|" GEOMETRIC_NAME
		"|" LUNDY_MEES_NAME "]\n"
		"\t-d [data file]\n"
		"\t-l [log level] where log level is one of [0|NONE] or "
		"[1|LOG] or [2|WARN] or [3|VERBOSE]\n"
		"\t-m [cg|fire|quickmin|sd]\n"
		"\t-p [apot]\n"
		"\t-v [v0] (in fJ)\n"
		"\t-A [number of annealing steps]\n"
		"\t-B [beta for langevin schedules]\n"
		"\t-C [acoll]\n"
		"\t-M [minimization step tolerance] (in fJ)\n"
		"\t-R [runtime (in LAMMPS units)]\n"
		"\t-S [langevin seed]\n"
		"\t-T [number of consecutive minimizations]\n"
		"\t-x [translation vector, x component (units of apot)]\n"
		"\t-y [translation vector, y component (units of apot)]\n"
		"\n"
		"Rotation mode options:\n"
		"\t-f [first angle]\n"
		"\t-s [angular step]\n"
		"\t-t [last angle]\n"
		"\n"
		"Translation mode options:\n"
		"\t-a [initial angle for translational simulation]\n",
		program);
	exit(EXIT_FAILURE);
}

int
main(int argc, char *argv[])
{
	Region *r = region_alloc(4 * 1024 * 1024);
	int opt;
	double fto = 6.0, ffrom = 0;
	char *filename = NULL;
	bool angle_specified = false, nstep_specified = false,
	     tx_specified = false, ty_specified = false;
	Metadata m = default_metadata;

	if (argc < 3) {
		usage(argv[0]);
	}
	if (strcmp(argv[1], "rotation") == 0 || strcmp(argv[1], "r") == 0) {
		m.mode = MODE_ANGLES;
	} else if (strcmp(argv[1], "translation") == 0
		|| strcmp(argv[1], "t") == 0) {
		m.mode = MODE_TRANSLATE;
	} else {
		DBG("Unknown mode: %s. Aborting.\n", argv[1]);
		region_free(r);
		exit(EXIT_FAILURE);
	}
	argv++;
	argc--;

	srand((unsigned int)time(NULL));
	while ((opt = getopt(argc, argv,
			"ha:b:d:c:f:l:m:n:p:s:t:v:x:y:A:B:C:M:R:S:T:"))
		!= -1) {
		String opt_string = make_string(strlen(optarg), optarg);
		switch (opt) {
		case 'a':
			angle_specified = true;
			if (!string_strtod(opt_string, &m.angle)) {
				DBG("Expected initial angle, found: %s. "
				    "Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			break;
		case 'b':
			if (!string_strtod(opt_string, &m.bond_coeff)) {
				DBG("Expected bond coefficient, found: %s. "
				    "Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			break;
		case 'd':
			m.data = opt_string;
			break;
		case 'c':
			if (strcmp(optarg, LINEAR_NAME) == 0) {
				m.cooling_fn = cooling_linear;
			} else if (strcmp(optarg, EXPONENTIAL_NAME) == 0) {
				m.cooling_fn = cooling_exponential;
			} else if (strcmp(optarg, GEOMETRIC_NAME) == 0) {
				m.cooling_fn = cooling_geometric;
			} else if (strcmp(optarg, LUNDY_MEES_NAME) == 0) {
				m.cooling_fn = cooling_lundy_mees;
			} else {
				DBG("Unknown cooling schedule %s, expected "
				    "one of " LINEAR_NAME
				    ", " EXPONENTIAL_NAME ", " GEOMETRIC_NAME
				    ", " LUNDY_MEES_NAME ". Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			break;
		case 'f':
			if (m.mode == MODE_TRANSLATE) {
				DBG("Unrecognized option 'f' in translation "
				    "mode.%s",
					"\n");
				region_free(r);
				exit(EXIT_FAILURE);
			}
			if (!string_strtod(opt_string, &ffrom)) {
				DBG("Expected first angle, "
				    "found: %s. Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			m.from = (int)floor(ffrom / m.step);
			break;
		case 'h':
			usage(argv[0]);
			break;
		case 'l':
			if (strcmp(optarg, "NONE") == 0
				|| strcmp(optarg, "0") == 0) {
				log_level = NO_LOG;
			} else if (strcmp(optarg, "WARN") == 0
				|| strcmp(optarg, "1") == 0) {
				log_level = WARN;
			} else if (strcmp(optarg, "LOG") == 0
				|| strcmp(optarg, "2") == 0) {
				log_level = LOG;
			} else if (strcmp(optarg, "VERBOSE") == 0
				|| strcmp(optarg, "3") == 0) {
				log_level = VERBOSE;
			} else {
				usage(argv[0]);
			}
			break;
		case 'm':
			if (strcmp(optarg, "fire") == 0
				|| strcmp(optarg, "cg") == 0
				|| strcmp(optarg, "quickmin") == 0
				|| strcmp(optarg, "sd") == 0) {
				m.min_style = opt_string;
			} else {
				usage(argv[0]);
			}
			break;
		case 'n':
			if (m.mode == MODE_ANGLES) {
				DBG("Unrecognized option 'n' in rotation "
				    "mode.%s",
					"\n");
				region_free(r);
				exit(EXIT_FAILURE);
			}
			if (!string_strtoi(opt_string, &m.nstep)) {
				DBG("Expected number of steps, defaulting to "
				    "100.%s",
					"\n");
			} else {
				nstep_specified = true;
			}
			break;
		case 'p':
			if (!string_strtod(opt_string, &m.apot)) {
				DBG("Expected apot, found: %s. "
				    "Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			break;
		case 's':
			if (m.mode == MODE_TRANSLATE) {
				DBG("Unrecognized option 's' in translation "
				    "mode.%s",
					"\n");
				region_free(r);
				exit(EXIT_FAILURE);
			}
			if (!string_strtod(opt_string, &m.step)) {
				DBG("Expected angular step, found: %s. "
				    "Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			m.to = (int)floor(fto / m.step);
			m.from = (int)floor(ffrom / m.step);
			break;
		case 't':
			if (m.mode == MODE_TRANSLATE) {
				DBG("Unrecognized option 't' in translation "
				    "mode.%s",
					"\n");
				region_free(r);
				exit(EXIT_FAILURE);
			}
			if (!string_strtod(opt_string, &fto)) {
				DBG("Expected last angle, "
				    "found: %s. Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			m.to = (int)floor(fto / m.step);
			break;
		case 'v':
			if (!string_strtod(opt_string, &m.v0)) {
				DBG("Expected substrate potential depth, "
				    "found: %s. Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			break;
		case 'x':
			tx_specified = true;
			if (!string_strtod(opt_string, &m.translation[0])) {
				DBG("Expected x component of translation "
				    "vector, found: %s. Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			break;
		case 'y':
			ty_specified = true;
			if (!string_strtod(opt_string, &m.translation[1])) {
				DBG("Expected y component of translation "
				    "vector, found: %s. Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			break;
		case 'A':
			if (!string_strtoi(opt_string, &m.annealing)) {
				DBG("Expected number of annealing steps, "
				    "found: %s. Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			break;
		case 'B':
			if (!string_strtod(opt_string, &m.beta)) {
				DBG("Expected cooling schedule parameter, "
				    "found: %s. Aborting\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			break;
		case 'C':
			if (!string_strtod(opt_string, &m.acol)) {
				DBG("Expected acoll, found: %s. Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			break;
		case 'M':
			if (!string_strtod(opt_string, &m.tolerance)) {
				DBG("Expected tolerance, found: %s. "
				    "Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			break;
		case 'R':
			if (!string_strtoi(opt_string, &m.runtime)) {
				DBG("Expected number of runtime timesteps, "
				    "found: %s. Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			break;
		case 'S':
			if (!string_strtoi(opt_string, &m.seed)) {
				DBG("Expected langevin seed, found: %s. "
				    "Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			break;
		case 'T':
			if (!string_strtoi(opt_string, &m.ntry)) {
				DBG("Expected number of tries, found: %s. "
				    "Aborting.\n",
					optarg);
				region_free(r);
				exit(EXIT_FAILURE);
			}
			break;
		default:
			usage(argv[0]);
		}
	}

	VERBOSE("bond_coeff = %.6e\n", m.bond_coeff);
	VERBOSE("data = " String_Fmt "\n", String_Arg(m.data));
	VERBOSE("from = %.1f\n", ffrom);
	VERBOSE("to = %.1f\n", fto);
	VERBOSE("potential v0 = %.6e fJ\n", m.v0);
	if (m.annealing > 0) {
		VERBOSE("annealing steps = %d\n", m.annealing);
		VERBOSE("annealing seed = %d\n", m.seed);
		VERBOSE("annealing runtime = %d\n", m.runtime);
		VERBOSE("annealing schedule = %s\n",
			cooling_schedule_name(m.cooling_fn));
	}

	filename = argv[optind];
	strip_slash(filename);
	if (m.mode == MODE_ANGLES) {
		m.nstep = m.to - m.from + 1;
		if (!angles(r, filename, &m)) {
			goto err;
		}
	} else if (m.mode == MODE_TRANSLATE) {
		if (!nstep_specified) {
			m.nstep = 100;
		}
		if (!tx_specified) {
			DBG("Specify the x component of the translation "
			    "vector with -x.%s",
				"\n");
			region_free(r);
			exit(EXIT_FAILURE);
		} else if (!ty_specified) {
			DBG("Specify the y component of the translation "
			    "vector with -y.%s",
				"\n");
			region_free(r);
			exit(EXIT_FAILURE);
		} else if (!angle_specified) {
			DBG("Specify an initial angle of rotation with -a.%s",
				"\n");
			region_free(r);
			exit(EXIT_FAILURE);
		}
		if (!translate(r, filename, &m)) {
			goto err;
		}
	}
	region_free(r);
	return 0;
err:
	fprintf(stderr,
		"Errors occurred, the simulation input might not be correct. "
		"Good luck.%s",
		"\n");
	region_free(r);
	exit(EXIT_FAILURE);
}

/*
 * Copyright ©️ 2023 Mario Forzanini <mf@marioforzanini.com>
 *
 * This file is part of my bachelor thesis.
 *
 * This file is free software: you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * This file is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this file. If not, see <https://www.gnu.org/licenses/>.
 *
 */