IBSimu 1.0.4
particleiterator.hpp
Go to the documentation of this file.
00001 
00005 /* Copyright (c) 2005-2010 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 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 //#define DEBUG_PARTICLE_ITERATOR 1
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         // Convert PP to Vec3D
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         // If inside solid, bracket for collision point
00149         if( (a = _pidata._g->inside( v2 )) >= 7 ) {
00150             K = _pidata._g->bracket_surface( a, v2, v1, vc );
00151         } else {
00152             return( true ); // No collision happened.
00153         }
00154 
00155         // Convert Vec3D to PP
00156         for( a = 0; a < PP::size(); a++ )
00157             status_x[a] = x2[a] + K*(x1[a]-x2[a]);
00158 
00159         // Remove all points from _traj after time status_x[0].
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         // Save last trajectory point and update status
00168         _traj.push_back( status_x );
00169         particle.set_status( PARTICLE_COLL );
00170         _end_coll++;
00171 
00172         return( false ); // Collision happened.
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         // Check if found edge at first encounter
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         // Mirror traj back to _xi
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         // Mirror rest of the coldata
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         // Mirror calculation point
00246         x2[2*a+1] = 2.0*xmirror - x2[2*a+1];
00247         x2[2*a+2] *= -1.0;
00248         
00249         // Coordinates changed, reset integrator
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         // Check for collisions with solids and advance coordinates i.
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         // Check for collisions/mirroring with simulation boundary. Here
00360         // coordinates i are already advanced to next mesh.
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                 // Save intersection coordinates
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                 //std::cout << "Found valid root: " << K << "\n";
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         // Construct trajectory representation
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         // Solve trajectory intersections
00449         for( size_t a = 0; a < PP::dim(); a++ ) {
00450 
00451             // Mesh number of x1 (start point)
00452             int i = (int)floor( (x1[2*a+1]-_pidata._g->origo(a))/_pidata._g->h() );
00453             
00454             // Search to negative(dj = -1) and positive (dj = +1) mesh directions
00455             for( int dj = -1; dj <= 1; dj += 2 ) {
00456                 int j = i;
00457                 if( dj == +1 )
00458                     j = i+1;
00459                 int Kcount;  // Solution counter
00460                 double K[3]; // Solution array
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; // No valid roots
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                     // Save roots to coldata
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; // limit numerical inaccuracy
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         // Limit trajectory advance to double the simulation box
00573         // If limitation done, force to linear interpolation
00574         if( limit_trajectory_advance( x1, x2 ) )
00575             force_linear = true;
00576 
00577         // Make coldata
00578         if( _polyint && !force_linear )
00579             build_coldata_poly( particle, x1, x2 );
00580         else
00581             build_coldata_linear( particle, x1, x2 );
00582 
00583         // No intersections, nothing to do
00584         if( _coldata.size() == 0 ) {
00585 #ifdef DEBUG_PARTICLE_ITERATOR
00586             std::cout << "No coldata\n";
00587 #endif
00588             return( true );
00589         }
00590 
00591         // Sort intersections in increasing time order
00592         sort( _coldata.begin(), _coldata.end() );
00593 
00594         // Starting mesh index
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         // Process intersection points
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             // Advance particle in mesh, check for possible collisions and
00615             // do mirroring.
00616             handle_trajectory_advance( particle, a, i, x2 );
00617 
00618             // Update space charge for one mesh.
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             // Clear coldata and exit if particle collided.
00633             if( particle.get_status() != PARTICLE_OK ) {
00634                 _coldata.clear();
00635                 return( false );
00636             }
00637 
00638             // Update last intersection point xi.
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         // Get acceleration at x2
00669         double dxdt[5];
00670         PP::get_derivatives( x2[0], &x2[1], dxdt, (void *)&_pidata );
00671 
00672         // Calculate crossover point assuming zero acceleration in
00673         // r-direction and constant acceleration in x-direction
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         // Mirror x2 to x3
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         // Handle step with linear interpolation to avoid going to r<=0
00698         if( !handle_trajectory( particle, x2, x3, true ) )
00699             return( false ); // Particle done
00700 
00701         // Save trajectory calculation points
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         // Next step not a continuation of previous one, reset
00709         // integrator
00710         gsl_odeiv_step_reset( _step );
00711         gsl_odeiv_evolve_reset( _evolve );
00712 
00713         // Continue iteration at mirrored point
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         // Check if inside solids of outside geometry.
00736         if( _pidata._g->inside( x.location() ) )
00737             return( false );
00738 
00739         // Check if particle on simulation geometry border and directed outwards
00740         /*
00741         for( size_t a = 0; a < PP::dim(); a++ ) {
00742             if( x[2*a+1] == _pidata._g->origo(a) && x[2*a+2] < 0.0 ) {
00743                 if( _mirror[2*a] ) {
00744                     x[2*a+2] *= -1.0;
00745 #ifdef DEBUG_PARTICLE_ITERATOR
00746         std::cout << "Mirroring to:\n";
00747         std::cout << "  x = " << x << "\n";
00748 #endif
00749                 } else {
00750                     return( false );
00751                 }
00752 
00753             } else if( x[2*a+1] == _pidata._g->max(a) & x[2*a+2] > 0.0 ) {
00754                 if( _mirror[2*a+1] ) {
00755                     x[2*a+2] *= -1.0;
00756 #ifdef DEBUG_PARTICLE_ITERATOR
00757         std::cout << "Mirroring to:\n";
00758         std::cout << "  x = " << x << "\n";
00759 #endif
00760                 } else {
00761                     return( false );
00762                 }
00763             }
00764         }
00765 
00766 #ifdef DEBUG_PARTICLE_ITERATOR
00767         std::cout << "Definition ok\n";
00768 #endif
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             //std::cout << "spd += " << dxdt[2*a]*dxdt[2*a] << "\n";
00780             spd += dxdt[2*a]*dxdt[2*a];
00781             //std::cout << "acc += " << dxdt[2*a+1]*dxdt[2*a+1] << "\n";
00782             acc += dxdt[2*a+1]*dxdt[2*a+1];
00783         }
00784         if( _pidata._g->geom_mode() == MODE_CYL ) {
00785             //std::cout << "MODE_CYL\n";
00786             //std::cout << "spd += " << x[3]*x[3]*x[5]*x[5] << "\n";
00787             spd += x[3]*x[3]*x[5]*x[5];
00788             //std::cout << "acc += " << x[3]*x[3]*dxdt[4]*dxdt[4] << "\n";
00789             acc += x[3]*x[3]*dxdt[4]*dxdt[4];
00790         }
00791         //std::cout << "spd = " << sqrt(spd) << "\n";
00792         //std::cout << "acc = " << sqrt(acc) << "\n";
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         // Initialize mirroring
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         // Initialize system of ordinary differential equations (ODE)
00844         _system.jacobian  = NULL;
00845         _system.params    = (void *)&_pidata;
00846         _system.function  = PP::get_derivatives;
00847         _system.dimension = PP::size()-1; // Time is not part of differential equation dimensions
00848 
00849         // Make scale
00850         // 2D:  x vx y vy
00851         // Cyl: x vx r vr omega
00852         // 3D:  x vx y vy z vz
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         // Initialize ODE solver
00862         _step    = gsl_odeiv_step_alloc( gsl_odeiv_step_rkck, _system.dimension );
00863         //_control = gsl_odeiv_control_standard_new( _epsabs, _epsrel, 1.0, 1.0 );
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         // Copy starting point to x and 
00915         PP x = particle->x();
00916 
00917         // Check particle definition
00918         if( !check_particle_definition( x ) ) {
00919             particle->set_status( PARTICLE_BADDEF );
00920             _end_baddef++;
00921             return;
00922         }
00923         particle->x() = x;
00924 
00925         // Reset trajectory and save first trajectory point.
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         // Reset integrator
00939         gsl_odeiv_step_reset( _step );
00940         gsl_odeiv_evolve_reset( _evolve );
00941         
00942         // Make initial guess for step size
00943         //std::cout << "Guess dt ------------------------------------------------\n";
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         // Iterate ODEs until maximum steps are done, time is used 
00958         // or particle collides.
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             // Take a step.
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]; // Reset time (this shouldn't be necessary - there 
00982                                   // is a bug in GSL-1.12, report has been sent)
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             // Check step count number and step size validity
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             // Handle collisions and space charge of step.
01007             if( !handle_trajectory( *particle, x, x2 ) ) {
01008                 x = x2;
01009                 break; // Particle done
01010             }
01011 
01012             // Check if particle mirroring is required to avoid 
01013             // singularity at symmetry axis.
01014             if( axis_mirror_required( x2 ) ) {
01015                 if( !handle_axis_mirror_step( *particle, x, x2 ) )
01016                     break; // Particle done
01017             }
01018 
01019             // Propagate coordinates
01020             x = x2;
01021 
01022             // Save trajectory point
01023             _traj.push_back( x2 );
01024             
01025             // Increase step count.
01026             nstp++;
01027         }
01028 
01029 #ifdef DEBUG_PARTICLE_ITERATOR
01030         std::cout << "\n*** Done stepping ***\n";
01031 #endif
01032 
01033         // Check if step count or time limited 
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         // Save step count
01043         _sum_steps += nstp;
01044 
01045         // Save trajectory of current particle
01046         if( _trajdiv != 0 && (particle-_first) % _trajdiv == 0 )
01047             particle->copy_trajectory( _traj );
01048 
01049         // Save last particle location
01050         particle->x() = x;
01051     }
01052 
01053 /*
01054     std::cout << "Kala\n";
01055 
01056     std::cout << _coldata.capacity() << "\n";
01057     std::cout << _traj.capacity() << "\n";
01058 */
01059 };
01060 
01061 
01062 
01063 
01064 #endif
01065 
01066 
01067 
01068 
01069 
01070 
01071 
01072 
01073 
01074 
01075 
01076 
01077 
01078 
01079