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 <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
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
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
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
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
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