%
% '07-literatura.tex' created: Fri Mar 28 16:44:29 CET 2003
%
% 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.
%
% $Platon: doc/diplomova-praca-rajo/07-prilohy.tex,v 1.1 2003/07/24 17:48:13 rajo Exp $
\def\appendixname{Príloha}
\appendix
\setAppendixHeader
\chapter{Zdrojový kód programu}
% don't remove 'hotGas.c' from next line, `insert-C-source` depends on it!
\beginCsource{hotGas.c}
/*
* 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
*/
\endCsource
% vim: ts=4
% vim600: fdl=0 fdm=marker fdc=3
Platon Group <platon@platon.org> http://platon.org/
|