00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
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
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
00078 #define MASS_U 1.66053873e-27
00079
00080 #define CHARGE_E 1.602176462e-19
00081
00082
00083
00084
00085
00086
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
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
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
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