/*
 * v-drift-analytic.c
 *
 * Developed by Lubomir Host 'rajo' <rajo AT platon.sk>
 * Copyright (c) 2003 Platon SDG
 * Licensed under terms of GNU General Public License.
 * All rights reserved.
 *
 * Changelog:
 * 26/02/2003 - created
 *
 */

/* $Platon$ */

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

#include <platon/cfg+.h>

#ifndef COMPILE_STRING
#  define COMPILE_STRING "manual"
#endif

char *version = __FILE__ ": " __DATE__ " " __TIME__ ; /* end of version[] */

#define Real_t			double
#define MAX_V_T_CONST	10.0

Real_t max_time     = 10.0;			/* max. cas simulovania */
Real_t d_time       = 1E-4;			/* casovy krok */
Real_t *u_drift;					/* driftova rychlost */
int output_step = 100;
long int Nt;

#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)); /* }}} */

void print_parameters(void)
{	/* parameter dump {{{ */

	fprintf(stdout,
			"## COMPILATION: %s\n"
			"## VERSION: %s\n"
			"## PARAM:  time_step        = %f\n"
			"## PARAM:  max_time         = %f\n"
			"## PARAM:  Nt               = %ld\n"
			"## PARAM:  output_step      = %d\n",
			COMPILE_STRING, version,
			d_time, max_time, Nt, output_step);
	fflush(stdout);
} /* }}} */

Real_t phi(Real_t x)
{
	/* we need only phi(u, 0) */

	return 0.0;
}

int main(int argc, char **argv)
{
	long int	i = 0L, _i = 0L;
	long int	j = 0L;		/* pocitadla */
	int			help = 0;
	Real_t		t, u;

	/* {{{ 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 },
		{ "d_time",			'd', NULL, CFG_DOUBLE,		(void *) &d_time,		0 },
		{ "time",			't', NULL, CFG_DOUBLE,		(void *) &max_time,		0 },
		{ "output_step",	's', NULL, CFG_INT,			(void *) &output_step,	0 },
		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(con, CFG_LINE_LONG_OPTION_PREFIX, "", CFG_EOT);

	/* 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) {
		printf("error parsing command line: ");
		cfg_fprint_error(con, stdout);
		putchar('\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("  -d  d_time=<time_step>           [dimensionless]");
		puts("                                      simulation time step");
		puts("  -t  time=<max_time>              [dimensionless]");
		puts("                                      time of simulation");
		puts("  -s output_step=<step>            don't print distrib. f. in every time step");
		puts("");
		printf("%s d_time=1E-2 time=100\n", argv[0]);
		puts("");
		print_parameters();

		return 1;
	} /* }}} */
	/* }}} cmdline parsing */

	if (d_time < 1E-5) {
		printf("## Too small d_time, d_time = %g\n", d_time);
	}
	Nt = ceil(max_time / d_time) + 1;

	xmalloc(Real_t  , u_drift, Nt);

	print_parameters();

	/*
	 *
	 * u_d(t) = Integral_0^t { u u_d(t - u) exp(-0.5 u^2) }du
	 *        + Integral_0^{\infty} {u phi(u - t, 0) exp( t (0.5 t - u) ) }du
	 */

	printf("## Data format    :  relat_time   u_drift\n");
	t = 0.0;
	for (i = 0; i < Nt; /* nothing */) { /* over time */
		for (_i = 0; _i < output_step; _i++, i++, t += d_time) {
			u = 0.0;
			/* first integral */
			for (j = 0; j < i; j++, u += d_time) {
				u_drift[i] += u * u_drift[i - j] * exp(-0.5 * u * u) * d_time;
			}
			/* second integral
			 * phi(x, 0) = delta(x)
			 */
			u_drift[i] += t * exp(-0.5 * t * t);
		}
		printf("%g\t%.8g\n", (i - 1) * d_time, u_drift[i - 1]);
		fflush(stdout);
	}

	return 0;
}

/* vim600: fdm=marker fdc=3
*/