particles.hpp

Go to the documentation of this file.
00001 
00005 /* Copyright (c) 2005-2009 Taneli Kalvas. All rights reserved.
00006  *
00007  * You can redistribute this software and/or modify it under the terms
00008  * of the GNU General Public License as published by the Free Software
00009  * Foundation; either version 2 of the License, or (at your option)
00010  * any later version.
00011  * 
00012  * This library is distributed in the hope that it will be useful, but
00013  * WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
00015  * General Public License for more details.
00016  * 
00017  * You should have received a copy of the GNU General Public License
00018  * along with this library (file "COPYING" included in the package);
00019  * if not, write to the Free Software Foundation, Inc., 51 Franklin
00020  * Street, Fifth Floor, Boston, MA 02110-1301 USA
00021  * 
00022  * If you have questions about your rights to use or distribute this
00023  * software, please contact Berkeley Lab's Technology Transfer
00024  * Department at TTD@lbl.gov. Other questions, comments and bug
00025  * reports should be sent directly to the author via email at
00026  * taneli.kalvas@jyu.fi.
00027  * 
00028  * NOTICE. This software was developed under partial funding from the
00029  * U.S.  Department of Energy.  As such, the U.S. Government has been
00030  * granted for itself and others acting on its behalf a paid-up,
00031  * nonexclusive, irrevocable, worldwide license in the Software to
00032  * reproduce, prepare derivative works, and perform publicly and
00033  * display publicly.  Beginning five (5) years after the date
00034  * permission to assert copyright is obtained from the U.S. Department
00035  * of Energy, and subject to any subsequent five (5) year renewals,
00036  * the U.S. Government is granted for itself and others acting on its
00037  * behalf a paid-up, nonexclusive, irrevocable, worldwide license in
00038  * the Software to reproduce, prepare derivative works, distribute
00039  * copies to the public, perform publicly and display publicly, and to
00040  * permit others to do so.
00041  */
00042 
00043 #ifndef PARTICLES_HPP
00044 #define PARTICLES_HPP 1
00045 
00046 
00047 #include <vector>
00048 #include <gsl/gsl_errno.h>
00049 #include "geometry.hpp"
00050 #include "scalarfield.hpp"
00051 #include "vectorfield.hpp"
00052 #include "efield.hpp"
00053 #include "vec3d.hpp"
00054 
00055 
00056 /* Integer error value that is supposed to diffed from internal GSL error values */
00057 #define IBSIMU_DERIV_ERROR 201
00058 
00059 
00067 enum particle_status_e {
00068     PARTICLE_OK = 0,
00069     PARTICLE_OUT,
00070     PARTICLE_COLL,
00071     PARTICLE_BADDEF,
00072     PARTICLE_TIME,
00073     PARTICLE_NSTP 
00074 };
00075 
00076 
00077 /* Atomic mass unit */
00078 #define MASS_U        1.66053873e-27
00079 /* Elementary charge */
00080 #define CHARGE_E      1.602176462e-19
00081 
00082 
00083 
00084 
00085 /* ************************************************************************************* *
00086  * Particle point classes                                                                *
00087  * ************************************************************************************* */
00088 
00089 
00094 class ParticlePBase
00095 {
00096 
00097 };
00098 
00099 
00105 class ParticleP2D : public ParticlePBase 
00106 {
00107     double         _x[5];   
00109 public:
00110 
00113     ParticleP2D() {}
00114 
00117     ParticleP2D( double t, double x, double vx, double y, double vy ) {
00118         _x[0] = t; _x[1] = x; _x[2] = vx; _x[3] = y; _x[4] = vy;
00119     }
00120 
00123     static geom_mode_e geom_mode() { return(MODE_2D); }
00124 
00127     static size_t dim() { return(2); }
00128 
00131     static size_t size() { return(5); }
00132 
00144     static int get_derivatives( double t, const double *x, double *dxdt, void *data );
00145 
00150     static int trajectory_intersections_at_plane( std::vector<ParticleP2D> &intsc, 
00151                                                   int crd, double val,
00152                                                   const ParticleP2D &x1, const ParticleP2D &x2 );
00153 
00156     Vec3D location() const { return( Vec3D( _x[1], _x[3], 0.0 ) ); }
00157 
00160     Vec3D velocity() const { return( Vec3D( _x[2], _x[4], 0.0 ) ); }
00161 
00164     double speed() { return( sqrt(_x[2]*_x[2] + _x[4]*_x[4]) ); }
00165 
00168     double &operator[]( int i ) { return( _x[i] ); }
00169 
00172     const double &operator[]( int i ) const { return( _x[i] ); }
00173 
00176     double &operator()( int i ) { return( _x[i] ); }
00177 
00180     const double &operator()( int i ) const { return( _x[i] ); }
00181 
00182     ParticleP2D operator+( const ParticleP2D &pp ) const { 
00183         ParticleP2D res;
00184         res[0] = _x[0] + pp[0];
00185         res[1] = _x[1] + pp[1];
00186         res[2] = _x[2] + pp[2];
00187         res[3] = _x[3] + pp[3];
00188         res[4] = _x[4] + pp[4];
00189         return( res );
00190     }
00191 
00192     ParticleP2D operator-( const ParticleP2D &pp ) const { 
00193         ParticleP2D res;
00194         res[0] = _x[0] - pp[0];
00195         res[1] = _x[1] - pp[1];
00196         res[2] = _x[2] - pp[2];
00197         res[3] = _x[3] - pp[3];
00198         res[4] = _x[4] - pp[4];
00199         return( res );
00200     }
00201 
00202     ParticleP2D operator*( double x ) const { 
00203         ParticleP2D res;
00204         res[0] = _x[0]*x;
00205         res[1] = _x[1]*x;
00206         res[2] = _x[2]*x;
00207         res[3] = _x[3]*x;
00208         res[4] = _x[4]*x;
00209         return( res );
00210     }
00211 
00212     friend ParticleP2D operator*( double x, const ParticleP2D &pp );
00213 };
00214 
00215 
00216 inline std::ostream &operator<<( std::ostream &os, const ParticleP2D &pp )
00217 {
00218     os << "("
00219        << std::setw(12) << pp(0) << ", "
00220        << std::setw(12) << pp(1) << ", "
00221        << std::setw(12) << pp(2) << ", "
00222        << std::setw(12) << pp(3) << ", "
00223        << std::setw(12) << pp(4) << ")";
00224     return( os );
00225 }
00226 
00227 
00228 inline ParticleP2D operator*( double x, const ParticleP2D &pp )
00229 {
00230     ParticleP2D res;
00231     res[0] = pp[0]*x;
00232     res[1] = pp[1]*x;
00233     res[2] = pp[2]*x;
00234     res[3] = pp[3]*x;
00235     res[4] = pp[4]*x;
00236     return( res );
00237 }
00238 
00239 
00245 class ParticlePCyl : public ParticlePBase 
00246 {
00247     double         _x[6];   
00249 public:
00250 
00253     ParticlePCyl() {}
00254 
00257     ParticlePCyl( double t, double x, double vx, double r, double vr, double w ) {
00258         _x[0] = t; _x[1] = x; _x[2] = vx; _x[3] = r; _x[4] = vr; _x[5] = w;
00259     }
00260 
00263     static geom_mode_e geom_mode() { return(MODE_CYL); }
00264 
00267     static size_t dim() { return(2); }
00268 
00271     static size_t size() { return(6); }
00272 
00289     static int get_derivatives( double t, const double *x, double *dxdt, void *data );
00290 
00295     static int trajectory_intersections_at_plane( std::vector<ParticlePCyl> &intsc, 
00296                                                   int crd, double val,
00297                                                   const ParticlePCyl &x1,
00298                                                   const ParticlePCyl &x2 );
00299 
00302     Vec3D location() const { return( Vec3D( _x[1], _x[3], 0.0 ) ); }
00303 
00306     Vec3D velocity() const { return( Vec3D( _x[2], _x[4], _x[5]*_x[3] ) ); }
00307 
00310     double speed() { return( sqrt(_x[2]*_x[2] + _x[4]*_x[4] + _x[3]*_x[3]*_x[5]*_x[5]) ); }
00311 
00314     double &operator[]( int i ) { return( _x[i] ); }
00315 
00318     const double &operator[]( int i ) const { return( _x[i] ); }
00319 
00322     double &operator()( int i ) { return( _x[i] ); }
00323 
00326     const double &operator()( int i ) const { return( _x[i] ); }
00327 
00328     ParticlePCyl operator+( const ParticlePCyl &pp ) const { 
00329         ParticlePCyl res;
00330         res[0] = _x[0] + pp[0];
00331         res[1] = _x[1] + pp[1];
00332         res[2] = _x[2] + pp[2];
00333         res[3] = _x[3] + pp[3];
00334         res[4] = _x[4] + pp[4];
00335         res[5] = _x[5] + pp[5];
00336         return( res );
00337     }
00338 
00339     ParticlePCyl operator-( const ParticlePCyl &pp ) const { 
00340         ParticlePCyl res;
00341         res[0] = _x[0] - pp[0];
00342         res[1] = _x[1] - pp[1];
00343         res[2] = _x[2] - pp[2];
00344         res[3] = _x[3] - pp[3];
00345         res[4] = _x[4] - pp[4];
00346         res[5] = _x[5] - pp[5];
00347         return( res );
00348     }
00349 
00350     ParticlePCyl operator*( double x ) const { 
00351         ParticlePCyl res;
00352         res[0] = _x[0]*x;
00353         res[1] = _x[1]*x;
00354         res[2] = _x[2]*x;
00355         res[3] = _x[3]*x;
00356         res[4] = _x[4]*x;
00357         res[5] = _x[5]*x;
00358         return( res );
00359     }
00360 
00361     friend ParticlePCyl operator*( double x, const ParticlePCyl &pp );
00362 };
00363 
00364 
00365 inline std::ostream &operator<<( std::ostream &os, const ParticlePCyl &pp )
00366 {
00367     os << "("
00368        << std::setw(12) << pp(0) << ", "
00369        << std::setw(12) << pp(1) << ", "
00370        << std::setw(12) << pp(2) << ", "
00371        << std::setw(12) << pp(3) << ", "
00372        << std::setw(12) << pp(4) << ", "
00373        << std::setw(12) << pp(5) << ")";
00374     return( os );
00375 }
00376 
00377 
00378 inline ParticlePCyl operator*( double x, const ParticlePCyl &pp )
00379 {
00380     ParticlePCyl res;
00381     res[0] = pp[0]*x;
00382     res[1] = pp[1]*x;
00383     res[2] = pp[2]*x;
00384     res[3] = pp[3]*x;
00385     res[4] = pp[4]*x;
00386     res[5] = pp[5]*x;
00387     return( res );
00388 }
00389 
00390 
00396 class ParticleP3D : public ParticlePBase 
00397 {
00398     double         _x[7];   
00400 public:
00401 
00404     ParticleP3D() {}
00405 
00408     ParticleP3D( double t, double x, double vx, double y, double vy, double z, double vz ) {
00409         _x[0] = t; _x[1] = x; _x[2] = vx; _x[3] = y; _x[4] = vy; _x[5] = z; _x[6] = vz;
00410     }
00411 
00414     static geom_mode_e geom_mode() { return(MODE_3D); }
00415 
00418     static size_t dim() { return(3); }
00419 
00422     static size_t size() { return(7); }
00423 
00437     static int get_derivatives( double t, const double *x, double *dxdt, void *data );
00438 
00443     static int trajectory_intersections_at_plane( std::vector<ParticleP3D> &intsc, 
00444                                                   int crd, double val,
00445                                                   const ParticleP3D &x1, const ParticleP3D &x2 );
00446 
00449     Vec3D location() const { return( Vec3D( _x[1], _x[3], _x[5] ) ); }
00450 
00453     Vec3D velocity() const { return( Vec3D( _x[2], _x[4], _x[6] ) ); }
00454 
00457     double speed() { return( sqrt(_x[2]*_x[2] + _x[4]*_x[4] + _x[6]*_x[6]) ); }
00458 
00461     double &operator[]( int i ) { return( _x[i] ); }
00462 
00465     const double &operator[]( int i ) const { return( _x[i] ); }
00466 
00469     double &operator()( int i ) { return( _x[i] ); }
00470 
00473     const double &operator()( int i ) const { return( _x[i] ); }
00474 
00475     ParticleP3D operator+( const ParticleP3D &pp ) const { 
00476         ParticleP3D res;
00477         res[0] = _x[0] + pp[0];
00478         res[1] = _x[1] + pp[1];
00479         res[2] = _x[2] + pp[2];
00480         res[3] = _x[3] + pp[3];
00481         res[4] = _x[4] + pp[4];
00482         res[5] = _x[5] + pp[5];
00483         res[6] = _x[6] + pp[6];
00484         return( res );
00485     }
00486 
00487     ParticleP3D operator-( const ParticleP3D &pp ) const { 
00488         ParticleP3D res;
00489         res[0] = _x[0] - pp[0];
00490         res[1] = _x[1] - pp[1];
00491         res[2] = _x[2] - pp[2];
00492         res[3] = _x[3] - pp[3];
00493         res[4] = _x[4] - pp[4];
00494         res[5] = _x[5] - pp[5];
00495         res[6] = _x[6] - pp[6];
00496         return( res );
00497     }
00498 
00499     ParticleP3D operator*( double x ) const { 
00500         ParticleP3D res;
00501         res[0] = _x[0]*x;
00502         res[1] = _x[1]*x;
00503         res[2] = _x[2]*x;
00504         res[3] = _x[3]*x;
00505         res[4] = _x[4]*x;
00506         res[5] = _x[5]*x;
00507         res[6] = _x[6]*x;
00508         return( res );
00509     }
00510 
00511     friend ParticleP3D operator*( double x, const ParticleP3D &pp );
00512 };
00513 
00514 
00515 inline std::ostream &operator<<( std::ostream &os, const ParticleP3D &pp )
00516 {
00517     os << "("
00518        << std::setw(12) << pp(0) << ", "
00519        << std::setw(12) << pp(1) << ", "
00520        << std::setw(12) << pp(2) << ", "
00521        << std::setw(12) << pp(3) << ", "
00522        << std::setw(12) << pp(4) << ", "
00523        << std::setw(12) << pp(5) << ", "
00524        << std::setw(12) << pp(6) << ")";
00525     return( os );
00526 }
00527 
00528 
00529 inline ParticleP3D operator*( double x, const ParticleP3D &pp )
00530 {
00531     ParticleP3D res;
00532     res[0] = pp[0]*x;
00533     res[1] = pp[1]*x;
00534     res[2] = pp[2]*x;
00535     res[3] = pp[3]*x;
00536     res[4] = pp[4]*x;
00537     res[5] = pp[5]*x;
00538     res[6] = pp[6]*x;
00539     return( res );
00540 }
00541 
00542 
00543 
00544 
00545 /* ************************************************************************************** *
00546  * Particle classes                                                                       *
00547  * ************************************************************************************** */
00548 
00549 
00554 class ParticleBase
00555 {
00556 protected:
00557 
00558     particle_status_e   _status;       
00559     double              _IQ;           
00570     double              _qm;           
00572     ParticleBase( double IQ, double q, double m ) 
00573         : _status(PARTICLE_OK), _qm(q/m) {
00574         if( _qm < 0 )
00575             _IQ = -fabs(IQ);
00576         else
00577             _IQ = fabs(IQ);
00578     }
00579 
00580     ~ParticleBase() {}
00581 
00582 public:
00583 
00586     particle_status_e get_status() { return( _status ); }
00587 
00590     void set_status( particle_status_e status ) { _status = status; }
00591 
00594     double IQ() const { return( _IQ ); }
00595 
00598     double qm() const { return( _qm ); }
00599 
00600 };
00601 
00602 
00611 template<class PP> class Particle : public ParticleBase
00612 {
00613     std::vector<PP>     _trajectory;   
00614     PP                  _x;            
00616 public:
00617 
00626     Particle( double IQ, double q, double m, const PP &x ) 
00627         : ParticleBase(IQ,q,m), _x(x) {}
00628 
00631     ~Particle() {}
00632 
00635     double &operator()( int i ) { return( _x(i) ); }
00636 
00639     const double &operator()( int i ) const { return( _x(i) ); }
00640 
00643     Vec3D location() const { return( _x.location() ); }
00644 
00647     Vec3D velocity() const { return( _x.velocity() ); }
00648 
00651     PP &x() { return( _x ); }
00652 
00655     const PP &x() const { return( _x ); }
00656 
00659     PP &traj( int i ) { return( _trajectory[i] ); }
00660 
00663     const PP &traj( int i ) const { return( _trajectory[i] ); }
00664 
00667     size_t traj_size( void ) const { return( _trajectory.size() ); }
00668 
00671     void add_trajectory_point( const PP &x ) { _trajectory.push_back( x ); }
00672 
00675     void copy_trajectory( const std::vector<PP> &traj ) { _trajectory = traj; }
00676 
00679     void clear_trajectory( void ) { _trajectory.clear(); }
00680 
00683     void debug_print( void ) const {
00684         size_t a;
00685         std::cout << "**Particle\n";
00686         if( _status == PARTICLE_OK )
00687             std::cout << "stat = PARTICLE_OK\n";
00688         else if( _status == PARTICLE_OUT )
00689             std::cout << "stat = PARTICLE_OUT\n";
00690         else if( _status == PARTICLE_COLL )
00691             std::cout << "stat = PARTICLE_COLL\n";
00692         else if( _status == PARTICLE_BADDEF )
00693             std::cout << "stat = PARTICLE_BADDEF\n";
00694         else if( _status == PARTICLE_TIME )
00695             std::cout << "stat = PARTICLE_TIME\n";
00696         else if( _status == PARTICLE_NSTP )
00697             std::cout << "stat = PARTICLE_NSTP\n";
00698         else
00699             std::cout << "status = Unknown\n";
00700         std::cout << "IQ   = " << _IQ << "\n";
00701         std::cout << "q/m  = " << _qm << "\n";
00702         std::cout << "x    = " << _x << "\n";
00703         std::cout << "Trajectory:\n";
00704         for( a = 0; a < _trajectory.size(); a++ )
00705             std::cout << "x[" << a << "] = " << _trajectory[a] << "\n";
00706 
00707         /*
00708         std::cout << "Trajectory:\n";
00709         for( a = 0; a < _trajectory.size(); a++ ) {
00710             std::cout << "x[" << a << "] = (";
00711             uint32_t b;
00712             const PP &tp = _trajectory[a];
00713             if( tp.size() > 0 ) {
00714                 for( b = 0; b < tp.size()-1; b++ )
00715                     std::cout << tp[b] << ", ";
00716                 std::cout << tp[b] << ")\n";
00717             } else {
00718                 std::cout << ")\n";
00719             }
00720         }
00721         */
00722     }
00723 };
00724 
00725 
00730 typedef Particle<ParticleP2D>  Particle2D;
00731 
00732 
00737 typedef Particle<ParticlePCyl> ParticleCyl;
00738 
00739 
00744 typedef Particle<ParticleP3D>  Particle3D;
00745 
00746 
00747 
00750 struct ParticleIteratorData {
00751     ScalarField          *_scharge;       
00752     const ScalarField    *_epot;          
00753     const Efield         *_efield;        
00754     const VectorField    *_bfield;        
00755     const Geometry       *_g;             
00756     double                _qm;            
00757     double                _phi_plasma;    
00761     ParticleIteratorData( ScalarField *scharge, const Efield *efield, 
00762                           const VectorField *bfield, const Geometry *g )
00763         : _scharge(scharge), _epot(NULL), _efield(efield), _bfield(bfield), 
00764           _g(g), _qm(0.0), _phi_plasma(0.0) {}
00765 };
00766 
00767 
00768 
00769 #endif
00770 
00771 
00772 
00773 
00774 
00775 

Generated on 18 Apr 2011 for IBSimu by  doxygen 1.6.1