Platon Technologies
not logged in Login Registration
EnglishSlovak
open source software development celebrating 10 years of open source development! Friday, March 29, 2024

File: [Platon] / doc / diplomova-praca-rajo / hotGas.c (download)

Revision 1.1, Thu Jul 24 17:50:07 2003 UTC (20 years, 8 months ago) by rajo

Diploma Thesis of Lubomir Host.
Title: Metoda Monte Carlo vo fyzike nizkoteplotnej plazmy. (Slovak language)

/*
 * 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://www.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
*/

Platon Group <platon@platon.org> http://platon.org/
Copyright © 2002-2006 Platon Group
Site powered by Metafox CMS
Go to Top