/*
* 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
*/
Platon Group <platon@platon.org> http://platon.org/
|