/*
 * Monte Carlo simulation of distribution function and drift speed
 * of ions in charge-exchange collision. One dimensional model.
 * Thermal motion of neutral gas taken into account.
 *
 * Developed by Lubomir Host 'rajo' <rajo AT platon.sk>
 * Licensed under terms of GNU General Public License.
 *
 * Compilation & executing:
 *
 * gcc -DDRIFT   -o hotGas-drift    hotGas.c -lm -lcfg+
 * gcc -DDISTRIB -o hotGas-distrib  hotGas.c -lm -lcfg+
 *
 * ./hotGas-drift N=1e6 time=8.0 d_time=0.01 T=0 lambda=3e-7 u=2.0 a=1e10
 *
 * libcfg+ library needed for compilation can be downloaded from
 *
 * http://freshmeat.net/projects/libcfgplus/
 * http://platon.sk/projects/libcfg+/
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <math.h>
#include <time.h>

#include <platon/cfg+.h>

/* type used for real numbers */
#define Real_t           double
/* constant used only for generation Mxw. distribution*/
#define MAX_V_T_CONST    4.5

/* DEBUG rutines, skip to line #140 {{{ */
#ifdef DEBUG

int debug          = 1;
int debug_prefixed = 0;

#define DEBUG_PREFIX "# @@ "

#define _DEBUG_MSG(_level, _msg...)                                \
	/* if _level = 0, print without prefix */                      \
	if (!(_level)) {                                               \
		if (debug_prefixed)                                        \
			fprintf(stderr, "\n");                                 \
		fprintf(stderr, _msg);                                     \
		debug_prefixed = 0;                                        \
	}                                                              \
	else if (debug >= (_level)) {                                  \
		/* first create whole string */                            \
		char debug_string[4096];                                   \
		snprintf(debug_string, 4096, _msg);                        \
		/* if string ends with '\n', prefix it */                  \
		if ((debug_string + strlen(debug_string)) !=               \
				(strrchr(debug_string, '\n') + 1)) {               \
			/* if something was already prefixed, terminate */     \
			if (debug_prefixed) {                                  \
				fprintf(stderr, "%s", debug_string);               \
			}                                                      \
			else {                                                 \
				fprintf(stderr, DEBUG_PREFIX "%.*s%s",             \
						(_level),   "               ",             \
						debug_string);                             \
				debug_prefixed = 1;                                \
			}                                                      \
		}                                                          \
		/* if string doesn't ends with '\n', print as is,      */  \
		/* but prefix it with DEBUG_PREFIX if not prefixed yet */  \
		else {                                                     \
			/* if previous output was prefixed, terminate it   */  \
			if (debug_prefixed) {                                  \
				fprintf(stderr, "\n");                             \
			}                                                      \
			fprintf(stderr, DEBUG_PREFIX "%.*s%s",                 \
					(_level),   "               ",                 \
					debug_string);                                 \
			/* and next output must be prefixed */                 \
			debug_prefixed = 0;                                    \
		}                                                          \
	}                                                              \

/* executes code if debug enabled */
#define _DEBUG(_level, _code)                                      \
	if (debug >= _level) {                                         \
		_code;                                                     \
	}

/* executes code if debug enabled */
#define _DEBUG_BEGIN(_level)                                       \
	if (debug >= _level) {
#define _DEBUG_END(_level)                                         \
	}

int dbg_rnd_idx = 0;
/* "random" speeds returned in DEBUG mode */
Real_t dbg_rnd_speed[] = {
	0.5, 1.5, 1.2, 0.3,
	1.7, 2.5, 1.0, 0.1,
	0.3, -1, -2, -3,
	-4, -5
};
Real_t dbg_rnd_time[]  = {
	2.2, 0.4, 0.6, 0.8, /* t = 4.0 */
	0.8, 0.1, 1.1, 0.5, /* t = 6.5 */
	1.5, 1, 1, 1
};

#else /* DEBUG not defined */

#define _DEBUG_MSG(_level, _msg...)    /* empty */
#define _DEBUG(_level, _code)          /* empty */
#define _DEBUG_BEGIN(_level)           \
	if (0) { /* executes never */
#define _DEBUG_END(_level)             \
	}

#endif /* }}} DEBUG rutines */

#ifndef COMPILE_STRING
#  define COMPILE_STRING "manual"
#endif

char *version = __FILE__ ": " __DATE__ " " __TIME__ /* version {{{ */
#ifdef DEBUG
", DEBUG defined"
#endif
#ifdef DISTRIB
", DISTRIB defined"
#endif
; /* i}}} end of version[] */

/* Variables declaration, default values {{{ */
long int N          = 100000;          /* number of ions*/
Real_t e_charge     = 1.602176462E-19; /* elementary charge [C] */
Real_t k_B          = 1.3806503E-23;   /* Boltzmann constant [J.K^-1] */
Real_t m_u          = 1.66053873E-27;  /* atomic mass constant [kg] */
Real_t T_g          = 0.0;             /* temperature of neutral gas in Kelvins */
Real_t m_g          = 4;               /* molekul mass [atomic mass units] He=4 He_2=8 etc.*/
Real_t lambda       = 0.001;           /* mean free path in meters */
Real_t acceleration = 1.0;             /* acceleration */
Real_t E            = 0.0;             /* intensity of electric field */
Real_t alpha        = 1.0;             /* ratio between energy from field and thermal energy */
Real_t beta         = 0.0;             /* random thermal speeds > beta * v_T */
Real_t start_speed  = 0.0;             /* startup speed of ions */
Real_t abs_acc      = 0.0;             /* absolute value of acceleration */
Real_t max_time     = 10.0;            /* max. time of simulation */
Real_t real_time    = 0.0;             /* real time computed from start of simulation */
Real_t d_time       = 1E-2;            /* time step */
Real_t v_T          = 0.0;             /* thermal speed v_T = sqrt( 2 * k_B * T_g / m_g) */
Real_t v_drift      = 0.0;             /* theoretical drift speed in time = \infty
								        *  v_drift = sqrt(2 * a * lambda / 3.1415) */

Real_t Dimens_x     = 0.0;             /* normalisation constant of position */
Real_t Dimens_v     = 0.0;             /* normalisation constant of speed */
Real_t Dimens_t     = 0.0;             /* normalisation constant of time */

#ifdef DISTRIB
long int **Fx_plus,  **Fv_plus;        /* distribution funkcions over x and over v */
long int **Fx_minus, **Fv_minus;       /* distribution funkcions for negative x, v */
long int Nt    = 100;                  /* # of points in time */
long int Nx    = 100;                  /* number of points in x/v axis */
long int idx   = 0;                    /* index into arrays with distribution function */
Real_t f_max_x = 5.0;                  /* x-limit for distribution function, dimensionless */
Real_t f_max_v = 5.0;                  /* v-limit for distribution function, dimensionless */
Real_t d_f_max_x;                      /* interval in distribution function */
Real_t d_f_max_v;                      /* interval in distribution function */
int output_step = 1;                   /* time step in output of distribution funkction */
#endif
/* }}} Variables declaration */

#define xmalloc(_type, _x, _count)                    /* {{{ */    \
	if ((_count) == 0)                                             \
		fprintf(stderr,                                            \
			"WARNING: memory allocation "                          \
			"(" #_type ") " #_x " of 0 counts ");                  \
		_x = (_type *) malloc((_count) * sizeof(_type));           \
		if (_x == NULL) {                                          \
			fprintf(stderr, "Can't alloc array " #_x               \
				"[%ld] of type (" #_type "), errno = %d - %s\n",   \
				(_count), errno, strerror(errno));                 \
				exit(errno);                                       \
		}                                                          \
	memset(_x, 0, (_count) * sizeof(_type)); /* }}} */

/* return signum of number */
Real_t signum(double x)
{ /* {{{ */
	if (x < 0.0)
		return -1.0;
	else
		return  1.0;
} /* }}} */

/* return random thermal speed */
Real_t random_thermal_speed(void)
{ /* {{{ */
	Real_t ret_v;       /* returned speed */
	Real_t ret_v2;      /* speed square */
	Real_t func;        /* distribution function */

#ifdef DEBUG
	_DEBUG_BEGIN(5)
		/* if debug is on, return speed from dbg_rnd_speed[] array */
		return dbg_rnd_speed[dbg_rnd_idx++];
	_DEBUG_END(5)
#endif

	if (v_T == 0.0)
		return 0.0;

	/*
	 * Maxwell distribution in 1-D model:
	 *
	 *                1                     v_x^2
	 * f_T(v_x) = --------------- exp ( - --------)
	 *             sqrt(\pi) v_T            v_T^2
	 *
	 * v_T = sqrt( 2 k_B T / m_g)
	 */

	/* von Neumann method */
	while (1) {
		ret_v  = -MAX_V_T_CONST + 2.0 * MAX_V_T_CONST * (rand() / (RAND_MAX + 1.0));
		ret_v2 = ret_v * ret_v;
#ifdef DISTRIB_FUNC_3D
		func = 2.7 * ret_v2 * exp(- ret_v2); /* 2.7 = exp(-1) */
#else
		func = exp( - ret_v2);
#endif
		if (fabs(ret_v) > beta
				&& ((1.0 * (rand() / (RAND_MAX + 1.0))) <= func)) {
			return v_T * ret_v;
		}
	}
} /* }}} */

/* return time to the next collison */
Real_t random_time(Real_t v)                               /* [LINE i] */
{ /* {{{ */
	Real_t rand_x, rand_time;
	Real_t z;

#ifdef DEBUG
	_DEBUG_BEGIN(5)
		/* if debug is on, return speed from dbg_rnd_speed[] array */
		_DEBUG_MSG(1, "random_time(%g) = rand_time = %g\n",
				v, dbg_rnd_time[dbg_rnd_idx]);
		return dbg_rnd_time[dbg_rnd_idx];
	_DEBUG_END(5)
#endif
	/* Distribution function:
	 *
	 * f(x) = (1 / lambda) * exp( - x / lambda)
	 *
	 * integral f(x):
	 *
	 * F(x) = 1 - exp( - x / lambda)
	 *
	 * inverse function:
	 *
	 * F^-1(y) = - lambda * ln(1 - y)
	 *
	 */

	/* inversion function method */
	z      = rand() / (RAND_MAX + 1.0);
	/* rand_x    = - lambda * log(1.0 - z); */
	rand_x = -lambda * log(z); /* optimalized */
	/*
	 * Analytical solution of quadr. equation:
	 * 0.5*a*t^2 + v_0*t - lambda = 0
	 * Special case sgn(a) != sgn(v) && v < sqrt(2*lambda*fabs(a))
	 */
	if (v == 0.0) {
		return sqrt(2 * rand_x / abs_acc);
	}
	if (signum(v) != signum(acceleration) && fabs(v) < sqrt(2 * rand_x * abs_acc))
		rand_time = (fabs(v) + sqrt(2 * abs_acc * rand_x - v * v)) / abs_acc;
	else
		rand_time = (2.0 * rand_x) /
			(fabs(v) + sqrt(v * v + 2 * acceleration * rand_x * signum(v)));
	_DEBUG_MSG(1, "random_time(%g) = rand_time = %g\n", v, rand_time);

	return rand_time;

} /* }}} */

/* print parameters after start of simulation */
void print_parameters(void)
{    /* parameter dump {{{ */

	fprintf(stdout,
			"## COMPILATION: %s\n"
			"## VERSION: %s\n"
			"## PARAM:  N                = %ld\n"
			"## PARAM:  time_step        = %f\n"
			"## PARAM:  max_time         = %f\n"
			"## PARAM:  acceleration     = %g\n"
			"## PARAM:  E                = %g\n"
			"## PARAM:  start_speed      = %g\n"
			"## PARAM:  alpha            = %g\n"
			"## PARAM:  beta             = %g\n"
#ifdef DISTRIB
			"## PARAM:  Nt               = %ld\n"
			"## PARAM:  Nx               = %ld\n"
			"## PARAM:  f_max_x          = %f\n"
			"## PARAM:  f_max_v          = %f\n"
			"## PARAM:  output_step      = %d\n"
#endif
			"## PARAM:  T_g              = %f\n"
			"## PARAM:  m_g              = %.10g\n"
			"## PARAM:  v_T              = %g\n"
			"## PARAM:  lambda           = %.5g\n"
			"## VALUE:  real_max_time    = %.10g\n"
			"##         sqrt(lambda/a)   = %.10g\n"
			"## VALUE:  real_d_time      = %.10g\n"
			"## VALUE:  k_B * T_g        = %.10g\n"
			"## VALUE:  lambda * q * E   = %.10g\n"    /*   a*m = q*E   */
			"## (%s field)\n",
			COMPILE_STRING, version,
			/* d_time and max_time must be recomputed, because its values are
			 * adjusted to Dimens_t */
			N, d_time / Dimens_t, max_time / Dimens_t,
			acceleration, E, start_speed / Dimens_v,
			alpha, beta,
#ifdef DISTRIB
			Nt, Nx, f_max_x, f_max_v, output_step,
#endif
			T_g, m_g, v_T, lambda,
			max_time, Dimens_t,
			d_time, k_B * T_g,
			(lambda * abs_acc * m_g),
			(lambda * abs_acc * m_g) < (2 * k_B * T_g) ? "low" : "STRONG");
	fflush(stdout);
} /* }}} */

int main(int argc, char **argv)
{
	double        i = 0;
	long int      j = 0L, k = 0L;    /* counters */
	Real_t        rand_time;
	Real_t        x, v, t, time_left;
#ifdef DISTRIB
	long int      l = 0L;            /* counters */
#else
	Real_t        *X;                /* x coordinate of ion */
	Real_t        *U;                /* drift speed */
	long int      Nb;
#endif
	int            help = 0;

	/* {{{ cmdline parsing */
	/* libcfg+ parsing context */
	CFG_CONTEXT con;

	/* Parsing return code */
	register int ret;

	/* Option set */
	struct cfg_option options[] = {
		{ "help",        'h', NULL, CFG_BOOL,   (void *) &help,         0 },
		{ "N",           'N', NULL, CFG_LONG,   (void *) &N,            0 },
#ifdef DISTRIB
		{ "Nx",          'X', NULL, CFG_LONG,   (void *) &Nx,           0 },
		{ "f_max_x",     'x', NULL, CFG_DOUBLE, (void *) &f_max_x,      0 },
		{ "f_max_v",     'v', NULL, CFG_DOUBLE, (void *) &f_max_v,      0 },
		{ "output_step", 's', NULL, CFG_INT,    (void *) &output_step,  0 },
#endif
		{ "d_time",      'd', NULL, CFG_DOUBLE, (void *) &d_time,       0 },
		{ "time",        't', NULL, CFG_DOUBLE, (void *) &max_time,     0 },
		{ "a",           'a', NULL, CFG_DOUBLE, (void *) &acceleration, 0 },
		{ "alpha",       'A', NULL, CFG_DOUBLE, (void *) &alpha,        0 },
		{ "beta",        'B', NULL, CFG_DOUBLE, (void *) &beta,         0 },
		{ "u",           'u', NULL, CFG_DOUBLE, (void *) &start_speed,  0 },
		{ "T",           'T', NULL, CFG_DOUBLE, (void *) &T_g,          0 },
		{ "mass",        'm', NULL, CFG_DOUBLE, (void *) &m_g,          0 },
		{ "lambda",      'l', NULL, CFG_DOUBLE, (void *) &lambda,       0 },
#ifdef DEBUG
		{ "debug",       'd', NULL, CFG_INT,    (void *) &debug,        0 },
#endif
		CFG_END_OF_LIST
	};

	/* Creating context */
	con = cfg_get_context(options);
	if (con == NULL) {
		puts("Not enough memory");
		return 1;
	}
	cfg_clear_property(con, CFG_LINE_LONG_OPTION_PREFIX);
	cfg_add_properties_type(con, CFG_LINE_LONG_OPTION_PREFIX, "--", "", NULL);

	/* Parse command line from second argument to end (NULL) */
	cfg_set_cmdline_context(con, 1, -1, argv);

	/* Parsing command line */
	ret = cfg_parse(con);
	if (ret != CFG_OK) {
		fprintf(stderr, "error parsing command line: ");
		cfg_fprint_error(con, stderr);
		fprintf(stderr, "\n");
		return ret < 0 ? -ret : ret;
	}

	/* Help screen {{{ */
	if (help) {
		puts("Usage:");
		printf("  %s [options]\n", argv[0]);
		puts("Options:");
		puts("  -h  help                         print this help screen");
		puts("  -N  N=<number>                   number of simulated ions");
		puts("  -d  d_time=<time_step>        [dimensionless]");
		puts("                                   simulation time step");
		puts("  -t  time=<max_time>           [dimensionless]");
		puts("                                   time of simulation");
		puts("  -a  a=<acceleration>          [meter.sec^-2]");
		puts("                                   acceleration of ions");
		puts("  -A  alpha=<koeficient>        [dimensionless]");
		puts("                                   alpha = acceleration * lambda / (v_T * v_T)");
		puts("  -B  beta=<koeficient>         [dimensionless]");
		puts("                                   random thermal speeds > beta * v_T");
		puts("  -u  u=<speed>                 [dimensionless]");
		puts("                                   startup speed of ions normalized to");
		puts("                                   sqrt(lambda*acceleration)");
		puts("  -T  T=<T_g>                   [Kelvin]");
		puts("                                   temperature of neutral gas");
		puts("  -m  mass=<m_g>                [atomic mass units]");
		puts("                                   mass of neutral gas molecule and ions (He=4)");
		puts("  -l  lambda=<mean_free_path>   [meters]");
		puts("                                   mean free path");
#ifdef DISTRIB
		puts("  -X  Nx=<pixels>               [count]");
		puts("                                   pixel per x axis in f(x,t) and/or f(v,t)");
		puts("  -x  f_max_x=<max_x>              range of x axis");
		puts("  -v  f_max_v=<max_v>              range of v axis");
		puts("  -s  output_step=<step>           don't print distrib. f. in every time step");
#endif
#ifdef DEBUG
		puts("  -d  debug=<level>                enable debug level <level>");
#endif
		puts("");
		printf("%s N=1e6 time=10 d_time=1e-2 T=10 lambda=3e-7 u=2.0 a=1e7\n", argv[0]);
		puts("");
	} /* }}} */
	/* }}} cmdline parsing */

	/* initialize random number generator */
#ifdef DEBUG
	srand(1);
#else
	srand(time(NULL));
#endif
	/* compute parameters from input parameters */
	m_g     *= m_u;                           /* from atomic mass units to kilograms */
	lambda   = fabs(lambda);                  /* make mean free path absolute */
	v_T      = sqrt( 2 * k_B * T_g / m_g);    /* thermal speed */
	if (T_g != 0.0)                           /* acceleration set from cmdline is overwrited */
		acceleration = alpha * v_T * v_T / lambda;
	abs_acc  = fabs(acceleration);            /* absolute value of acceleration */
	E        = m_g * acceleration / e_charge; /* intensity of electric field */

	/* normalisation constants */
	Dimens_x = lambda;
	Dimens_t = sqrt(lambda / abs_acc);
	if (acceleration == 0.0)
		Dimens_v = 1.0; /* TRIVIAL case */
	else
		Dimens_v = sqrt(abs_acc * lambda);

	/* theoretical drift speed int time = oo */
	v_drift = signum(acceleration) * sqrt(2 / M_PI) * Dimens_v;

	/* recompute parameters from dimensionless to real */
	start_speed *= Dimens_v;
	max_time    *= Dimens_t;
	d_time      *= Dimens_t;

#ifdef DISTRIB
	Nt = ceil(max_time / d_time) + 1;
	d_f_max_x = f_max_x * Dimens_v * max_time / Nx;
	d_f_max_v = f_max_v * Dimens_v / Nx;

	xmalloc(long int *, Fx_plus, Nt);
	xmalloc(long int *, Fv_plus, Nt);
	xmalloc(long int *, Fx_minus, Nt);
	xmalloc(long int *, Fv_minus, Nt);
	for (j = 0; j < Nt; j++) {
		xmalloc(long int, Fx_plus[j], Nx);
		xmalloc(long int, Fv_plus[j], Nx);
		xmalloc(long int, Fx_minus[j], Nx);
		xmalloc(long int, Fv_minus[j], Nx);
	}
#else
	Nb = ceil(max_time / d_time);
	xmalloc(Real_t, X, Nb);
	xmalloc(Real_t, U, Nb);
#endif
	print_parameters();
	if (help)
		return 1;

	/* loop over all ions, one by one ... */               /* [LINE a] */
	for (i = 0; i < N; i++) {
		/* If start_speed is not set, ions must start with thermal speed */
		if (start_speed == 0.0)
			v = random_thermal_speed();
		else
			v = start_speed;
		t = x = 0.0;
		/* index must be reseted, because we are
		 * going from start with every ion
		 */
		j = 0;

		real_time = d_time;
		rand_time = random_time(v);
		time_left = d_time;

		/* move with ion while real time is < max. time */
		while (real_time < max_time) {
			if ((t + rand_time) > max_time)
				rand_time = max_time - t;
			/* FIRST part of algoritmus */                 /* [LINE c] */
			/* run while rand_time > time step d_time */
			while ((t + rand_time) > real_time) {
				_DEBUG_MSG(2, "X.1.X");
				x    += ((1.0 / 2.0) * acceleration * time_left + v) * time_left;
				v    += acceleration * time_left;
				j++;
#ifdef DISTRIB /* Distribution function or drift speed {{{ */
				/* compute distribution function if DISTRIB defined ... */
				if (x < 0) {
					idx = floorl(-x / d_f_max_x);
					if (idx < Nx) {
						Fx_minus[j][idx] += 1;
					}
				}
				else {
					idx = floorl(x / d_f_max_x);
					/* fprintf(stderr, "# x_idx = %ld\n", idx); */
					if (idx < Nx) {
						/* fprintf(stderr, "x j = %ld\n", j); */
						Fx_plus[j][idx] += 1;
					}
				}
				if (v < 0) {
					idx = floorl(-v / d_f_max_v);
					if (idx < Nx) {
						Fv_minus[j][idx] += 1;
					}
				}
				else {
					idx = floorl(v / d_f_max_v);
					/* fprintf(stderr, "# v_idx = %ld\n", idx); */
					if (idx < Nx) {
						/* fprintf(stderr, "v j = %ld\n", j); */
						Fv_plus[j][idx] += 1;
					}
				}
#else
				/* ... otherwise compute drift sped */
				X[j] += x;
				U[j] += v;
#endif /* }}} */
				_DEBUG_MSG(1, "%g    %g    %g\n", real_time, x, v);
				_DEBUG_MSG(2, "YYYY t = %g real_time = %g rand_time = %g time_left = %g\n",
						t, real_time, rand_time, time_left);
				real_time += d_time;    /* increase real time */
				t         += time_left;
				rand_time  = rand_time - time_left;        /* [LINE h] */
				/* time to next point "time_left" from SECOND part "X.2.X"
				 * is taken into account, run next loop with step d_time
				 * time_left is recomputed after FIRST part of algoritmus
				 */
				time_left = d_time;
			}                                              /* [LINE d] */

			time_left = t + rand_time - real_time + d_time;

			/* If intervals between 2 collision are < time step, we need many
			 * of these intervals for one time step
			 */
			/* SECOND part of algoritmus */                /* [LINE e] */
			do {
				_DEBUG_MSG(2, "X.2.X");
				x   += ((1.0 / 2.0) * acceleration * time_left + v) * time_left;
				t   += time_left;
				/* v   += acceleration * time_left; */
				v    = random_thermal_speed();

				rand_time = random_time(v);
				time_left = rand_time;

				_DEBUG_MSG(2, "ZZZZ t = %g real_time = %g rand_time = %g time_left = %g\n",
						t, real_time, rand_time, time_left);
			} while ((t + rand_time) < real_time           /* [LINE f] */
					&& real_time < max_time);
			_DEBUG_MSG(2, "t = %g real_time = %g rand_time = %g time_left = %g\n",
					t, real_time, rand_time, time_left);

			time_left = real_time - t;                     /* [LINE g] */
		} /* end of loop over time */
	} /* end of loop over ions */                          /* [LINE b] */


	/* Data output {{{ */
	/* Print data to stdout with full description,
	 * parameters are already dumped
	 */
#ifdef DISTRIB
	for (k = 1; k < Nt; k += output_step) {
		fprintf(stdout,
				"\n## Data format    : 1:real_time 2:relat_time  "
				"3:x   4:F(x>0,t)   5:F(x<0,t)   "
				"6:v   7:F(v>0,t)   8:F(v<0,t)\n");
		for (l = 0; l < Nx; l++)
			fprintf(stdout,
					"%.5g""\t%.5g"  "\t%.5g""\t%.8g""\t%.8g"  "\t%.5g""\t%.8g""\t%.8g\n",
					k * d_time, k * d_time / Dimens_t,
					l * d_f_max_x, (1.0 * Fx_plus[k][l]) / N, (1.0 * Fx_minus[k][l]) / N,
					l * d_f_max_v, (1.0 * Fv_plus[k][l]) / N, (1.0 * Fv_minus[k][l]) / N);
		fprintf(stdout, "\n"); /* line separator for gnuplot */
	}
#else
	fprintf(stdout,
			"\n## Data format    : 1:real_time   2:relat_time   "
			"3:(x / time / N)   4:(v_drift / N)   5:x   6:u_drift\n");
	fprintf(stdout,
			"%.5g" "\t%.5g" "\t%.8g" "\t%.8g" "\t%.8g" "\t%.8g\n",
			0.0, 0.0, start_speed, start_speed, 0.0, start_speed / Dimens_v);
	fflush(stdout);
	for (k = 1; k < Nb; k++) {
		fprintf(stdout,
				"%.5g" "\t%.5g" "\t%.8g" "\t%.8g" "\t%.8g" "\t%.8g\n",
				k * d_time, k * d_time / Dimens_t,
				X[k] / N / k / d_time, U[k] / N, X[k] / N,
				U[k] / N / Dimens_v);
	}
#endif
	/* }}} */

	return 0;
}

/* Modeline for text editor Vim:
 * vim: ts=4
 * vim600: fdm=marker fdc=3
*/