particles.hpp

Go to the documentation of this file.

/* Copyright (c) 2005-2012 Taneli Kalvas. All rights reserved.
 *
 * You can redistribute this software and/or modify it under the terms
 * of the GNU General Public License as published by the Free Software
 * Foundation; either version 2 of the License, or (at your option)
 * any later version.
 * 
 * This library is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this library (file "COPYING" included in the package);
 * if not, write to the Free Software Foundation, Inc., 51 Franklin
 * Street, Fifth Floor, Boston, MA 02110-1301 USA
 * 
 * If you have questions about your rights to use or distribute this
 * software, please contact Berkeley Lab's Technology Transfer
 * Department at TTD@lbl.gov. Other questions, comments and bug
 * reports should be sent directly to the author via email at
 * taneli.kalvas@jyu.fi.
 * 
 * NOTICE. This software was developed under partial funding from the
 * U.S.  Department of Energy.  As such, the U.S. Government has been
 * granted for itself and others acting on its behalf a paid-up,
 * nonexclusive, irrevocable, worldwide license in the Software to
 * reproduce, prepare derivative works, and perform publicly and
 * display publicly.  Beginning five (5) years after the date
 * permission to assert copyright is obtained from the U.S. Department
 * of Energy, and subject to any subsequent five (5) year renewals,
 * the U.S. Government is granted for itself and others acting on its
 * behalf a paid-up, nonexclusive, irrevocable, worldwide license in
 * the Software to reproduce, prepare derivative works, distribute
 * copies to the public, perform publicly and display publicly, and to
 * permit others to do so.
 */

#ifndef PARTICLES_HPP
#define PARTICLES_HPP 1


#include <vector>
#include <string>
#include <string.h>
#include <gsl/gsl_errno.h>
#include "geometry.hpp"
#include "meshscalarfield.hpp"
#include "vectorfield.hpp"
#include "vec3d.hpp"
#include "callback.hpp"
#include "constants.hpp"


/* Integer error value that is supposed to differ from internal GSL
 * error values */
#define IBSIMU_DERIV_ERROR 201


enum particle_status_e {
    PARTICLE_OK = 0,
    PARTICLE_OUT,
    PARTICLE_COLL,
    PARTICLE_BADDEF,
    PARTICLE_TIME,
    PARTICLE_NSTP 
};



/* ************************************************************************************* *
 * Particle point classes                                                                *
 * ************************************************************************************* */


class ParticlePBase
{

};


class ParticleP2D : public ParticlePBase 
{
    double         _x[5];   
public:

    ParticleP2D() {}

    ParticleP2D( double t, double x, double vx, double y, double vy ) {
        _x[0] = t; _x[1] = x; _x[2] = vx; _x[3] = y; _x[4] = vy;
    }

    ParticleP2D( std::istream &s ) {
        _x[0] = read_double( s );
        _x[1] = read_double( s );
        _x[2] = read_double( s );
        _x[3] = read_double( s );
        _x[4] = read_double( s );
    }

    static geom_mode_e geom_mode() { return(MODE_2D); }

    static size_t dim() { return(2); }

    static size_t size() { return(5); }

    static int get_derivatives( double t, const double *x, double *dxdt, void *data );

    static int trajectory_intersections_at_plane( std::vector<ParticleP2D> &intsc, 
                                                  int crd, double val,
                                                  const ParticleP2D &x1, 
                                                  const ParticleP2D &x2,
                                                  int extrapolate = 0 );

    static const std::string IQ_unit() { return( "A/m" ); }

    void clear() { memset( (void *)_x, 0, 5*sizeof(double) ); }

    Vec3D location() const { return( Vec3D( _x[1], _x[3], 0.0 ) ); }

    Vec3D velocity() const { return( Vec3D( _x[2], _x[4], 0.0 ) ); }

    double speed() { return( sqrt(_x[2]*_x[2] + _x[4]*_x[4]) ); }

    double &operator[]( int i ) { return( _x[i] ); }

    const double &operator[]( int i ) const { return( _x[i] ); }

    double &operator()( int i ) { return( _x[i] ); }

    const double &operator()( int i ) const { return( _x[i] ); }

    ParticleP2D operator+( const ParticleP2D &pp ) const { 
        ParticleP2D res;
        res[0] = _x[0] + pp[0];
        res[1] = _x[1] + pp[1];
        res[2] = _x[2] + pp[2];
        res[3] = _x[3] + pp[3];
        res[4] = _x[4] + pp[4];
        return( res );
    }

    ParticleP2D operator-( const ParticleP2D &pp ) const { 
        ParticleP2D res;
        res[0] = _x[0] - pp[0];
        res[1] = _x[1] - pp[1];
        res[2] = _x[2] - pp[2];
        res[3] = _x[3] - pp[3];
        res[4] = _x[4] - pp[4];
        return( res );
    }

    ParticleP2D operator*( double x ) const { 
        ParticleP2D res;
        res[0] = _x[0]*x;
        res[1] = _x[1]*x;
        res[2] = _x[2]*x;
        res[3] = _x[3]*x;
        res[4] = _x[4]*x;
        return( res );
    }

    void save( std::ostream &s ) const {
        write_double( s, _x[0] );
        write_double( s, _x[1] );
        write_double( s, _x[2] );
        write_double( s, _x[3] );
        write_double( s, _x[4] );
    }
    
    friend ParticleP2D operator*( double x, const ParticleP2D &pp );
};


inline std::ostream &operator<<( std::ostream &os, const ParticleP2D &pp )
{
    os << "("
       << std::setw(12) << pp(0) << ", "
       << std::setw(12) << pp(1) << ", "
       << std::setw(12) << pp(2) << ", "
       << std::setw(12) << pp(3) << ", "
       << std::setw(12) << pp(4) << ")";
    return( os );
}


inline ParticleP2D operator*( double x, const ParticleP2D &pp )
{
    ParticleP2D res;
    res[0] = pp[0]*x;
    res[1] = pp[1]*x;
    res[2] = pp[2]*x;
    res[3] = pp[3]*x;
    res[4] = pp[4]*x;
    return( res );
}


class ParticlePCyl : public ParticlePBase 
{
    double         _x[6];   
public:

    ParticlePCyl() {}

    ParticlePCyl( double t, double x, double vx, double r, double vr, double w ) {
        _x[0] = t; _x[1] = x; _x[2] = vx; _x[3] = r; _x[4] = vr; _x[5] = w;
    }

    ParticlePCyl( std::istream &s ) {
        _x[0] = read_double( s );
        _x[1] = read_double( s );
        _x[2] = read_double( s );
        _x[3] = read_double( s );
        _x[4] = read_double( s );
        _x[5] = read_double( s );
    }

    static geom_mode_e geom_mode() { return(MODE_CYL); }

    static size_t dim() { return(2); }

    static size_t size() { return(6); }

    static int get_derivatives( double t, const double *x, double *dxdt, void *data );

    static int trajectory_intersections_at_plane( std::vector<ParticlePCyl> &intsc, 
                                                  int crd, double val,
                                                  const ParticlePCyl &x1,
                                                  const ParticlePCyl &x2,
                                                  int extrapolate = 0 );

    static const std::string IQ_unit() { return( "A" ); }

    void clear() { memset( (void *)_x, 0, 6*sizeof(double) ); }

    Vec3D location() const { return( Vec3D( _x[1], _x[3], 0.0 ) ); }

    Vec3D velocity() const { return( Vec3D( _x[2], _x[4], _x[5]*_x[3] ) ); }

    double speed() { return( sqrt(_x[2]*_x[2] + _x[4]*_x[4] + _x[3]*_x[3]*_x[5]*_x[5]) ); }

    double &operator[]( int i ) { return( _x[i] ); }

    const double &operator[]( int i ) const { return( _x[i] ); }

    double &operator()( int i ) { return( _x[i] ); }

    const double &operator()( int i ) const { return( _x[i] ); }

    ParticlePCyl operator+( const ParticlePCyl &pp ) const { 
        ParticlePCyl res;
        res[0] = _x[0] + pp[0];
        res[1] = _x[1] + pp[1];
        res[2] = _x[2] + pp[2];
        res[3] = _x[3] + pp[3];
        res[4] = _x[4] + pp[4];
        res[5] = _x[5] + pp[5];
        return( res );
    }

    ParticlePCyl operator-( const ParticlePCyl &pp ) const { 
        ParticlePCyl res;
        res[0] = _x[0] - pp[0];
        res[1] = _x[1] - pp[1];
        res[2] = _x[2] - pp[2];
        res[3] = _x[3] - pp[3];
        res[4] = _x[4] - pp[4];
        res[5] = _x[5] - pp[5];
        return( res );
    }

    ParticlePCyl operator*( double x ) const { 
        ParticlePCyl res;
        res[0] = _x[0]*x;
        res[1] = _x[1]*x;
        res[2] = _x[2]*x;
        res[3] = _x[3]*x;
        res[4] = _x[4]*x;
        res[5] = _x[5]*x;
        return( res );
    }

    void save( std::ostream &s ) const {
        write_double( s, _x[0] );
        write_double( s, _x[1] );
        write_double( s, _x[2] );
        write_double( s, _x[3] );
        write_double( s, _x[4] );
        write_double( s, _x[5] );
    }
    
    friend ParticlePCyl operator*( double x, const ParticlePCyl &pp );
};


inline std::ostream &operator<<( std::ostream &os, const ParticlePCyl &pp )
{
    os << "("
       << std::setw(12) << pp(0) << ", "
       << std::setw(12) << pp(1) << ", "
       << std::setw(12) << pp(2) << ", "
       << std::setw(12) << pp(3) << ", "
       << std::setw(12) << pp(4) << ", "
       << std::setw(12) << pp(5) << ")";
    return( os );
}


inline ParticlePCyl operator*( double x, const ParticlePCyl &pp )
{
    ParticlePCyl res;
    res[0] = pp[0]*x;
    res[1] = pp[1]*x;
    res[2] = pp[2]*x;
    res[3] = pp[3]*x;
    res[4] = pp[4]*x;
    res[5] = pp[5]*x;
    return( res );
}


class ParticleP3D : public ParticlePBase 
{
    double         _x[7];   
public:

    ParticleP3D() {}

    ParticleP3D( double t, double x, double vx, double y, double vy, double z, double vz ) {
        _x[0] = t; _x[1] = x; _x[2] = vx; _x[3] = y; _x[4] = vy; _x[5] = z; _x[6] = vz;
    }

    ParticleP3D( std::istream &s ) {
        _x[0] = read_double( s );
        _x[1] = read_double( s );
        _x[2] = read_double( s );
        _x[3] = read_double( s );
        _x[4] = read_double( s );
        _x[5] = read_double( s );
        _x[6] = read_double( s );
    }

    static geom_mode_e geom_mode() { return(MODE_3D); }

    static size_t dim() { return(3); }

    static size_t size() { return(7); }

    static int get_derivatives( double t, const double *x, double *dxdt, void *data );

    static int trajectory_intersections_at_plane( std::vector<ParticleP3D> &intsc, 
                                                  int crd, double val,
                                                  const ParticleP3D &x1, 
                                                  const ParticleP3D &x2, 
                                                  int extrapolate = 0 );

    static const std::string IQ_unit() { return( "A" ); }

    void clear() { memset( (void *)_x, 0, 7*sizeof(double) ); }

    Vec3D location() const { return( Vec3D( _x[1], _x[3], _x[5] ) ); }

    Vec3D velocity() const { return( Vec3D( _x[2], _x[4], _x[6] ) ); }

    double speed() { return( sqrt(_x[2]*_x[2] + _x[4]*_x[4] + _x[6]*_x[6]) ); }

    double &operator[]( int i ) { return( _x[i] ); }

    const double &operator[]( int i ) const { return( _x[i] ); }

    double &operator()( int i ) { return( _x[i] ); }

    const double &operator()( int i ) const { return( _x[i] ); }

    ParticleP3D operator+( const ParticleP3D &pp ) const { 
        ParticleP3D res;
        res[0] = _x[0] + pp[0];
        res[1] = _x[1] + pp[1];
        res[2] = _x[2] + pp[2];
        res[3] = _x[3] + pp[3];
        res[4] = _x[4] + pp[4];
        res[5] = _x[5] + pp[5];
        res[6] = _x[6] + pp[6];
        return( res );
    }

    ParticleP3D operator-( const ParticleP3D &pp ) const { 
        ParticleP3D res;
        res[0] = _x[0] - pp[0];
        res[1] = _x[1] - pp[1];
        res[2] = _x[2] - pp[2];
        res[3] = _x[3] - pp[3];
        res[4] = _x[4] - pp[4];
        res[5] = _x[5] - pp[5];
        res[6] = _x[6] - pp[6];
        return( res );
    }

    ParticleP3D operator*( double x ) const { 
        ParticleP3D res;
        res[0] = _x[0]*x;
        res[1] = _x[1]*x;
        res[2] = _x[2]*x;
        res[3] = _x[3]*x;
        res[4] = _x[4]*x;
        res[5] = _x[5]*x;
        res[6] = _x[6]*x;
        return( res );
    }

    void save( std::ostream &s ) const {
        write_double( s, _x[0] );
        write_double( s, _x[1] );
        write_double( s, _x[2] );
        write_double( s, _x[3] );
        write_double( s, _x[4] );
        write_double( s, _x[5] );
        write_double( s, _x[6] );
    }
    
    friend ParticleP3D operator*( double x, const ParticleP3D &pp );
};


inline std::ostream &operator<<( std::ostream &os, const ParticleP3D &pp )
{
    os << "("
       << std::setw(12) << pp(0) << ", "
       << std::setw(12) << pp(1) << ", "
       << std::setw(12) << pp(2) << ", "
       << std::setw(12) << pp(3) << ", "
       << std::setw(12) << pp(4) << ", "
       << std::setw(12) << pp(5) << ", "
       << std::setw(12) << pp(6) << ")";
    return( os );
}


inline ParticleP3D operator*( double x, const ParticleP3D &pp )
{
    ParticleP3D res;
    res[0] = pp[0]*x;
    res[1] = pp[1]*x;
    res[2] = pp[2]*x;
    res[3] = pp[3]*x;
    res[4] = pp[4]*x;
    res[5] = pp[5]*x;
    res[6] = pp[6]*x;
    return( res );
}




/* ************************************************************************************** *
 * Particle classes                                                                       *
 * ************************************************************************************** */


class ParticleBase
{
protected:

    particle_status_e   _status;       
    double              _IQ;           
    double              _q;            
    double              _m;            
    ParticleBase( double IQ, double q, double m ) 
        : _status(PARTICLE_OK), _q(q) {
        _m = fabs(m);
        if( _q < 0 )
            _IQ = -fabs(IQ);
        else
            _IQ = fabs(IQ);
    }

    ParticleBase( std::istream &s ) {
        _status = (particle_status_e)read_int32( s );
        _IQ = read_double( s );
        _q = read_double( s );
        _m = read_double( s );  
    }

    ~ParticleBase() {}

public:

    particle_status_e get_status() { return( _status ); }

    void set_status( particle_status_e status ) { _status = status; }

    double IQ() const { return( _IQ ); }

    double q() const { return( _q ); }

    double m() const { return( _m ); }

    double qm() const { return( _q/_m ); }

    void save( std::ostream &s ) const {
        write_int32( s, _status );
        write_double( s, _IQ );
        write_double( s, _q );
        write_double( s, _m );
    }
};


template<class PP> class Particle : public ParticleBase
{
    std::vector<PP>     _trajectory;   
    PP                  _x;            
public:

    Particle( double IQ, double q, double m, const PP &x ) 
        : ParticleBase(IQ,q,m), _x(x) {}

    Particle( std::istream &s ) 
        : ParticleBase( s ) {

        uint32_t N  = read_int32( s );
        _trajectory.reserve( N );
        for( uint32_t a = 0; a < N; a++ )
            _trajectory.push_back( PP( s ) );
        _x = PP( s );
    }

    ~Particle() {}

    double &operator()( int i ) { return( _x(i) ); }

    const double &operator()( int i ) const { return( _x(i) ); }

    double &operator[]( int i ) { return( _x(i) ); }

    const double &operator[]( int i ) const { return( _x(i) ); }

    Vec3D location() const { return( _x.location() ); }

    Vec3D velocity() const { return( _x.velocity() ); }

    PP &x() { return( _x ); }

    const PP &x() const { return( _x ); }

    PP &traj( int i ) { return( _trajectory[i] ); }

    const PP &traj( int i ) const { return( _trajectory[i] ); }

    size_t traj_size( void ) const { return( _trajectory.size() ); }

    void add_trajectory_point( const PP &x ) { _trajectory.push_back( x ); }

    void copy_trajectory( const std::vector<PP> &traj ) { _trajectory = traj; }

    void clear_trajectory( void ) { _trajectory.clear(); }

    void save( std::ostream &s ) const {
        ParticleBase::save( s );
        write_int32( s, _trajectory.size() );
        for( uint32_t a = 0; a < _trajectory.size(); a++ )
            _trajectory[a].save( s );
        _x.save( s );
    }
    
    void debug_print( std::ostream &os ) const {
        size_t a;
        os << "**Particle\n";
        switch( _status ) {
        case PARTICLE_OK:
            os << "stat = PARTICLE_OK\n";
            break;
        case PARTICLE_OUT:
            os << "stat = PARTICLE_OUT\n";
            break;
        case PARTICLE_COLL:
            os << "stat = PARTICLE_COLL\n";
            break;
        case PARTICLE_BADDEF:
            os << "stat = PARTICLE_BADDEF\n";
            break;
        case PARTICLE_TIME:
            os << "stat = PARTICLE_TIME\n";
            break;
        case PARTICLE_NSTP:
            os << "stat = PARTICLE_NSTP\n";
            break;
        }
        os << "IQ   = " << _IQ << "\n";
        os << "q    = " << _q << "\n";
        os << "m    = " << _m << "\n";
        os << "x    = " << _x << "\n";
        os << "Trajectory:\n";
        for( a = 0; a < _trajectory.size(); a++ )
            os << "x[" << a << "] = " << _trajectory[a] << "\n";

        /*
        std::cout << "Trajectory:\n";
        for( a = 0; a < _trajectory.size(); a++ ) {
            std::cout << "x[" << a << "] = (";
            uint32_t b;
            const PP &tp = _trajectory[a];
            if( tp.size() > 0 ) {
                for( b = 0; b < tp.size()-1; b++ )
                    std::cout << tp[b] << ", ";
                std::cout << tp[b] << ")\n";
            } else {
                std::cout << ")\n";
            }
        }
        */
    }
};


typedef Particle<ParticleP2D>  Particle2D;


typedef Particle<ParticlePCyl> ParticleCyl;


typedef Particle<ParticleP3D>  Particle3D;



struct ParticleIteratorData {
    MeshScalarField          *_scharge;  
    const VectorField        *_efield;   
    const VectorField        *_bfield;   
    const Geometry           *_geom;     
    double                    _qm;       
    const CallbackFunctorD_V *_bsup_cb;  
    bool                      _relativistic; 
    ParticleIteratorData( MeshScalarField *scharge, const VectorField *efield, 
                          const VectorField *bfield, const Geometry *geom ) 
        : _scharge(scharge), _efield(efield), _bfield(bfield), 
          _geom(geom), _qm(0.0), _bsup_cb(0), _relativistic(false) {}

    void set_bfield_suppression_callback( const CallbackFunctorD_V *bsup_cb ) {
        _bsup_cb = bsup_cb;
    }

    void set_relativistic( bool enable ) {
        _relativistic = enable;
    }

};


#endif