particles.hpp
Go to the documentation of this file.
00001 
00005 /* Copyright (c) 2005-2011 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 <string>
00049 #include <gsl/gsl_errno.h>
00050 #include "geometry.hpp"
00051 #include "scalarfield.hpp"
00052 #include "vectorfield.hpp"
00053 #include "vec3d.hpp"
00054 #include "callback.hpp"
00055 
00056 
00057 /* Integer error value that is supposed to diffed from internal GSL error values */
00058 #define IBSIMU_DERIV_ERROR 201
00059 
00060 
00068 enum particle_status_e {
00069     PARTICLE_OK = 0,
00070     PARTICLE_OUT,
00071     PARTICLE_COLL,
00072     PARTICLE_BADDEF,
00073     PARTICLE_TIME,
00074     PARTICLE_NSTP 
00075 };
00076 
00077 
00078 
00079 /* ************************************************************************************* *
00080  * Particle point classes                                                                *
00081  * ************************************************************************************* */
00082 
00083 
00088 class ParticlePBase
00089 {
00090 
00091 };
00092 
00093 
00099 class ParticleP2D : public ParticlePBase 
00100 {
00101     double         _x[5];   
00103 public:
00104 
00107     ParticleP2D() {}
00108 
00111     ParticleP2D( double t, double x, double vx, double y, double vy ) {
00112         _x[0] = t; _x[1] = x; _x[2] = vx; _x[3] = y; _x[4] = vy;
00113     }
00114 
00117     ParticleP2D( std::istream &s ) {
00118         _x[0] = read_double( s );
00119         _x[1] = read_double( s );
00120         _x[2] = read_double( s );
00121         _x[3] = read_double( s );
00122         _x[4] = read_double( s );
00123     }
00124 
00127     static geom_mode_e geom_mode() { return(MODE_2D); }
00128 
00131     static size_t dim() { return(2); }
00132 
00135     static size_t size() { return(5); }
00136 
00148     static int get_derivatives( double t, const double *x, double *dxdt, void *data );
00149 
00154     static int trajectory_intersections_at_plane( std::vector<ParticleP2D> &intsc, 
00155                                                   int crd, double val,
00156                                                   const ParticleP2D &x1, 
00157                                                   const ParticleP2D &x2,
00158                                                   int extrapolate = 0 );
00159 
00164     static const std::string IQ_unit() { return( "A/m" ); }
00165 
00168     Vec3D location() const { return( Vec3D( _x[1], _x[3], 0.0 ) ); }
00169 
00172     Vec3D velocity() const { return( Vec3D( _x[2], _x[4], 0.0 ) ); }
00173 
00176     double speed() { return( sqrt(_x[2]*_x[2] + _x[4]*_x[4]) ); }
00177 
00180     double &operator[]( int i ) { return( _x[i] ); }
00181 
00184     const double &operator[]( int i ) const { return( _x[i] ); }
00185 
00188     double &operator()( int i ) { return( _x[i] ); }
00189 
00192     const double &operator()( int i ) const { return( _x[i] ); }
00193 
00194     ParticleP2D operator+( const ParticleP2D &pp ) const { 
00195         ParticleP2D res;
00196         res[0] = _x[0] + pp[0];
00197         res[1] = _x[1] + pp[1];
00198         res[2] = _x[2] + pp[2];
00199         res[3] = _x[3] + pp[3];
00200         res[4] = _x[4] + pp[4];
00201         return( res );
00202     }
00203 
00204     ParticleP2D operator-( const ParticleP2D &pp ) const { 
00205         ParticleP2D res;
00206         res[0] = _x[0] - pp[0];
00207         res[1] = _x[1] - pp[1];
00208         res[2] = _x[2] - pp[2];
00209         res[3] = _x[3] - pp[3];
00210         res[4] = _x[4] - pp[4];
00211         return( res );
00212     }
00213 
00214     ParticleP2D operator*( double x ) const { 
00215         ParticleP2D res;
00216         res[0] = _x[0]*x;
00217         res[1] = _x[1]*x;
00218         res[2] = _x[2]*x;
00219         res[3] = _x[3]*x;
00220         res[4] = _x[4]*x;
00221         return( res );
00222     }
00223 
00226     void save( std::ostream &s ) const {
00227         write_double( s, _x[0] );
00228         write_double( s, _x[1] );
00229         write_double( s, _x[2] );
00230         write_double( s, _x[3] );
00231         write_double( s, _x[4] );
00232     }
00233     
00234     friend ParticleP2D operator*( double x, const ParticleP2D &pp );
00235 };
00236 
00237 
00238 inline std::ostream &operator<<( std::ostream &os, const ParticleP2D &pp )
00239 {
00240     os << "("
00241        << std::setw(12) << pp(0) << ", "
00242        << std::setw(12) << pp(1) << ", "
00243        << std::setw(12) << pp(2) << ", "
00244        << std::setw(12) << pp(3) << ", "
00245        << std::setw(12) << pp(4) << ")";
00246     return( os );
00247 }
00248 
00249 
00250 inline ParticleP2D operator*( double x, const ParticleP2D &pp )
00251 {
00252     ParticleP2D res;
00253     res[0] = pp[0]*x;
00254     res[1] = pp[1]*x;
00255     res[2] = pp[2]*x;
00256     res[3] = pp[3]*x;
00257     res[4] = pp[4]*x;
00258     return( res );
00259 }
00260 
00261 
00267 class ParticlePCyl : public ParticlePBase 
00268 {
00269     double         _x[6];   
00271 public:
00272 
00275     ParticlePCyl() {}
00276 
00279     ParticlePCyl( double t, double x, double vx, double r, double vr, double w ) {
00280         _x[0] = t; _x[1] = x; _x[2] = vx; _x[3] = r; _x[4] = vr; _x[5] = w;
00281     }
00282 
00285     ParticlePCyl( std::istream &s ) {
00286         _x[0] = read_double( s );
00287         _x[1] = read_double( s );
00288         _x[2] = read_double( s );
00289         _x[3] = read_double( s );
00290         _x[4] = read_double( s );
00291         _x[5] = read_double( s );
00292     }
00293 
00296     static geom_mode_e geom_mode() { return(MODE_CYL); }
00297 
00300     static size_t dim() { return(2); }
00301 
00304     static size_t size() { return(6); }
00305 
00322     static int get_derivatives( double t, const double *x, double *dxdt, void *data );
00323 
00328     static int trajectory_intersections_at_plane( std::vector<ParticlePCyl> &intsc, 
00329                                                   int crd, double val,
00330                                                   const ParticlePCyl &x1,
00331                                                   const ParticlePCyl &x2,
00332                                                   int extrapolate = 0 );
00333 
00338     static const std::string IQ_unit() { return( "A" ); }
00339 
00342     Vec3D location() const { return( Vec3D( _x[1], _x[3], 0.0 ) ); }
00343 
00346     Vec3D velocity() const { return( Vec3D( _x[2], _x[4], _x[5]*_x[3] ) ); }
00347 
00350     double speed() { return( sqrt(_x[2]*_x[2] + _x[4]*_x[4] + _x[3]*_x[3]*_x[5]*_x[5]) ); }
00351 
00354     double &operator[]( int i ) { return( _x[i] ); }
00355 
00358     const double &operator[]( int i ) const { return( _x[i] ); }
00359 
00362     double &operator()( int i ) { return( _x[i] ); }
00363 
00366     const double &operator()( int i ) const { return( _x[i] ); }
00367 
00368     ParticlePCyl operator+( const ParticlePCyl &pp ) const { 
00369         ParticlePCyl res;
00370         res[0] = _x[0] + pp[0];
00371         res[1] = _x[1] + pp[1];
00372         res[2] = _x[2] + pp[2];
00373         res[3] = _x[3] + pp[3];
00374         res[4] = _x[4] + pp[4];
00375         res[5] = _x[5] + pp[5];
00376         return( res );
00377     }
00378 
00379     ParticlePCyl operator-( const ParticlePCyl &pp ) const { 
00380         ParticlePCyl res;
00381         res[0] = _x[0] - pp[0];
00382         res[1] = _x[1] - pp[1];
00383         res[2] = _x[2] - pp[2];
00384         res[3] = _x[3] - pp[3];
00385         res[4] = _x[4] - pp[4];
00386         res[5] = _x[5] - pp[5];
00387         return( res );
00388     }
00389 
00390     ParticlePCyl operator*( double x ) const { 
00391         ParticlePCyl res;
00392         res[0] = _x[0]*x;
00393         res[1] = _x[1]*x;
00394         res[2] = _x[2]*x;
00395         res[3] = _x[3]*x;
00396         res[4] = _x[4]*x;
00397         res[5] = _x[5]*x;
00398         return( res );
00399     }
00400 
00403     void save( std::ostream &s ) const {
00404         write_double( s, _x[0] );
00405         write_double( s, _x[1] );
00406         write_double( s, _x[2] );
00407         write_double( s, _x[3] );
00408         write_double( s, _x[4] );
00409         write_double( s, _x[5] );
00410     }
00411     
00412     friend ParticlePCyl operator*( double x, const ParticlePCyl &pp );
00413 };
00414 
00415 
00416 inline std::ostream &operator<<( std::ostream &os, const ParticlePCyl &pp )
00417 {
00418     os << "("
00419        << std::setw(12) << pp(0) << ", "
00420        << std::setw(12) << pp(1) << ", "
00421        << std::setw(12) << pp(2) << ", "
00422        << std::setw(12) << pp(3) << ", "
00423        << std::setw(12) << pp(4) << ", "
00424        << std::setw(12) << pp(5) << ")";
00425     return( os );
00426 }
00427 
00428 
00429 inline ParticlePCyl operator*( double x, const ParticlePCyl &pp )
00430 {
00431     ParticlePCyl res;
00432     res[0] = pp[0]*x;
00433     res[1] = pp[1]*x;
00434     res[2] = pp[2]*x;
00435     res[3] = pp[3]*x;
00436     res[4] = pp[4]*x;
00437     res[5] = pp[5]*x;
00438     return( res );
00439 }
00440 
00441 
00447 class ParticleP3D : public ParticlePBase 
00448 {
00449     double         _x[7];   
00451 public:
00452 
00455     ParticleP3D() {}
00456 
00459     ParticleP3D( double t, double x, double vx, double y, double vy, double z, double vz ) {
00460         _x[0] = t; _x[1] = x; _x[2] = vx; _x[3] = y; _x[4] = vy; _x[5] = z; _x[6] = vz;
00461     }
00462 
00465     ParticleP3D( std::istream &s ) {
00466         _x[0] = read_double( s );
00467         _x[1] = read_double( s );
00468         _x[2] = read_double( s );
00469         _x[3] = read_double( s );
00470         _x[4] = read_double( s );
00471         _x[5] = read_double( s );
00472         _x[6] = read_double( s );
00473     }
00474 
00477     static geom_mode_e geom_mode() { return(MODE_3D); }
00478 
00481     static size_t dim() { return(3); }
00482 
00485     static size_t size() { return(7); }
00486 
00500     static int get_derivatives( double t, const double *x, double *dxdt, void *data );
00501 
00506     static int trajectory_intersections_at_plane( std::vector<ParticleP3D> &intsc, 
00507                                                   int crd, double val,
00508                                                   const ParticleP3D &x1, 
00509                                                   const ParticleP3D &x2, 
00510                                                   int extrapolate = 0 );
00511 
00516     static const std::string IQ_unit() { return( "A" ); }
00517 
00520     Vec3D location() const { return( Vec3D( _x[1], _x[3], _x[5] ) ); }
00521 
00524     Vec3D velocity() const { return( Vec3D( _x[2], _x[4], _x[6] ) ); }
00525 
00528     double speed() { return( sqrt(_x[2]*_x[2] + _x[4]*_x[4] + _x[6]*_x[6]) ); }
00529 
00532     double &operator[]( int i ) { return( _x[i] ); }
00533 
00536     const double &operator[]( int i ) const { return( _x[i] ); }
00537 
00540     double &operator()( int i ) { return( _x[i] ); }
00541 
00544     const double &operator()( int i ) const { return( _x[i] ); }
00545 
00546     ParticleP3D operator+( const ParticleP3D &pp ) const { 
00547         ParticleP3D res;
00548         res[0] = _x[0] + pp[0];
00549         res[1] = _x[1] + pp[1];
00550         res[2] = _x[2] + pp[2];
00551         res[3] = _x[3] + pp[3];
00552         res[4] = _x[4] + pp[4];
00553         res[5] = _x[5] + pp[5];
00554         res[6] = _x[6] + pp[6];
00555         return( res );
00556     }
00557 
00558     ParticleP3D operator-( const ParticleP3D &pp ) const { 
00559         ParticleP3D res;
00560         res[0] = _x[0] - pp[0];
00561         res[1] = _x[1] - pp[1];
00562         res[2] = _x[2] - pp[2];
00563         res[3] = _x[3] - pp[3];
00564         res[4] = _x[4] - pp[4];
00565         res[5] = _x[5] - pp[5];
00566         res[6] = _x[6] - pp[6];
00567         return( res );
00568     }
00569 
00570     ParticleP3D operator*( double x ) const { 
00571         ParticleP3D res;
00572         res[0] = _x[0]*x;
00573         res[1] = _x[1]*x;
00574         res[2] = _x[2]*x;
00575         res[3] = _x[3]*x;
00576         res[4] = _x[4]*x;
00577         res[5] = _x[5]*x;
00578         res[6] = _x[6]*x;
00579         return( res );
00580     }
00581 
00584     void save( std::ostream &s ) const {
00585         write_double( s, _x[0] );
00586         write_double( s, _x[1] );
00587         write_double( s, _x[2] );
00588         write_double( s, _x[3] );
00589         write_double( s, _x[4] );
00590         write_double( s, _x[5] );
00591         write_double( s, _x[6] );
00592     }
00593     
00594     friend ParticleP3D operator*( double x, const ParticleP3D &pp );
00595 };
00596 
00597 
00598 inline std::ostream &operator<<( std::ostream &os, const ParticleP3D &pp )
00599 {
00600     os << "("
00601        << std::setw(12) << pp(0) << ", "
00602        << std::setw(12) << pp(1) << ", "
00603        << std::setw(12) << pp(2) << ", "
00604        << std::setw(12) << pp(3) << ", "
00605        << std::setw(12) << pp(4) << ", "
00606        << std::setw(12) << pp(5) << ", "
00607        << std::setw(12) << pp(6) << ")";
00608     return( os );
00609 }
00610 
00611 
00612 inline ParticleP3D operator*( double x, const ParticleP3D &pp )
00613 {
00614     ParticleP3D res;
00615     res[0] = pp[0]*x;
00616     res[1] = pp[1]*x;
00617     res[2] = pp[2]*x;
00618     res[3] = pp[3]*x;
00619     res[4] = pp[4]*x;
00620     res[5] = pp[5]*x;
00621     res[6] = pp[6]*x;
00622     return( res );
00623 }
00624 
00625 
00626 
00627 
00628 /* ************************************************************************************** *
00629  * Particle classes                                                                       *
00630  * ************************************************************************************** */
00631 
00632 
00637 class ParticleBase
00638 {
00639 protected:
00640 
00641     particle_status_e   _status;       
00642     double              _IQ;           
00653     double              _q;            
00654     double              _m;            
00656     ParticleBase( double IQ, double q, double m ) 
00657         : _status(PARTICLE_OK), _q(q) {
00658         _m = fabs(m);
00659         if( _q < 0 )
00660             _IQ = -fabs(IQ);
00661         else
00662             _IQ = fabs(IQ);
00663     }
00664 
00667     ParticleBase( std::istream &s ) {
00668         _status = (particle_status_e)read_int32( s );
00669         _IQ = read_double( s );
00670         _q = read_double( s );
00671         _m = read_double( s );  
00672     }
00673 
00674     ~ParticleBase() {}
00675 
00676 public:
00677 
00680     particle_status_e get_status() { return( _status ); }
00681 
00684     void set_status( particle_status_e status ) { _status = status; }
00685 
00691     double IQ() const { return( _IQ ); }
00692 
00695     double q() const { return( _q ); }
00696 
00699     double m() const { return( _m ); }
00700 
00703     double qm() const { return( _q/_m ); }
00704 
00707     void save( std::ostream &s ) const {
00708         write_int32( s, _status );
00709         write_double( s, _IQ );
00710         write_double( s, _q );
00711         write_double( s, _m );
00712     }
00713 };
00714 
00715 
00724 template<class PP> class Particle : public ParticleBase
00725 {
00726     std::vector<PP>     _trajectory;   
00727     PP                  _x;            
00729 public:
00730 
00739     Particle( double IQ, double q, double m, const PP &x ) 
00740         : ParticleBase(IQ,q,m), _x(x) {}
00741 
00744     Particle( std::istream &s ) 
00745         : ParticleBase( s ) {
00746 
00747         uint32_t N  = read_int32( s );
00748         _trajectory.reserve( N );
00749         for( uint32_t a = 0; a < N; a++ )
00750             _trajectory.push_back( PP( s ) );
00751         _x = PP( s );
00752     }
00753 
00756     ~Particle() {}
00757 
00760     double &operator()( int i ) { return( _x(i) ); }
00761 
00764     const double &operator()( int i ) const { return( _x(i) ); }
00765 
00768     double &operator[]( int i ) { return( _x(i) ); }
00769 
00772     const double &operator[]( int i ) const { return( _x(i) ); }
00773 
00776     Vec3D location() const { return( _x.location() ); }
00777 
00780     Vec3D velocity() const { return( _x.velocity() ); }
00781 
00784     PP &x() { return( _x ); }
00785 
00788     const PP &x() const { return( _x ); }
00789 
00792     PP &traj( int i ) { return( _trajectory[i] ); }
00793 
00796     const PP &traj( int i ) const { return( _trajectory[i] ); }
00797 
00800     size_t traj_size( void ) const { return( _trajectory.size() ); }
00801 
00804     void add_trajectory_point( const PP &x ) { _trajectory.push_back( x ); }
00805 
00808     void copy_trajectory( const std::vector<PP> &traj ) { _trajectory = traj; }
00809 
00812     void clear_trajectory( void ) { _trajectory.clear(); }
00813 
00816     void save( std::ostream &s ) const {
00817         ParticleBase::save( s );
00818         write_int32( s, _trajectory.size() );
00819         for( uint32_t a = 0; a < _trajectory.size(); a++ )
00820             _trajectory[a].save( s );
00821         _x.save( s );
00822     }
00823     
00826     void debug_print( std::ostream &os ) const {
00827         size_t a;
00828         os << "**Particle\n";
00829         switch( _status ) {
00830         case PARTICLE_OK:
00831             os << "stat = PARTICLE_OK\n";
00832             break;
00833         case PARTICLE_OUT:
00834             os << "stat = PARTICLE_OUT\n";
00835             break;
00836         case PARTICLE_COLL:
00837             os << "stat = PARTICLE_COLL\n";
00838             break;
00839         case PARTICLE_BADDEF:
00840             os << "stat = PARTICLE_BADDEF\n";
00841             break;
00842         case PARTICLE_TIME:
00843             os << "stat = PARTICLE_TIME\n";
00844             break;
00845         case PARTICLE_NSTP:
00846             os << "stat = PARTICLE_NSTP\n";
00847             break;
00848         }
00849         os << "IQ   = " << _IQ << "\n";
00850         os << "q    = " << _q << "\n";
00851         os << "m    = " << _m << "\n";
00852         os << "x    = " << _x << "\n";
00853         os << "Trajectory:\n";
00854         for( a = 0; a < _trajectory.size(); a++ )
00855             os << "x[" << a << "] = " << _trajectory[a] << "\n";
00856 
00857         /*
00858         std::cout << "Trajectory:\n";
00859         for( a = 0; a < _trajectory.size(); a++ ) {
00860             std::cout << "x[" << a << "] = (";
00861             uint32_t b;
00862             const PP &tp = _trajectory[a];
00863             if( tp.size() > 0 ) {
00864                 for( b = 0; b < tp.size()-1; b++ )
00865                     std::cout << tp[b] << ", ";
00866                 std::cout << tp[b] << ")\n";
00867             } else {
00868                 std::cout << ")\n";
00869             }
00870         }
00871         */
00872     }
00873 };
00874 
00875 
00880 typedef Particle<ParticleP2D>  Particle2D;
00881 
00882 
00887 typedef Particle<ParticlePCyl> ParticleCyl;
00888 
00889 
00894 typedef Particle<ParticleP3D>  Particle3D;
00895 
00896 
00897 
00900 struct ParticleIteratorData {
00901     ScalarField              *_scharge;  
00902     const VectorField        *_efield;   
00903     const VectorField        *_bfield;   
00904     const Geometry           *_geom;     
00905     double                    _qm;       
00906     const CallbackFunctorD_V *_bsup_cb;  
00908     ParticleIteratorData( ScalarField *scharge, const VectorField *efield, 
00909                           const VectorField *bfield, const Geometry *geom ) 
00910         : _scharge(scharge), _efield(efield), _bfield(bfield), 
00911           _geom(geom), _qm(0.0), _bsup_cb(0) {}
00912 
00915     void set_bfield_suppression_callback( const CallbackFunctorD_V *bsup_cb ) {
00916         _bsup_cb = bsup_cb;
00917     }
00918 
00919 };
00920 
00921 
00922 
00923 #endif
00924