IBSimu 1.0.4
|
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