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 PARTICLEITERATOR_HPP
00044 #define PARTICLEITERATOR_HPP 1
00045
00046
00047 #include <vector>
00048 #include <iostream>
00049 #include <algorithm>
00050 #include <iomanip>
00051 #include <gsl/gsl_odeiv.h>
00052 #include <gsl/gsl_poly.h>
00053 #include "geometry.hpp"
00054 #include "particles.hpp"
00055 #include "efield.hpp"
00056 #include "scalarfield.hpp"
00057 #include "scharge.hpp"
00058 #include "scheduler.hpp"
00059 #include "polysolver.hpp"
00060
00061
00062
00063
00064
00067 enum particle_iterator_type_e {
00068 PARTICLE_ITERATOR_ADAPTIVE = 0,
00069 PARTICLE_ITERATOR_FIXED_STEP_LEN
00070 };
00071
00072
00078 template <class PP> class ParticleIterator {
00079
00082 struct ColData {
00083 PP _x;
00084 int _dir;
00087 ColData( PP x, int dir ) : _x(x), _dir(dir) {}
00088
00089 bool operator<( const ColData &cd ) const {
00090 return( _x[0] < cd._x[0] );
00091 }
00092 };
00093
00094 gsl_odeiv_system _system;
00095 gsl_odeiv_step *_step;
00096 gsl_odeiv_control *_control;
00097 gsl_odeiv_evolve *_evolve;
00099 particle_iterator_type_e _type;
00101 bool _polyint;
00102 double _epsabs;
00103 double _epsrel;
00104 uint32_t _maxsteps;
00105 double _maxt;
00106 uint32_t _trajdiv;
00108 bool _mirror[6];
00110 Particle<PP> *_first;
00111 ParticleIteratorData _pidata;
00113 PP _xi;
00115 std::vector<PP> _traj;
00116 std::vector<ColData> _coldata;
00118 uint32_t _end_time;
00119 uint32_t _end_step;
00120 uint32_t _end_out;
00122 uint32_t _end_coll;
00124 uint32_t _end_baddef;
00125 uint32_t _sum_steps;
00136 bool check_collision( Particle<PP> &particle, const PP &x1, const PP &x2, PP &status_x ) {
00137
00138 size_t a;
00139 double K;
00140 Vec3D v1, v2, vc;
00141
00142
00143 for( a = 0; a < (PP::size()-1)/2; a++ ) {
00144 v1[a] = x1[2*a+1];
00145 v2[a] = x2[2*a+1];
00146 }
00147
00148
00149 if( (a = _pidata._g->inside( v2 )) >= 7 ) {
00150 K = _pidata._g->bracket_surface( a, v2, v1, vc );
00151 } else {
00152 return( true );
00153 }
00154
00155
00156 for( a = 0; a < PP::size(); a++ )
00157 status_x[a] = x2[a] + K*(x1[a]-x2[a]);
00158
00159
00160 for( a = _traj.size()-1; a > 0; a-- ) {
00161 if( _traj[a][0] > status_x[0] )
00162 _traj.pop_back();
00163 else
00164 break;
00165 }
00166
00167
00168 _traj.push_back( status_x );
00169 particle.set_status( PARTICLE_COLL );
00170 _end_coll++;
00171
00172 return( false );
00173 }
00174
00175
00183 void handle_mirror( size_t c, int i[3], size_t a, int border, PP &x2 ) {
00184
00185 #ifdef DEBUG_PARTICLE_ITERATOR
00186 std::cout << " handle_mirror( c = " << c
00187 << ", i = (" << i[0] << ", " << i[1] << ", " << i[2]
00188 << "), a = " << a << ", border = " << border
00189 << ")\n";
00190 #endif
00191
00192 double xmirror;
00193 if( border < 0 ) {
00194 xmirror = _pidata._g->origo(a);
00195 i[a] = -i[a]-1;
00196 } else {
00197 xmirror = _pidata._g->max(a);
00198 i[a] = 2*_pidata._g->size(a)-i[a]-3;
00199 }
00200
00201 #ifdef DEBUG_PARTICLE_ITERATOR
00202 std::cout << " xmirror = " << xmirror << "\n";
00203 std::cout << " i = (" << i[0] << ", " << i[1] << ", " << i[2] << ")\n";
00204 std::cout << " xi = " << _xi << "\n";
00205 #endif
00206
00207
00208 bool caught_at_boundary = false;
00209 if( _coldata[c]._dir == border*((int)a+1) &&
00210 ( i[a] == 0 || i[a] == (int)_pidata._g->size(a)-2 ) ) {
00211 caught_at_boundary = true;
00212 #ifdef DEBUG_PARTICLE_ITERATOR
00213 std::cout << " caught_at_boundary\n";
00214 #endif
00215 }
00216
00217
00218 if( caught_at_boundary ) {
00219 _traj.push_back( _coldata[c]._x );
00220 } else {
00221 for( int b = _traj.size()-1; b > 0; b-- ) {
00222 if( _traj[b][0] >= _xi[0] ) {
00223
00224 #ifdef DEBUG_PARTICLE_ITERATOR
00225 std::cout << " mirroring traj[" << b << "] = " << _traj[b] << "\n";
00226 #endif
00227 _traj[b][2*a+1] = 2.0*xmirror - _traj[b][2*a+1];
00228 _traj[b][2*a+2] *= -1.0;
00229 } else
00230 break;
00231 }
00232 }
00233
00234
00235 for( size_t b = c; b < _coldata.size(); b++ ) {
00236 if( (size_t)abs(_coldata[b]._dir) == a+1 )
00237 _coldata[b]._dir *= -1;
00238 _coldata[b]._x[2*a+1] = 2.0*xmirror - _coldata[b]._x[2*a+1];
00239 _coldata[b]._x[2*a+2] *= -1.0;
00240 }
00241
00242 if( caught_at_boundary )
00243 _traj.push_back( _coldata[c]._x );
00244
00245
00246 x2[2*a+1] = 2.0*xmirror - x2[2*a+1];
00247 x2[2*a+2] *= -1.0;
00248
00249
00250 gsl_odeiv_step_reset( _step );
00251 gsl_odeiv_evolve_reset( _evolve );
00252 }
00253
00254
00255 void handle_collision( Particle<PP> &particle, size_t c, PP &status_x ) {
00256
00257 #ifdef DEBUG_PARTICLE_ITERATOR
00258 std::cout << " handle_collision()\n";
00259 #endif
00260
00261 _traj.push_back( _coldata[c]._x );
00262 status_x = _coldata[c]._x;
00263 particle.set_status( PARTICLE_OUT );
00264 _end_out++;
00265 }
00266
00267
00276 bool handle_trajectory_advance( Particle<PP> &particle, size_t c, int i[3], PP &x2 ) {
00277
00278
00279 if( PP::dim() == 2 ) {
00280 if( _coldata[c]._dir == -1 ) {
00281 if( ( abs(_pidata._g->mesh(i[0], i[1] )) >= 7 ||
00282 abs(_pidata._g->mesh(i[0], i[1]+1)) >= 7 ) &&
00283 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00284 return( false );
00285 i[0]--;
00286 } else if( _coldata[c]._dir == +1 ) {
00287 if( ( abs(_pidata._g->mesh(i[0]+1,i[1] )) >= 7 ||
00288 abs(_pidata._g->mesh(i[0]+1,i[1]+1)) >= 7 ) &&
00289 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00290 return( false );
00291 i[0]++;
00292 } else if( _coldata[c]._dir == -2 ) {
00293 if( ( abs(_pidata._g->mesh(i[0], i[1] )) >= 7 ||
00294 abs(_pidata._g->mesh(i[0]+1,i[1] )) >= 7 ) &&
00295 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00296 return( false );
00297 i[1]--;
00298 } else {
00299 if( ( abs(_pidata._g->mesh(i[0], i[1]+1)) >= 7 ||
00300 abs(_pidata._g->mesh(i[0]+1,i[1]+1)) >= 7 ) &&
00301 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00302 return( false );
00303 i[1]++;
00304 }
00305 } else if( PP::dim() == 3 ) {
00306 if( _coldata[c]._dir == -1 ) {
00307 if( ( abs(_pidata._g->mesh(i[0], i[1], i[2] )) >= 7 ||
00308 abs(_pidata._g->mesh(i[0], i[1]+1,i[2] )) >= 7 ||
00309 abs(_pidata._g->mesh(i[0], i[1], i[2]+1)) >= 7 ||
00310 abs(_pidata._g->mesh(i[0], i[1]+1,i[2]+1)) >= 7 ) &&
00311 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00312 return( false );
00313 i[0]--;
00314 } else if( _coldata[c]._dir == +1 ) {
00315 if( ( abs(_pidata._g->mesh(i[0]+1,i[1], i[2] )) >= 7 ||
00316 abs(_pidata._g->mesh(i[0]+1,i[1]+1,i[2] )) >= 7 ||
00317 abs(_pidata._g->mesh(i[0]+1,i[1], i[2]+1)) >= 7 ||
00318 abs(_pidata._g->mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) &&
00319 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00320 return( false );
00321 i[0]++;
00322 } else if( _coldata[c]._dir == -2 ) {
00323 if( ( abs(_pidata._g->mesh(i[0], i[1],i[2] )) >= 7 ||
00324 abs(_pidata._g->mesh(i[0]+1,i[1],i[2] )) >= 7 ||
00325 abs(_pidata._g->mesh(i[0], i[1],i[2]+1)) >= 7 ||
00326 abs(_pidata._g->mesh(i[0]+1,i[1],i[2]+1)) >= 7 ) &&
00327 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00328 return( false );
00329 i[1]--;
00330 } else if( _coldata[c]._dir == +2 ) {
00331 if( ( abs(_pidata._g->mesh(i[0], i[1]+1,i[2] )) >= 7 ||
00332 abs(_pidata._g->mesh(i[0]+1,i[1]+1,i[2] )) >= 7 ||
00333 abs(_pidata._g->mesh(i[0], i[1]+1,i[2]+1)) >= 7 ||
00334 abs(_pidata._g->mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) &&
00335 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00336 return( false );
00337 i[1]++;
00338 } else if( _coldata[c]._dir == -3 ) {
00339 if( ( abs(_pidata._g->mesh(i[0], i[1], i[2] )) >= 7 ||
00340 abs(_pidata._g->mesh(i[0]+1,i[1], i[2] )) >= 7 ||
00341 abs(_pidata._g->mesh(i[0], i[1]+1,i[2])) >= 7 ||
00342 abs(_pidata._g->mesh(i[0]+1,i[1]+1,i[2])) >= 7 ) &&
00343 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00344 return( false );
00345 i[2]--;
00346 } else {
00347 if( ( abs(_pidata._g->mesh(i[0], i[1], i[2]+1)) >= 7 ||
00348 abs(_pidata._g->mesh(i[0]+1,i[1], i[2]+1)) >= 7 ||
00349 abs(_pidata._g->mesh(i[0], i[1]+1,i[2]+1)) >= 7 ||
00350 abs(_pidata._g->mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) &&
00351 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00352 return( false );
00353 i[2]++;
00354 }
00355 } else {
00356 throw( Error( ERROR_LOCATION, "unsupported dimension number" ) );
00357 }
00358
00359
00360
00361 for( size_t a = 0; a < PP::dim(); a++ ) {
00362
00363 if( i[a] < 0 ) {
00364 if( _mirror[2*a] )
00365 handle_mirror( c, i, a, -1, x2 );
00366 else {
00367 handle_collision( particle, c, x2 );
00368 return( false );
00369 }
00370 } else if( i[a] >= (_pidata._g->size(a)-1) ) {
00371 if( _mirror[2*a+1] )
00372 handle_mirror( c, i, a, +1, x2 );
00373 else {
00374 handle_collision( particle, c, x2 );
00375 return( false );
00376 }
00377 }
00378 }
00379
00380 return( true );
00381 }
00382
00391 void build_coldata_linear( Particle<PP> &particle, const PP &x1, const PP &x2 ) {
00392
00393 int a1, a2;
00394
00395 for( size_t a = 0; a < PP::dim(); a++ ) {
00396
00397 a1 = (int)floor( (x1[2*a+1]-_pidata._g->origo(a))/_pidata._g->h() );
00398 a2 = (int)floor( (x2[2*a+1]-_pidata._g->origo(a))/_pidata._g->h() );
00399 if( a1 > a2 ) {
00400 int a = a2;
00401 a2 = a1;
00402 a1 = a;
00403 }
00404
00405 for( int b = a1+1; b <= a2; b++ ) {
00406
00407
00408 double K = (b*_pidata._g->h() + _pidata._g->origo(a) - x1[2*a+1]) /
00409 (x2[2*a+1] - x1[2*a+1]);
00410 if( K < 0.0 ) K = 0.0;
00411 else if( K > 1.0 ) K = 1.0;
00412
00413
00414 if( x2[2*a+1] > x1[2*a+1] )
00415 _coldata.push_back( ColData( x1 + (x2-x1)*K, a+1 ) );
00416 else
00417 _coldata.push_back( ColData( x1 + (x2-x1)*K, -a-1 ) );
00418 }
00419 }
00420
00421 #ifdef DEBUG_PARTICLE_ITERATOR
00422 std::cout << " Coldata linear built\n";
00423 #endif
00424 }
00425
00434 void build_coldata_poly( Particle<PP> &particle, const PP &x1, const PP &x2 ) {
00435
00436 #ifdef DEBUG_PARTICLE_ITERATOR
00437 std::cout << "Building coldata using polynomial interpolation\n";
00438 #endif
00439
00440
00441 TrajectoryRep1D traj[PP::dim()];
00442 for( size_t a = 0; a < PP::dim(); a++ ) {
00443 traj[a].construct( x2[0]-x1[0],
00444 x1[2*a+1], x1[2*a+2],
00445 x2[2*a+1], x2[2*a+2] );
00446 }
00447
00448
00449 for( size_t a = 0; a < PP::dim(); a++ ) {
00450
00451
00452 int i = (int)floor( (x1[2*a+1]-_pidata._g->origo(a))/_pidata._g->h() );
00453
00454
00455 for( int dj = -1; dj <= 1; dj += 2 ) {
00456 int j = i;
00457 if( dj == +1 )
00458 j = i+1;
00459 int Kcount;
00460 double K[3];
00461 while( 1 ) {
00462 double val = _pidata._g->origo(a) + _pidata._g->h() * j;
00463 #ifdef DEBUG_PARTICLE_ITERATOR
00464 std::cout << " Searching intersections at coord(" << a << ") = " << val << "\n";
00465 #endif
00466 Kcount = traj[a].solve( K, val );
00467 if( Kcount == 0 )
00468 break;
00469
00470 #ifdef DEBUG_PARTICLE_ITERATOR
00471 std::cout << " Found " << Kcount << " valid roots: ";
00472 for( int p = 0; p < Kcount; p++ )
00473 std::cout << K[p] << " ";
00474 std::cout << "\n";
00475 #endif
00476
00477
00478 for( int b = 0; b < Kcount; b++ ) {
00479 PP xcol;
00480 double x, v;
00481 xcol(0) = x1[0] + K[b]*(x2[0]-x1[0]);
00482 for( size_t c = 0; c < PP::dim(); c++ ) {
00483 traj[c].coord( x, v, K[b] );
00484 if( a == c )
00485 xcol[2*c+1] = val;
00486 else
00487 xcol[2*c+1] = x;
00488 xcol[2*c+2] = v;
00489 }
00490 if( _pidata._g->geom_mode() == MODE_CYL )
00491 xcol[5] = x1[5] + K[b]*(x2[5]-x1[5]);
00492 if( xcol[2*a+2] >= 0.0 )
00493 _coldata.push_back( ColData( xcol, a+1 ) );
00494 else
00495 _coldata.push_back( ColData( xcol, -a-1 ) );
00496 }
00497
00498 j += dj;
00499 }
00500 }
00501 }
00502
00503 #ifdef DEBUG_PARTICLE_ITERATOR
00504 std::cout << " Coldata built\n";
00505 #endif
00506 }
00507
00513 bool limit_trajectory_advance( const PP &x1, PP &x2 ) {
00514
00515 bool touched = false;
00516
00517 for( size_t a = 0; a < PP::dim(); a++ ) {
00518
00519 double lim1 = _pidata._g->origo(a) -
00520 (_pidata._g->size(a)-1)*_pidata._g->h();
00521 double lim2 = _pidata._g->origo(a) +
00522 2*(_pidata._g->size(a)-1)*_pidata._g->h();
00523
00524 if( x2[2*a+1] < lim1 ) {
00525
00526 double K = (lim1 - x1[2*a+1]) / (x2[2*a+1] - x1[2*a+1]);
00527 x2 = x1 + K*(x2-x1);
00528 touched = true;
00529 #ifdef DEBUG_PARTICLE_ITERATOR
00530 std::cout << "Limiting step to:\n";
00531 std::cout << " x2: " << x2 << "\n";
00532 #endif
00533 } else if(x2[2*a+1] > lim2 ) {
00534
00535 double K = (lim2 - x1[2*a+1]) / (x2[2*a+1] - x1[2*a+1]);
00536 x2 = x1 + K*(x2-x1);
00537 touched = true;
00538 #ifdef DEBUG_PARTICLE_ITERATOR
00539 std::cout << "Limiting step to:\n";
00540 std::cout << " x2: " << x2 << "\n";
00541 #endif
00542 }
00543 }
00544
00545 return( touched );
00546 }
00547
00563 bool handle_trajectory( Particle<PP> &particle, const PP &x1, PP &x2,
00564 bool force_linear=false ) {
00565
00566 #ifdef DEBUG_PARTICLE_ITERATOR
00567 std::cout << "Handle trajectory from x1 to x2:\n";
00568 std::cout << " x1: " << x1 << "\n";
00569 std::cout << " x2: " << x2 << "\n";
00570 #endif
00571
00572
00573
00574 if( limit_trajectory_advance( x1, x2 ) )
00575 force_linear = true;
00576
00577
00578 if( _polyint && !force_linear )
00579 build_coldata_poly( particle, x1, x2 );
00580 else
00581 build_coldata_linear( particle, x1, x2 );
00582
00583
00584 if( _coldata.size() == 0 ) {
00585 #ifdef DEBUG_PARTICLE_ITERATOR
00586 std::cout << "No coldata\n";
00587 #endif
00588 return( true );
00589 }
00590
00591
00592 sort( _coldata.begin(), _coldata.end() );
00593
00594
00595 int i[3] = {0, 0, 0};
00596 for( size_t cdir = 0; cdir < PP::dim(); cdir++ )
00597 i[cdir] = (int)floor( (x1[2*cdir+1]-_pidata._g->origo(cdir))/_pidata._g->h() );
00598
00599
00600 #ifdef DEBUG_PARTICLE_ITERATOR
00601 std::cout << "Process coldata points:\n";
00602 #endif
00603 for( size_t a = 0; a < _coldata.size(); a++ ) {
00604
00605 #ifdef DEBUG_PARTICLE_ITERATOR
00606 std::cout << " Coldata " << std::setw(4) << a << ": "
00607 << _coldata[a]._x << ", "
00608 << std::setw(3) << i[0] << " "
00609 << std::setw(3) << i[1] << " "
00610 << std::setw(3) << i[2] << " "
00611 << std::setw(3) << _coldata[a]._dir << "\n";
00612 #endif
00613
00614
00615
00616 handle_trajectory_advance( particle, a, i, x2 );
00617
00618
00619 if( _pidata._scharge )
00620 scharge_add_from_trajectory( *_pidata._scharge, particle.IQ(),
00621 _xi, _coldata[a]._x );
00622
00623 #ifdef DEBUG_PARTICLE_ITERATOR
00624 if( particle.get_status() == PARTICLE_OUT ) {
00625 std::cout << " Particle out\n";
00626 std::cout << " x = " << x2 << "\n";
00627 } else if( particle.get_status() == PARTICLE_COLL ) {
00628 std::cout << " Particle collided\n";
00629 std::cout << " x = " << x2 << "\n";
00630 }
00631 #endif
00632
00633 if( particle.get_status() != PARTICLE_OK ) {
00634 _coldata.clear();
00635 return( false );
00636 }
00637
00638
00639 _xi = _coldata[a]._x;
00640 }
00641
00642 #ifdef DEBUG_PARTICLE_ITERATOR
00643 std::cout << "Coldata done\n";
00644 #endif
00645 _coldata.clear();
00646 return( true );
00647 }
00648
00649
00652 bool axis_mirror_required( const PP &x2 ) {
00653 return( _pidata._g->geom_mode() == MODE_CYL &&
00654 x2[4] < 0.0 &&
00655 x2[3] <= 0.01*_pidata._g->h() &&
00656 x2[3]*fabs(x2[5]) <= 1.0e-9*fabs(x2[4]) );
00657
00658 }
00659
00660
00666 bool handle_axis_mirror_step( Particle<PP> &particle, const PP &x1, PP &x2 ) {
00667
00668
00669 double dxdt[5];
00670 PP::get_derivatives( x2[0], &x2[1], dxdt, (void *)&_pidata );
00671
00672
00673
00674 double dt = -x2[3]/x2[4];
00675 PP xc;
00676 xc[0] = x2[0]+dt;
00677 xc[1] = x2[1]+(x2[2]+0.5*dxdt[1]*dt)*dt;
00678 xc[2] = x2[2];
00679 xc[3] = x2[3]+x2[4]*dt;
00680 xc[4] = x2[4];
00681 xc[5] = x2[5];
00682
00683
00684 PP x3 = 2*xc - x2;
00685 x3[3] *= -1.0;
00686 x3[4] *= -1.0;
00687 x3[5] *= -1.0;
00688
00689 #ifdef DEBUG_PARTICLE_ITERATOR
00690 std::cout << "Particle mirror:\n";
00691 std::cout << " x1: " << x1 << "\n";
00692 std::cout << " x2: " << x2 << "\n";
00693 std::cout << " xc: " << xc << "\n";
00694 std::cout << " x3: " << x3 << "\n";
00695 #endif
00696
00697
00698 if( !handle_trajectory( particle, x2, x3, true ) )
00699 return( false );
00700
00701
00702 _traj.push_back( x2 );
00703 _traj.push_back( xc );
00704 xc[4] *= -1.0;
00705 xc[5] *= -1.0;
00706 _traj.push_back( xc );
00707
00708
00709
00710 gsl_odeiv_step_reset( _step );
00711 gsl_odeiv_evolve_reset( _evolve );
00712
00713
00714 x2 = x3;
00715 return( true );
00716 }
00717
00728 bool check_particle_definition( PP &x ) {
00729
00730 #ifdef DEBUG_PARTICLE_ITERATOR
00731 std::cout << "Particle defined at:\n";
00732 std::cout << " x = " << x << "\n";
00733 #endif
00734
00735
00736 if( _pidata._g->inside( x.location() ) )
00737 return( false );
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771 return( true );
00772 }
00773
00774 double calculate_dt( const PP &x, const double *dxdt ) {
00775
00776 double spd = 0.0, acc = 0.0;
00777
00778 for( size_t a = 0; a < (PP::size()-1)/2; a++ ) {
00779
00780 spd += dxdt[2*a]*dxdt[2*a];
00781
00782 acc += dxdt[2*a+1]*dxdt[2*a+1];
00783 }
00784 if( _pidata._g->geom_mode() == MODE_CYL ) {
00785
00786
00787 spd += x[3]*x[3]*x[5]*x[5];
00788
00789 acc += x[3]*x[3]*dxdt[4]*dxdt[4];
00790 }
00791
00792
00793 spd = _pidata._g->h() / sqrt(spd);
00794 acc = sqrt( 2.0*_pidata._g->h() / sqrt(acc) );
00795
00796 return( spd < acc ? spd : acc );
00797 }
00798
00799 public:
00800
00827 ParticleIterator( particle_iterator_type_e type, double epsabs, double epsrel,
00828 bool polyint, uint32_t maxsteps, double maxt,
00829 uint32_t trajdiv, bool mirror[6], ScalarField *scharge, const Efield *efield,
00830 const VectorField *bfield, const Geometry *g, Particle<PP> *first )
00831 : _type(type), _polyint(polyint), _epsabs(epsabs), _epsrel(epsrel), _maxsteps(maxsteps), _maxt(maxt),
00832 _trajdiv(trajdiv), _first(first), _pidata(scharge,efield,bfield,g), _end_time(0),
00833 _end_step(0), _end_out(0), _end_coll(0), _end_baddef(0), _sum_steps(0) {
00834
00835
00836 _mirror[0] = mirror[0];
00837 _mirror[1] = mirror[1];
00838 _mirror[2] = mirror[2];
00839 _mirror[3] = mirror[3];
00840 _mirror[4] = mirror[4];
00841 _mirror[5] = mirror[5];
00842
00843
00844 _system.jacobian = NULL;
00845 _system.params = (void *)&_pidata;
00846 _system.function = PP::get_derivatives;
00847 _system.dimension = PP::size()-1;
00848
00849
00850
00851
00852
00853 double scale_abs[PP::size()-1];
00854 for( uint32_t a = 0; a < (uint32_t)PP::size()-2; a+=2 ) {
00855 scale_abs[a+0] = 1.0;
00856 scale_abs[a+1] = 1.0e6;
00857 }
00858 if( _pidata._g->geom_mode() == MODE_CYL )
00859 scale_abs[4] = 1.0;
00860
00861
00862 _step = gsl_odeiv_step_alloc( gsl_odeiv_step_rkck, _system.dimension );
00863
00864 _control = gsl_odeiv_control_scaled_new( _epsabs, _epsrel, 1.0, 1.0, scale_abs, PP::size()-1 );
00865 _evolve = gsl_odeiv_evolve_alloc( _system.dimension );
00866 }
00867
00868
00871 ~ParticleIterator() {
00872 gsl_odeiv_evolve_free( _evolve );
00873 gsl_odeiv_control_free( _control );
00874 gsl_odeiv_step_free( _step );
00875 }
00876
00877
00880 void enable_nsimp_plasma_threshold( const ScalarField *epot, double phi_plasma ) {
00881 _pidata._epot = epot;
00882 _pidata._phi_plasma = phi_plasma;
00883 }
00884
00885
00893 void get_statistics( uint32_t &end_time, uint32_t &end_step, uint32_t &end_out,
00894 uint32_t &end_coll, uint32_t &end_baddef, uint32_t &sum_steps ) {
00895 end_time = _end_time;
00896 end_step = _end_step;
00897 end_out = _end_out;
00898 end_coll = _end_coll;
00899 end_baddef = _end_baddef;
00900 sum_steps = _sum_steps;
00901 }
00902
00911 void operator()( Particle<PP> *particle,
00912 Scheduler<ParticleIterator<PP>,Particle<PP>,Error> &scheduler ) {
00913
00914
00915 PP x = particle->x();
00916
00917
00918 if( !check_particle_definition( x ) ) {
00919 particle->set_status( PARTICLE_BADDEF );
00920 _end_baddef++;
00921 return;
00922 }
00923 particle->x() = x;
00924
00925
00926 _traj.clear();
00927 _traj.push_back( x );
00928 #ifdef DEBUG_PARTICLE_ITERATOR
00929 std::cout << x[0] << " "
00930 << x[1] << " "
00931 << x[2] << " "
00932 << x[3] << " "
00933 << x[4] << "\n";
00934 #endif
00935 _pidata._qm = particle->qm();
00936 _xi = x;
00937
00938
00939 gsl_odeiv_step_reset( _step );
00940 gsl_odeiv_evolve_reset( _evolve );
00941
00942
00943
00944 double dxdt[PP::size()-1];
00945 PP::get_derivatives( 0.0, &x[1], dxdt, (void *)&_pidata );
00946 double dt = calculate_dt( x, dxdt );
00947
00948 #ifdef DEBUG_PARTICLE_ITERATOR
00949 std::cout << "dxdt = ";
00950 for( size_t a = 0; a < PP::size()-1; a++ )
00951 std::cout << dxdt[a] << " ";
00952 std::cout << "\n";
00953 std::cout << "dt = " << dt << "\n";
00954 std::cout << "*** Starting iteration\n";
00955 #endif
00956
00957
00958
00959 PP x2;
00960 size_t nstp = 0;
00961 while( nstp < _maxsteps && x[0] < _maxt ) {
00962
00963 #ifdef DEBUG_PARTICLE_ITERATOR
00964 std::cout << "\n*** Step ***\n";
00965 std::cout << " x = " << x2 << "\n";
00966 std::cout << " dt = " << dt << " (proposed)\n";
00967 #endif
00968
00969
00970 x2 = x;
00971
00972 while( nstp < _maxsteps ) {
00973 int retval = gsl_odeiv_evolve_apply( _evolve, _control, _step, &_system,
00974 &x2[0], _maxt, &dt, &x2[1] );
00975 if( retval == IBSIMU_DERIV_ERROR ) {
00976 #ifdef DEBUG_PARTICLE_ITERATOR
00977 std::cout << "Step rejected\n";
00978 std::cout << " x2 = " << x2 << "\n";
00979 std::cout << " dt = " << dt << "\n";
00980 #endif
00981 x2[0] = x[0];
00982
00983 dt *= 0.5;
00984 nstp++;
00985 continue;
00986 } else if( retval == GSL_SUCCESS ) {
00987 break;
00988 } else {
00989 throw( Error( ERROR_LOCATION, "gsl_odeiv_evolve_apply failed" ) );
00990 }
00991 }
00992
00993
00994 if( nstp >= _maxsteps )
00995 break;
00996 if( x2[0] == x[0] )
00997 throw( Error( ERROR_LOCATION, "too small step size" ) );
00998
00999 #ifdef DEBUG_PARTICLE_ITERATOR
01000 std::cout << "Step accepted from x1 to x2:\n";
01001 std::cout << " dt = " << dt << " (taken)\n";
01002 std::cout << " x1 = " << x << "\n";
01003 std::cout << " x2 = " << x2 << "\n";
01004 #endif
01005
01006
01007 if( !handle_trajectory( *particle, x, x2 ) ) {
01008 x = x2;
01009 break;
01010 }
01011
01012
01013
01014 if( axis_mirror_required( x2 ) ) {
01015 if( !handle_axis_mirror_step( *particle, x, x2 ) )
01016 break;
01017 }
01018
01019
01020 x = x2;
01021
01022
01023 _traj.push_back( x2 );
01024
01025
01026 nstp++;
01027 }
01028
01029 #ifdef DEBUG_PARTICLE_ITERATOR
01030 std::cout << "\n*** Done stepping ***\n";
01031 #endif
01032
01033
01034 if( nstp == _maxsteps ) {
01035 particle->set_status( PARTICLE_NSTP );
01036 _end_step++;
01037 } else if( x[0] >= _maxt ) {
01038 particle->set_status( PARTICLE_TIME );
01039 _end_time++;
01040 }
01041
01042
01043 _sum_steps += nstp;
01044
01045
01046 if( _trajdiv != 0 && (particle-_first) % _trajdiv == 0 )
01047 particle->copy_trajectory( _traj );
01048
01049
01050 particle->x() = x;
01051 }
01052
01053
01054
01055
01056
01057
01058
01059 };
01060
01061
01062
01063
01064 #endif
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079