source: ZHANGProjects/ICOSIM/CPP/trunk/source/Particle.cc @ 2

#include <iostream>
#include <vector>
#include <string>
#include <cmath>
#include "Particle.h"
using namespace std;


//==================================Constructors, destructor====================================

Particle::Particle(double x, double dx, double y, double dy, double deltaE, double t, double mass, double charge, double moment)
    : Ap0(mass), Zp0(charge), dpoporiginal(moment)
{
    coordonnees[0][0] = x;
    coordonnees[0][1] = dx;
    coordonnees[0][2] = y;
    coordonnees[0][3] = dy;
    coordonnees[0][4] = deltaE;
    coordonnees[0][5] = t;
    inabs = 1;
    coordonnees[1][0] = 0;
    coordonnees[1][1] = 0;
    coordonnees[1][2] = 0;
    coordonnees[1][3] = 0;
    coordonnees[1][4] = 0;
    coordonnees[1][5] = 0;
}



Particle::Particle(const Particle& obj)
    : Ap0(obj.Ap0),
      Zp0(obj.Zp0),
      in(obj.in),
      inabs(obj.inabs),
      nrevhitp(obj.nrevhitp),
      dpoporiginal(obj.dpoporiginal),
      dt(obj.dt)

{
    coordonnees[0][0] = obj.coordonnees[0][0];
    coordonnees[0][1] = obj.coordonnees[0][1];
    coordonnees[0][2] = obj.coordonnees[0][2];
    coordonnees[0][3] = obj.coordonnees[0][3];
    coordonnees[0][4] = obj.coordonnees[0][4];
    coordonnees[0][5] = obj.coordonnees[0][5];
    identification = obj.getidentification();
}



void Particle::afficheCoordonnees()
{

    for (int j(0); j < 2; ++j) {
        for (int i(0); i < 6; ++i) {
            cout << "Here is the " << i + 1 << "th coordinate of the particle at the " << j + 1 << "th interpolation point: " << coordonnees[j][i] << endl;
        }
    }

}



void Particle::afficheFirstCoordonnees()
{
    for (int i(0); i < 6; ++i) {
        cout << "Here is the " << i + 1 << "th coordinate of the particle:" << coordonnees[0][i] << endl;
    }
}


int Particle::incog(double a, double x, double& gin, double& gim, double& gip)
{

    double xam, r, s, ga, t0;
    int k;

    if ((a < 0.0) || (x < 0)) {
        return 1;
    }
    xam = -x + a * log(x);
    if ((xam > 700) || ( a > 170.0)) {
        return 1;
    }
    if (x == 0.0) {
        gin = 0.0;
        gim = exp(gamma(a));
        gip = 0.0;
        return 0;
    }
    if (x <= 1.0 + a) {
        s = 1.0 / a;
        r = s;
        for (k = 1; k <= 60; k++) {
            r *= x / (a + k);
            s += r;
            if (fabs(r / s) < 1e-15) {
                break;
            }
        }
        gin = exp(xam) * s;
        ga = exp(gamma(a));
        gip = gin / ga;
        gim = ga - gin;
    } else {
        t0 = 0.0;
        for (k = 60; k >= 1; k--) {
            t0 = (k - a) / (1.0 + k / (x + t0));
        }
        gim = exp(xam) / (x + t0);
        ga = exp(gamma(a));
        gin = ga - gim;
        gip = 1.0 - gim / ga;
    }
    return 0;
}


int Particle::getidentification() const
{
    return identification;
}


void Particle::setidentification(int id)
{
    identification = id;
}


Particle& Particle::operator=(Particle const& source)
{
    Ap0 = source.Ap0;
    Zp0 = source.Zp0;
    inabs = source.inabs;
    for (int i(0); i < 2; ++i) {
        for (int j(0); j < 6; ++j) {
            coordonnees[i][j] = source.coordonnees[i][j];
        }
    }
    nrevhitp = source.nrevhitp;
    in = source.in;
    identification = source.getidentification();

    return *this;
}


double Particle::random()
{
    double num;
    num = rand() / double(RAND_MAX);
    return -num;
    //return -0.5;//this lign can be uncomment to fix the randomness for tests (see the file randomness.txt)
}


void Particle::genpartdist(const int& i0, const int& n, const string& type, const double& r1r2skin, const double& emx, const double& emy, const double& sigdpp, const double& bx, const double& ax, const double& dx, const double& dpx, const double& by, const double& ay, const double& dy, const double& dpy, const double& nsigi)
{

    int nh;
    double nr, sum, Rx, Ry;
    vector < vector < double > > h1, h;
    vector < double > hsum1, hsum, xh, yh, xsh, ysh;

    nh = (int)ceil(n / 0.29);
    nr = 0;

    //after this loop:
    //h is 4 x nr, where nr>=n
    //length of all columns is less than 1
    //elements of h uniformly distributed between -1 and 1

    while (nr < n) {

        hsum.clear();
        hsum1.clear();
        h.clear();
        h1.clear();

        for (int i(0); i < 4; ++i) {
            h.push_back(vector <double> (nh));
            for (int j(0); j < nh; ++j) {
                h[i][j] = (-2 * random() - 1); //uniformly distributed numbers in [-1,1]
            }
        }

        for (int k(0); k < nh; ++k) {
            sum = 0;
            for (int l(0); l < 4; ++l) {
                sum = sum + (h[k][l] * h[k][l]);
            }

            hsum.push_back(sum);
        }

        int incr(0);

        for (int i(0); i < nh; ++i) {
            if (hsum[i] < 1) {
                h1.push_back(vector < double > (4));
                for (int l(0); l < 4; ++l) {
                    h1[incr][l] = h[i][l];
                }
                hsum1.push_back(hsum[i]);
                ++incr;
            }
        }

        nr = hsum1.size();

    }//end while

    for (int i(0); i < nr; ++i) {
        hsum1[i] = sqrt(hsum1[i]);
    }

    h.clear();
    hsum.clear();

    for (int j(0); j < n; ++j) {
        h.push_back(h1[j]);
        hsum.push_back(hsum1[j]);
    }

    if (type == "kv") {
        for (int i(0); i < n; ++i) {
            for (int j(0); j < 4; ++j) {
                h[i][j] = h[i][j] / hsum[i]; //normalizes the columns of h
            }
        }

        Rx = sqrt(emx);
        Ry = sqrt(emy);

        for (int k(0); k < n; ++k) {
            xh.push_back(Rx * h[k][0]);
            yh.push_back(Ry * h[k][1]);
            xsh.push_back(Rx * h[k][2]);
            ysh.push_back(Ry * h[k][3]);
        }

    } //end if

    else if (type == "waterbag") {
        //xh uniformly distributed on [-sqrt(emx) sqrt(emx)]
        //similarly yh, xsh, ysh
        //BUT such that length([xh/sqrt(emx) ...])>1
        Rx = sqrt(emx);
        Ry = sqrt(emy);

        for (int k(0); k < n; ++k) {
            xh.push_back(Rx * h[k][0]);
            yh.push_back(Ry * h[k][1]);
            xsh.push_back(Rx * h[k][2]);
            ysh.push_back(Ry * h[k][3]);
        }
    }

    else if (type == "r1r2") {
        double R2 , rv;
        R2 = sqrt(r1r2skin);
        rv = pow(((R2 * R2 * R2 * R2 - 1) * (-random() + 1 / (R2 * R2 * R2 * R2 - 1))), 0.25);

        for (int i(0); i < n; ++i) {
            for (int j(0); j < 4; ++j) {
                h[i][j] = rv * h[i][j] / hsum[i];
            }
        }

        Rx = sqrt(emx);
        Ry = sqrt(emy);

        for (int k(0); k < n; ++k) {
            xh.push_back(Rx * h[k][0]);
            yh.push_back(Ry * h[k][1]);
            xsh.push_back(Rx * h[k][2]);
            ysh.push_back(Ry * h[k][3]);
        }
    } else if (type == "gauss") {

        int p(3);
        vector <double> rh, fh;
        int step, temp, count;
        double gin1, gim1, gip1, gin2, gim2, gip2, randnumber, rv;

        rh.push_back(0);

        step = (int)(nsigi / 0.001) + 1;

        for (int u(0); u < step; ++u) {
            rh.push_back(rh[rh.size() - 1] + 0.001);
        }
        temp = incog((1 + p) / 2, 0.5 * nsigi * nsigi, gin2, gim2, gip2);

        for (int u(0); u <= step; ++u) {
            temp = incog((1 + p) / 2, 0.5 * rh[u] * rh[u], gin1, gim1, gip1);
            fh.push_back(gin1 / gin2);
        }

        randnumber = -random();

        count = 0;
        for (int k(0); k < step; ++k) {
            if ((randnumber > fh[k]) && (randnumber < fh[k + 1])) {
                break;
            } else {
                ++ count;
            }
        }

        rv = interp1(fh[count], rh[count], fh[count + 1], rh[count + 1], randnumber);

        for (int i(0); i < n; ++i) {
            for (int j(0); j < 4; ++j) {
                h[i][j] = rv * h[i][j] / hsum[i];
            }
        }

        Rx = sqrt(emx);
        Ry = sqrt(emy);

        for (int k(0); k < n; ++k) {
            xh.push_back(Rx * h[k][0]);
            yh.push_back(Ry * h[k][1]);
            xsh.push_back(Rx * h[k][2]);
            ysh.push_back(Ry * h[k][3]);
        }

    } else {

        cerr << "GENPARTDIST ERROR: Distribution type " << type << " not implemented !" << endl;
    }


    //generating a normal random variable using Box Muller method


    this->coordonnees[0][4] = sigdpp * sqrt(-2 * log(-random())) * cos(2 * M_PI * (-random()));
    this->coordonnees[0][0] = sqrt(bx) * xh[0] + dx * this->coordonnees[0][4];
    this->coordonnees[0][1] = -ax / sqrt(bx) * xh[0] + xsh[0] / sqrt(bx) + dpx * this->coordonnees[0][4];
    this->coordonnees[0][2] = sqrt(by) * yh[0] + dy * this->coordonnees[0][4];
    this->coordonnees[0][3] = -ay / sqrt(by) * yh[0] + ysh[0] / sqrt(by) + dpy * this->coordonnees[0][4];
    this->dt = 0;
    this->dpoporiginal = this->coordonnees[0][4];

    //afficheFirstCoordonnees();


}


double Particle::interp1(double x1, double y1, double x2, double y2, double x)
{

    double d1(x1 - x2);
    double d2(y1 - y2);

    double y;

    y = (x - x1) * (d2 / d1) + y1;

    return y;
}