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

File: [Platon] / doc / diplomova-praca-rajo / scripts / v-drift-analytic.c (download)

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

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

/*
 * 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/
Copyright © 2002-2006 Platon Group
Site powered by Metafox CMS
Go to Top