ParallelTTracker.cpp 126 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
// ------------------------------------------------------------------------
// $RCSfile: ParallelTTracker.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.2.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: ParallelTTracker
//   The visitor class for tracking particles with time as independent
//   variable.
//
// ------------------------------------------------------------------------
//
// $Date: 2004/11/12 20:10:11 $
// $Author: adelmann $
//
// ------------------------------------------------------------------------

kraus's avatar
kraus committed
20 21
#include "Algorithms/ParallelTTracker.h"

gsell's avatar
gsell committed
22 23 24 25 26 27
#include <cfloat>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <string>
28
#include <limits>
adelmann's avatar
adelmann committed
29
#include <cmath>
gsell's avatar
gsell committed
30

31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
#include "Algorithms/PartPusher.h"
#include "AbsBeamline/AlignWrapper.h"
#include "AbsBeamline/BeamBeam.h"
#include "AbsBeamline/Collimator.h"
#include "AbsBeamline/Corrector.h"
#include "AbsBeamline/Diagnostic.h"
#include "AbsBeamline/Drift.h"
#include "AbsBeamline/ElementBase.h"
#include "AbsBeamline/Lambertson.h"
#include "AbsBeamline/Marker.h"
#include "AbsBeamline/Monitor.h"
#include "AbsBeamline/Multipole.h"
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/TravelingWave.h"
#include "AbsBeamline/RFQuadrupole.h"
#include "AbsBeamline/SBend.h"
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/ParallelPlate.h"
#include "AbsBeamline/CyclotronValley.h"
#include "Beamlines/Beamline.h"
#include "Lines/Sequence.h"
56 57
//--------- Added by Xiaoying Pang 04/22/2014 ---------------
#include "Solvers/CSRWakeFunction.hh"
gsell's avatar
gsell committed
58 59 60 61

#include "AbstractObjects/OpalData.h"

#include "BasicActions/Option.h"
kraus's avatar
kraus committed
62 63
#include "Utilities/OpalOptions.h"
#include "Utilities/Options.h"
gsell's avatar
gsell committed
64 65

#include "Distribution/Distribution.h"
66
#include "ValueDefinitions/RealVariable.h"
gsell's avatar
gsell committed
67 68
#include "Utilities/Timer.h"
#include "Utilities/OpalException.h"
69
#include "Solvers/SurfacePhysicsHandler.hh"
gsell's avatar
gsell committed
70 71 72 73 74 75 76 77 78 79
#include "Structure/BoundaryGeometry.h"

class PartData;

using namespace std;
using namespace OPALTimer;

ParallelTTracker::ParallelTTracker(const Beamline &beamline,
                                   const PartData &reference,
                                   bool revBeam,
adelmann's avatar
adelmann committed
80 81
                                   bool revTrack,
				   size_t N):
82 83 84 85 86 87 88 89 90 91 92 93
Tracker(beamline, reference, revBeam, revTrack),
itsBunch(NULL),
itsDataSink_m(NULL),
bgf_m(NULL),
itsOpalBeamline_m(),
lineDensity_m(),
RefPartR_zxy_m(0.0),
RefPartP_zxy_m(0.0),
RefPartR_suv_m(0.0),
RefPartP_suv_m(0.0),
globalEOL_m(false),
wakeStatus_m(false),
94 95
//--------- Added by Xiaoying Pang 04/22/2014 ---------------
wakeFunction_m(NULL),
96 97 98 99
surfaceStatus_m(false),
secondaryFlg_m(false),
mpacflg_m(true),
nEmissionMode_m(false),
100
zStop_m(),
101 102 103 104
scaleFactor_m(1.0),
vscaleFactor_m(scaleFactor_m),
recpGamma_m(1.0),
rescale_coeff_m(1.0),
105 106
dtCurrentTrack_m(0.0),
dtAllTracks_m(),
107
surfaceEmissionStop_m(-1),
108
specifiedNPart_m(N),
109 110 111 112 113 114 115
minStepforReBin_m(-1),
minBinEmitted_m(std::numeric_limits<size_t>::max()),
repartFreq_m(-1),
lastVisited_m(-1),
numRefs_m(-1),
gunSubTimeSteps_m(-1),
emissionSteps_m(std::numeric_limits<unsigned int>::max()),
116
localTrackSteps_m(),
117 118 119 120
maxNparts_m(0),
numberOfFieldEmittedParticles_m(std::numeric_limits<size_t>::max()),
bends_m(0),
numParticlesInSimulation_m(0),
kraus's avatar
kraus committed
121
space_orientation_m(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0),
122 123 124 125 126 127
timeIntegrationTimer1_m(IpplTimings::getTimer("TIntegration1")),
timeIntegrationTimer2_m(IpplTimings::getTimer("TIntegration2")),
timeFieldEvaluation_m(IpplTimings::getTimer("Fieldeval")),
BinRepartTimer_m(IpplTimings::getTimer("Binaryrepart")),
WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
Nimpact_m(0),
128 129
SeyNum_m(0.0),
sphys_m(NULL){
gsell's avatar
gsell committed
130 131 132 133 134 135 136 137
}

ParallelTTracker::ParallelTTracker(const Beamline &beamline,
                                   PartBunch &bunch,
                                   DataSink &ds,
                                   const PartData &reference,
                                   bool revBeam,
                                   bool revTrack,
138
                                   const std::vector<unsigned long long> &maxSteps,
139
                                   const std::vector<double> &zstop,
adelmann's avatar
adelmann committed
140
                                   int timeIntegrator,
141
                                   const std::vector<double> &dt,
adelmann's avatar
adelmann committed
142
				   size_t N):
143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162
Tracker(beamline, reference, revBeam, revTrack),
itsBunch(&bunch),
itsDataSink_m(&ds),
bgf_m(NULL),
itsOpalBeamline_m(),
lineDensity_m(),
RefPartR_zxy_m(0.0),
RefPartP_zxy_m(0.0),
RefPartR_suv_m(0.0),
RefPartP_suv_m(0.0),
globalEOL_m(false),
wakeStatus_m(false),
surfaceStatus_m(false),
secondaryFlg_m(false),
mpacflg_m(true),
nEmissionMode_m(false),
scaleFactor_m(itsBunch->getdT() * Physics::c),
vscaleFactor_m(scaleFactor_m),
recpGamma_m(1.0),
rescale_coeff_m(1.0),
163
dtCurrentTrack_m(0.0),
164
surfaceEmissionStop_m(-1),
165
specifiedNPart_m(N),
166 167 168 169 170 171 172 173 174 175 176
minStepforReBin_m(-1),
minBinEmitted_m(std::numeric_limits<size_t>::max()),
repartFreq_m(-1),
lastVisited_m(-1),
numRefs_m(-1),
gunSubTimeSteps_m(-1),
emissionSteps_m(numeric_limits<unsigned int>::max()),
maxNparts_m(0),
numberOfFieldEmittedParticles_m(numeric_limits<size_t>::max()),
bends_m(0),
numParticlesInSimulation_m(0),
kraus's avatar
kraus committed
177
space_orientation_m(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0),
178 179 180 181 182 183 184 185
timeIntegrationTimer1_m(IpplTimings::getTimer("TIntegration1")),
timeIntegrationTimer2_m(IpplTimings::getTimer("TIntegration2")),
timeFieldEvaluation_m(IpplTimings::getTimer("Fieldeval")),
BinRepartTimer_m(IpplTimings::getTimer("Binaryrepart")),
WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
timeIntegrator_m(timeIntegrator),
Nimpact_m(0),
SeyNum_m(0.0) {
186

187 188 189 190 191 192 193 194 195 196
    for (std::vector<unsigned long long>::const_iterator it = maxSteps.begin(); it != maxSteps.end(); ++ it) {
        localTrackSteps_m.push(*it);
    }
    for (std::vector<double>::const_iterator it = dt.begin(); it != dt.end(); ++ it) {
        dtAllTracks_m.push(*it);
    }
    for (std::vector<double>::const_iterator it = zstop.begin(); it != zstop.end(); ++ it) {
        zStop_m.push(*it);
    }

197
    //    itsBeamline = dynamic_cast<Beamline*>(beamline.clone());
198

gsell's avatar
gsell committed
199 200
}

201 202 203 204 205 206 207
#ifdef HAVE_AMR_SOLVER
ParallelTTracker::ParallelTTracker(const Beamline &beamline,
                                   PartBunch &bunch,
                                   DataSink &ds,
                                   const PartData &reference,
                                   bool revBeam,
                                   bool revTrack,
208
                                   const std::vector<unsigned long long> &maxSteps,
209
                                   const std::vector<double> &zstop,
210
                                   int timeIntegrator,
211
                                   const std::vector<double> &dt,
212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233
				   size_t N,
				   Amr* amrptr_in):
Tracker(beamline, reference, revBeam, revTrack),
itsBunch(&bunch),
itsDataSink_m(&ds),
bgf_m(NULL),
itsOpalBeamline_m(),
lineDensity_m(),
RefPartR_zxy_m(0.0),
RefPartP_zxy_m(0.0),
RefPartR_suv_m(0.0),
RefPartP_suv_m(0.0),
globalEOL_m(false),
wakeStatus_m(false),
surfaceStatus_m(false),
secondaryFlg_m(false),
mpacflg_m(true),
nEmissionMode_m(false),
scaleFactor_m(itsBunch->getdT() * Physics::c),
vscaleFactor_m(scaleFactor_m),
recpGamma_m(1.0),
rescale_coeff_m(1.0),
234
dtCurrentTrack_m(0.0),
235
surfaceEmissionStop_m(-1),
236
specifiedNPart_m(N),
237 238 239 240 241 242 243 244 245 246 247
minStepforReBin_m(-1),
minBinEmitted_m(std::numeric_limits<size_t>::max()),
repartFreq_m(-1),
lastVisited_m(-1),
numRefs_m(-1),
gunSubTimeSteps_m(-1),
emissionSteps_m(numeric_limits<unsigned int>::max()),
maxNparts_m(0),
numberOfFieldEmittedParticles_m(numeric_limits<size_t>::max()),
bends_m(0),
numParticlesInSimulation_m(0),
kraus's avatar
kraus committed
248
space_orientation_m(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0),
249 250 251 252 253 254 255 256 257
timeIntegrationTimer1_m(IpplTimings::getTimer("TIntegration1")),
timeIntegrationTimer2_m(IpplTimings::getTimer("TIntegration2")),
timeFieldEvaluation_m(IpplTimings::getTimer("Fieldeval")),
BinRepartTimer_m(IpplTimings::getTimer("Binaryrepart")),
WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
timeIntegrator_m(timeIntegrator),
Nimpact_m(0),
SeyNum_m(0.0),
amrptr(amrptr_in){
258 259 260 261 262 263 264 265 266 267

    for (std::vector<unsigned long long>::const_iterator it = maxSteps.begin(); it != maxSteps.end(); ++ it) {
        localTrackSteps_m.push(*it);
    }
    for (std::vector<double>::const_iterator it = dt.begin(); it != dt.end(); ++ it) {
        dtAllTracks_m.push(*it);
    }
    for (std::vector<double>::const_iterator it = zstop.begin(); it != zstop.end(); ++ it) {
        zStop_m.push(*it);
    }
268 269
}
#endif
gsell's avatar
gsell committed
270 271

ParallelTTracker::~ParallelTTracker() {
272

gsell's avatar
gsell committed
273 274
}

275 276 277
void ParallelTTracker::applyEntranceFringe(double angle, double curve,
                                           const BMultipoleField &field, double scale) {
}
278 279


280 281 282
void ParallelTTracker::applyExitFringe(double angle, double curve,
                                       const BMultipoleField &field, double scale) {
}
283 284


285
void ParallelTTracker::checkCavity(double s, Component *& comp, double & cavity_start_pos) {
286

287 288 289 290 291 292 293 294 295 296 297
    comp = NULL;
    for(FieldList::iterator fit = cavities_m.begin(); fit != cavities_m.end(); ++ fit) {
        if((fit != currently_ap_cavity_m)
           && ((*fit).getStart() <= s) && (s <= (*fit).getEnd())) {
            comp = (*fit).getElement().get();
            cavity_start_pos = (*fit).getStart();
            currently_ap_cavity_m = fit;
            return;
        }
    }
}
298

299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317
void ParallelTTracker::updateRFElement(std::string elName, double maxPhi) {
    /**
     The maximum phase is added to the nominal phase of
     the element. This is done on all nodes except node 0 where
     the Autophase took place.
     */
    double phi  = 0.0;
    for(FieldList::iterator fit = cavities_m.begin(); fit != cavities_m.end(); ++ fit) {
        if((*fit).getElement()->getName() == elName) {
            if((*fit).getElement()->getType() == "TravelingWave") {
                phi  =  static_cast<TravelingWave *>((*fit).getElement().get())->getPhasem();
                phi += maxPhi;
                static_cast<TravelingWave *>((*fit).getElement().get())->updatePhasem(phi);
            } else {
                phi  = static_cast<RFCavity *>((*fit).getElement().get())->getPhasem();
                phi += maxPhi;
                static_cast<RFCavity *>((*fit).getElement().get())->updatePhasem(phi);
            }
        }
gsell's avatar
gsell committed
318 319 320
    }
}

321

322 323 324 325 326
void ParallelTTracker::updateAllRFElements(double phiShift) {
    /**
     All RF-Elements gets updated, where the phiShift is the
     global phase shift in units of seconds.
     */
327
    Inform msg("ParallelTTracker ");
328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356
    double phi = 0;
    double freq = 0.0;
    const double RADDEG = 1.0 / Physics::pi * 180.0;
    //   Inform m ("updateALLRFElements ",INFORM_ALL_NODES);
    msg << "\n-------------------------------------------------------------------------------------\n";
    for(FieldList::iterator fit = cavities_m.begin(); fit != cavities_m.end(); ++ fit) {
        if(fit != cavities_m.begin())
            msg << "\n";
        if((*fit).getElement()->getType() == "TravelingWave") {
            freq = static_cast<TravelingWave *>((*fit).getElement().get())->getFrequencym();
            phi = static_cast<TravelingWave *>((*fit).getElement().get())->getPhasem();
            msg << (*fit).getElement()->getName()
            << ": phi= phi_nom + phi_maxE + global phase shift= " << (phi*RADDEG)-(phiShift*freq*RADDEG) << " degree, "
            << "(global phase shift= " << -phiShift *freq *RADDEG << " degree)\n";
            phi -= (phiShift * freq);
            static_cast<TravelingWave *>((*fit).getElement().get())->updatePhasem(phi);
        } else {
            freq = static_cast<RFCavity *>((*fit).getElement().get())->getFrequencym();
            phi = static_cast<RFCavity *>((*fit).getElement().get())->getPhasem();
            msg << (*fit).getElement()->getName()
            << ": phi= phi_nom + phi_maxE + global phase shift= " << (phi*RADDEG)-(phiShift*freq*RADDEG) << " degree, "
            << "global phase shift= " << -phiShift *freq *RADDEG << " degree\n";
            phi -= (phiShift * freq);
            static_cast<RFCavity *>((*fit).getElement().get())->updatePhasem(phi);
        }
    }
    msg << "-------------------------------------------------------------------------------------\n"
	<< endl;
}
357 358


359 360
FieldList ParallelTTracker::executeAutoPhaseForSliceTracker() {
    Inform msg("executeAutoPhaseForSliceTracker ");
361

362
    double gPhaseSave;
363

364 365
    gPhaseSave = OpalData::getInstance()->getGlobalPhaseShift();
    OpalData::getInstance()->setGlobalPhaseShift(0.0);
366

367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382
    itsBeamline_m.accept(*this);
    // make sure that no monitor has overlap with two tracks
    FieldList monitors = itsOpalBeamline_m.getElementByType("Monitor");
    for(FieldList::iterator it = monitors.begin(); it != monitors.end(); ++ it) {
        double zbegin, zend;
        it->getElement()->getDimensions(zbegin, zend);
        if(zbegin < zStop_m.front() && zend >= zStop_m.front()) {
            msg << "\033[0;31m"
            << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
            << "% Removing '" << it->getElement()->getName() << "' since it resides in two tracks.   %\n"
            << "% Please adjust zstop or place your monitor at a different position to prevent this. %\n "
            << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
            << "\033[0m"
            << endl;
            static_cast<Monitor *>(it->getElement().get())->moveBy(-zend - 0.001);
            itsOpalBeamline_m.removeElement(it->getElement()->getName());
gsell's avatar
gsell committed
383
        }
384 385
    }
    itsOpalBeamline_m.prepareSections();
386

387 388 389 390
    cavities_m = itsOpalBeamline_m.getElementByType("RFCavity");
    currently_ap_cavity_m = cavities_m.end();
    FieldList travelingwaves = itsOpalBeamline_m.getElementByType("TravelingWave");
    cavities_m.merge(travelingwaves, ClassicField::SortAsc);
391

392 393 394 395 396 397 398 399 400 401
    int tag = 101;
    int Parent = 0;
    Vector_t iniR(0.0);
    Vector_t iniP(0.0, 0.0, 1E-6);
    PID_t id;
    Ppos_t r, p, x;
    ParticleAttrib<double> q, dt;
    ParticleAttrib<int> bin;
    ParticleAttrib<long> ls;
    ParticleAttrib<short> ptype;
402

403
    double zStop = itsOpalBeamline_m.calcBeamlineLenght();
404

405
    msg << "Preparation done zstop= " << zStop << endl;
406

407 408 409
    if(Ippl::myNode() == 0)
      itsBunch->create(1);
    itsBunch->update();
410

411 412 413 414 415 416 417
    if(Ippl::myNode() == 0) {
      itsBunch->R[0] = iniR;
      itsBunch->P[0] = iniP;
      itsBunch->Bin[0] = 0;
      itsBunch->Q[0] = itsBunch->getChargePerParticle();
      itsBunch->PType[0] = 0;
      itsBunch->LastSection[0] = 0;
418

419 420
      RefPartP_suv_m = iniP;
      RefPartR_suv_m = iniR;
gsell's avatar
gsell committed
421
    }
422

423 424 425 426 427 428 429 430
    updateSpaceOrientation(false);
    executeAutoPhase(Options::autoPhase, zStop);
    if(Ippl::myNode() == 0) {
      RefPartP_suv_m = itsBunch->P[0];
      // need to rebuild for updateAllRFElements
      cavities_m = itsOpalBeamline_m.getElementByType("RFCavity");
      travelingwaves = itsOpalBeamline_m.getElementByType("TravelingWave");
      cavities_m.merge(travelingwaves, ClassicField::SortAsc);
431

432 433 434 435
      // now send all max phases and names of the cavities to
      // all the other nodes for updating.
      Message *mess = new Message();
      putMessage(*mess, OpalData::getInstance()->getNumberOfMaxPhases());
436

437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458
      for(std::vector<MaxPhasesT>::iterator it = OpalData::getInstance()->getFirstMaxPhases(); it < OpalData::getInstance()->getLastMaxPhases(); it++) {
	putMessage(*mess, (*it).first);
	putMessage(*mess, (*it).second);
      }
      Ippl::Comm->broadcast_all(mess, tag);

      delete mess;
    } else {
      // receive max phases and names and update the structures
      int nData = 0;
      Message *mess = Ippl::Comm->receive_block(Parent, tag);
      getMessage(*mess, nData);
      for(int i = 0; i < nData; i++) {
	std::string elName;
	double maxPhi;
	getMessage(*mess, elName);
	getMessage(*mess, maxPhi);
	updateRFElement(elName, maxPhi);
	OpalData::getInstance()->setMaxPhase(elName, maxPhi);
      }

      delete mess;
gsell's avatar
gsell committed
459
    }
460

461 462 463
    if(Ippl::myNode() == 0)
      itsBunch->destroy(1, 0);
    itsBunch->update();
464

465 466 467
    OpalData::getInstance()->setGlobalPhaseShift(gPhaseSave);
    return cavities_m;
}
468 469


470 471
void ParallelTTracker::executeAutoPhase(int numRefs, double zStop) {
    Inform msg("Autophasing ");
472

473
    const double RADDEG = 180.0 / Physics::pi;
474

475
    Vector_t rmin, rmax;
476

477 478 479
    std::queue<double> dtAutoPhasing(dtAllTracks_m);
    std::queue<double> maxSPosAutoPhasing(zStop_m);
    std::queue<unsigned long long> maxStepsAutoPhasing(localTrackSteps_m);
480

481
    maxSPosAutoPhasing.back() = zStop;
482

483 484 485 486
    size_t step = 0;
    double tSave = itsBunch->getT();
    double dTSave = itsBunch->getdT();
    double scaleFactorSave = scaleFactor_m;
487

488 489 490 491 492 493
    const int dtfraction = 2;
    double newDt = itsBunch->getdT() / dtfraction;
    itsBunch->setdT(newDt);
    itsBunch->dt = newDt;
    scaleFactor_m = newDt * Physics::c;
    vscaleFactor_m = Vector_t(scaleFactor_m);
494

495
    bends_m = 0;  //fixme: AP and bends will not work at the moment
496

497
    BorisPusher pusher(itsReference);
498

499 500 501 502 503 504 505
    for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) {
        long &l = itsBunch->LastSection[i];
        l = -1;
        itsOpalBeamline_m.getSectionIndexAt(itsBunch->R[i], l);
        itsBunch->ResetLocalCoordinateSystem(i, itsOpalBeamline_m.getOrientation(l),
                                             itsOpalBeamline_m.getSectionStart(l));
    }
506

507 508
    RefPartR_suv_m = RefPartR_zxy_m = rmin = rmax = itsBunch->R[0];
    RefPartP_suv_m = RefPartP_zxy_m = itsBunch->P[0];
509

510 511 512 513 514 515
    /* Activate all elements which influence the particles when the simulation starts;
     * mark all elements which are already past.
     *
     * Increase margin from 3.*c*dt to 10.*c*dt to prevent that fieldmaps are accessed
     * before they are allocated when increasing the timestep in the gun.
     */
516

517
    double margin = 10. * RefPartP_suv_m(2) * scaleFactor_m / sqrt(1.0 + dot(RefPartP_suv_m, RefPartP_suv_m));
518

519
    margin = 0.01 > margin ? 0.01 : margin;
520

521
    itsOpalBeamline_m.switchElements(rmin(2) - margin, rmax(2) + margin, true);
522

523 524
    double cavity_start = 0.0;
    Component *cavity = NULL;
525

526 527 528 529 530 531 532
    while (maxStepsAutoPhasing.size() > 0) {
        maxStepsAutoPhasing.front() = maxStepsAutoPhasing.front() * dtfraction + step;
        newDt = dtAutoPhasing.front() / dtfraction;
        itsBunch->setdT(newDt);
        itsBunch->dt = newDt;
        scaleFactor_m = newDt * Physics::c;
        vscaleFactor_m = Vector_t(scaleFactor_m);
533

534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614
        msg << "\n"
            << "**********************************************\n"
            << "new set of dt, zstop and max number of time steps\n"
            << "**********************************************\n\n"
            << "start at t = " << itsBunch->getT() << " [s], "
            << "zstop at z = " << maxSPosAutoPhasing.front() << " [m],\n "
            << "dt = " << itsBunch->dt[0] << " [s], "
            << "step = " << step << ", "
            << "R =  " << itsBunch->R[0] << " [m]" << endl;

        for(; step < maxStepsAutoPhasing.front(); ++step) {

            if (itsBunch->getLocalNum() != 0) {
                // let's do a drifting step to probe if the particle will reach element in next step
                Vector_t R_drift = itsBunch->R[0] + itsBunch->P[0] / sqrt(1.0 + dot(itsBunch->P[0],
                                                                                    itsBunch->P[0])) * vscaleFactor_m;
                checkCavity(R_drift[2], cavity, cavity_start);
            }
            else
                cavity = NULL;

            if(cavity != NULL) {
                double orig_phi = 0.0;
                double Phimax = 0.0, Emax = 0.0;
                double PhiAstra;
                //////
                const double beta = sqrt(1. - 1 / (itsBunch->P[0](2) * itsBunch->P[0](2) + 1.));
                const double tErr  = (cavity_start - itsBunch->R[0](2)) / (Physics::c * beta);

                bool apVeto;

                INFOMSG("Found " << cavity->getName()
                        << " at " << itsBunch->R[0](2) << " [m], "
                        << "step  " << step << ", "
                        << "t= " << itsBunch->getT() << " [s],\n"
                        << "E= " << getEnergyMeV(itsBunch->P[0]) << " [MeV]\n"
                        << "start phase scan ... " << endl);


                INFOMSG("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
                if(cavity->getType() == "TravelingWave") {
                    orig_phi = static_cast<TravelingWave *>(cavity)->getPhasem();
                    apVeto = static_cast<TravelingWave *>(cavity)->getAutophaseVeto();
                    if(apVeto) {
                        msg << " ----> APVETO -----> "
                            << static_cast<TravelingWave *>(cavity)->getName() <<  endl;
                        Phimax = orig_phi;
                    }
                    INFOMSG(cavity->getName() << ", "
                            << "start Ekin= " << getEnergyMeV(itsBunch->P[0]) << " MeV, "
                            << "t= " << itsBunch->getT() << " s, "
                            << "phi= " << orig_phi << ", " << endl;);

                    if(!apVeto) {
                        TravelingWave *element = static_cast<TravelingWave *>(cavity);
                        Phimax = element->getAutoPhaseEstimate(getEnergyMeV(itsBunch->P[0]),
                                                               itsBunch->getT() + tErr,
                                                               itsReference.getQ(),
                                                               itsReference.getM() * 1e-6);
                    }
                } else {
                    orig_phi = static_cast<RFCavity *>(cavity)->getPhasem();
                    apVeto = static_cast<RFCavity *>(cavity)->getAutophaseVeto();
                    if(apVeto) {
                        msg << " ----> APVETO -----> "
                            << static_cast<RFCavity *>(cavity)->getName() << endl;
                        Phimax = orig_phi;
                    }
                    INFOMSG(cavity->getName() << ", "
                            << "start Ekin= " << getEnergyMeV(itsBunch->P[0]) << " MeV, "
                            << "t= " << itsBunch->getT() << " s, "
                            << "phi= " << orig_phi << ", " << endl;);

                    if(!apVeto) {
                        RFCavity *element = static_cast<RFCavity *>(cavity);
                        Phimax = element->getAutoPhaseEstimate(getEnergyMeV(itsBunch->P[0]),
                                                               itsBunch->getT() + tErr,
                                                               itsReference.getQ(),
                                                               itsReference.getM() * 1e-6);
                    }
                }
615

616 617 618 619
                double Phiini = Phimax;
                double phi = Phiini;
                double dphi = Physics::pi / 360.0;
                int j = -1;
620

621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637
                double E = APtrack(cavity, cavity_start, phi);
                if(!apVeto) {
                    msg << "Did APtrack with phi= " << phi << " result E= " << E << endl;
                    //                INFOMSG("do fine scan around effective max energy (" << E << " MeV)" << ", dphi= " << dphi << endl;);
                    do {
                        j ++;
                        Emax = E;
                        Phiini = phi;
                        phi -= dphi;
                        //                    INFOMSG("try phi= " << phi << " rad -> DEkin= ";);
                        E = APtrack(cavity, cavity_start, phi);
                        //                    if(E > Emax) {
                        //                        INFOMSG(E - Emax << " MeV: accepted" << endl;);
                        //  } else {
                        //                        INFOMSG(E - Emax << " MeV: rejected" << " E= " << E << " Emax= " << Emax << endl;);
                        // }
                    } while(E > Emax);
638

639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683
                    if(j == 0) {
                        phi = Phiini;
                        E = Emax;
                        j = -1;
                        do {
                            j ++;
                            Emax = E;
                            Phiini = phi;
                            phi += dphi;
                            //                        INFOMSG("try phi= " << phi << " rad -> DEkin= ";);
                            E = APtrack(cavity, cavity_start, phi);
                            //                        if(E > Emax) {
                            //  INFOMSG(E - Emax << " MeV: accepted" << endl;);
                            //  } else {
                            //  INFOMSG(E - Emax << " MeV: rejected" << endl;);
                            // }
                        } while(E > Emax);
                    }
                    for(int refinement_level = 0; refinement_level < numRefs; refinement_level ++) {
                        dphi /= 2.;
                        //                    INFOMSG("refinement level: " << refinement_level + 1 << ", dphi= " << dphi << endl;);
                        phi = Phiini - dphi;
                        // INFOMSG("try phi= " << phi << " rad -> DEkin= ";);
                        E = APtrack(cavity, cavity_start, phi);
                        if(E > Emax) {
                            //  INFOMSG(E - Emax << " MeV: accepted" << endl;);
                            Phiini = phi;
                            Emax = E;
                        } else {
                            // INFOMSG(E - Emax << " MeV: rejected" << endl;);
                            phi = Phiini + dphi;
                            // INFOMSG("try phi= " << phi << " rad -> DEkin= ";);
                            E = APtrack(cavity, cavity_start, phi);
                            if(E > Emax) {
                                //  INFOMSG(E - Emax << " MeV: accepted" << endl;);
                                Phiini = phi;
                                Emax = E;
                            }
                            //else {
                            //  INFOMSG(E - Emax << " MeV: rejected" << endl;);
                            //}
                        }
                    }
                    Phimax = Phiini;
                    phi = Phimax + orig_phi;
684

685 686 687 688 689 690 691 692 693 694 695 696 697
                    INFOMSG("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \n");
                } else {
                    msg << "Tracking with phi= " << orig_phi << " result E= " << E << endl;
                    phi = orig_phi;
                    Emax = E;
                    E = APtrack(cavity, cavity_start, phi - Physics::pi);
                    msg << "Tracking with phi= " << orig_phi - Physics::pi << " result E= " << E << endl;
                    if (E > Emax) {
                        phi = orig_phi - Physics::pi;
                        Phimax = phi;
                        Emax = E;
                    }
                }
698

699 700 701 702 703
                if(cavity->getType() == "TravelingWave") {
                    static_cast<TravelingWave *>(cavity)->updatePhasem(phi);
                } else {
                    static_cast<RFCavity *>(cavity)->updatePhasem(phi);
                }
704

705 706
                PhiAstra = (Phimax * RADDEG) + 90.0;
                PhiAstra -= floor(PhiAstra / 360.) * 360.;
707

708 709 710
                msg << cavity->getName() << "_phi= "  << Phimax << " rad / "
                    << Phimax *RADDEG <<  " deg, AstraPhi= " << PhiAstra << " deg,\n"
                    << "E= " << Emax << " (MeV), " << "phi_nom= " << orig_phi *RADDEG << endl;
711

712
                OpalData::getInstance()->setMaxPhase(cavity->getName(), Phimax);
gsell's avatar
gsell committed
713
            }
714

715 716 717 718
            // --- doOneStep(pusher);
            // --- timeIntegration1(pusher);
            IpplTimings::startTimer(timeIntegrationTimer1_m);
            double t = itsBunch->getT();
719

720 721
            // increase margin from 3.*c*dt to 10.*c*dt to prevent that fieldmaps are accessed
            // before they are allocated when increasing the timestep in the gun.
722

723 724
            switchElements(10.0);
            itsOpalBeamline_m.resetStatus();
725

726 727 728 729 730

            for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) {
                itsBunch->R[i] /= Vector_t(Physics::c * itsBunch->getdT());
                pusher.push(itsBunch->R[i], itsBunch->P[i], itsBunch->dt[i]);
                itsBunch->R[i] *= Vector_t(Physics::c * itsBunch->getdT());
gsell's avatar
gsell committed
731
            }
732 733
            IpplTimings::stopTimer(timeIntegrationTimer1_m);
            // --- timeIntegration1(pusher);
734

735 736
            itsBunch->Ef = Vector_t(0.0);
            itsBunch->Bf = Vector_t(0.0);
737

738 739 740 741 742 743 744 745 746 747 748
            // --- computeExternalFields();
            IpplTimings::startTimer(timeFieldEvaluation_m);
            for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) {
                long ls = itsBunch->LastSection[i];
                itsOpalBeamline_m.getSectionIndexAt(itsBunch->R[i], ls);
                if(ls != itsBunch->LastSection[i]) {
                    if(!itsOpalBeamline_m.section_is_glued_to(itsBunch->LastSection[i], ls)) {
                        itsBunch->ResetLocalCoordinateSystem(i, itsOpalBeamline_m.getOrientation(ls), itsOpalBeamline_m.getSectionStart(ls));
                    }
                    itsBunch->LastSection[i] = ls;
                }
749

750 751 752 753 754 755 756 757 758
                Vector_t externalE, externalB;
                const unsigned long rtv = itsOpalBeamline_m.getFieldAt(i, itsBunch->R[i], ls, itsBunch->getT() + itsBunch->dt[i] / 2., externalE, externalB);
                if (rtv)
                    ;
                itsBunch->Ef[i] += externalE;
                itsBunch->Bf[i] += externalB;
            }
            IpplTimings::stopTimer(timeFieldEvaluation_m);
            // --- computeExternalFields
759

760 761 762 763
            // --- timeIntegration2(pusher);
            IpplTimings::startTimer(timeIntegrationTimer2_m);
            double recpgamma = 1.0 / sqrt(1.0 + dot(RefPartP_suv_m, RefPartP_suv_m));
            RefPartR_zxy_m += RefPartP_zxy_m * recpgamma / 2. * scaleFactor_m;
764

765 766 767 768
            // kickParticlesAutophase(pusher);
            for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) {
                pusher.kick(itsBunch->R[i], itsBunch->P[i], itsBunch->Ef[i], itsBunch->Bf[i], itsBunch->dt[i]);
            }
769

770
            itsBunch->calcBeamParametersLight();
771

772 773 774 775 776 777 778 779 780 781
            for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) {
                itsBunch->R[i] /= Vector_t(Physics::c * itsBunch->getdT());
                pusher.push(itsBunch->R[i], itsBunch->P[i], itsBunch->dt[i]);
                itsBunch->R[i] *= Vector_t(Physics::c * itsBunch->getdT());
            }
            t += itsBunch->getdT();
            itsBunch->setT(t);
            IpplTimings::stopTimer(timeIntegrationTimer2_m);
            // --- timeIntegration2(pusher);
            // --- doOneStep(pusher);
782

783
            double sposRef	= 0.0;
784

785
            if(!(step % 5000))	{
786

787 788
                if (itsBunch->getLocalNum() != 0)
                    sposRef = itsBunch->R[0](2);
789

790
                reduce(sposRef,sposRef,OpAddAssign());
791

792 793
                if(sposRef > maxSPosAutoPhasing.front())
                    maxStepsAutoPhasing.front() = step;
794 795


796 797
                INFOMSG("step = " << step << ", spos = " << sposRef << " [m], t= " << itsBunch->getT() << " [s], "
                        << "E= " << itsBunch->get_meanEnergy() << " [MeV] " << endl);
gsell's avatar
gsell committed
798 799
            }
        }
800

801 802 803
        dtAutoPhasing.pop();
        maxSPosAutoPhasing.pop();
        maxStepsAutoPhasing.pop();
gsell's avatar
gsell committed
804 805
    }

806 807 808 809 810 811
    FieldList sbends = itsOpalBeamline_m.getElementByType("SBend");
    for (FieldList::iterator it = sbends.begin(); it != sbends.end(); ++ it) {
        SBend* bend = static_cast<SBend*>(it->getElement().get());
        bend->resetReinitializeFlag();
        bend->resetRecalcRefTrajFlag();
    }
gsell's avatar
gsell committed
812

813 814 815 816 817
    FieldList rbends = itsOpalBeamline_m.getElementByType("RBend");
    for (FieldList::iterator it = rbends.begin(); it != rbends.end(); ++ it) {
        RBend* bend = static_cast<RBend*>(it->getElement().get());
        bend->resetReinitializeFlag();
        bend->resetRecalcRefTrajFlag();
818
    }
gsell's avatar
gsell committed
819

820 821 822 823
    scaleFactor_m = scaleFactorSave;
    itsBunch->setT(tSave);
    itsBunch->setdT(dTSave);
}
824

825 826 827
double ParallelTTracker::APtrack(Component *cavity, double cavity_start_pos, const double &phi) const {
    double beta = std::max(sqrt(1. - 1 / (itsBunch->P[0](2) * itsBunch->P[0](2) + 1.)), 0.0001);
    double tErr  = (cavity_start_pos - itsBunch->R[0](2)) / (Physics::c * beta);
828

829
    //    INFOMSG("beta = " << beta << " tErr = " << tErr << endl;);
830

831 832 833 834 835 836 837 838 839 840 841 842 843
    double finalMomentum = 0.0;
    if(cavity->getType() == "TravelingWave") {
        TravelingWave *tws = static_cast<TravelingWave *>(cavity);
        tws->updatePhasem(phi);
        std::pair<double, double> pe = tws->trackOnAxisParticle(itsBunch->P[0](2),
                                                                itsBunch->getT() + tErr,
                                                                itsBunch->dt[0],
                                                                itsBunch->getQ(),
                                                                itsBunch->getM() * 1e-6);
        finalMomentum = pe.first;
    } else {
        RFCavity *rfc = static_cast<RFCavity *>(cavity);
        rfc->updatePhasem(phi);
844

845 846 847 848 849 850 851 852 853 854 855
        std::pair<double, double> pe = rfc->trackOnAxisParticle(itsBunch->P[0](2),
                                                                itsBunch->getT() + tErr,
                                                                itsBunch->dt[0],
                                                                itsBunch->getQ(),
                                                                itsBunch->getM() * 1e-6);
        finalMomentum = pe.first;
    }
    double finalGamma = sqrt(1.0 + finalMomentum * finalMomentum);
    double finalKineticEnergy = (finalGamma - 1.0) * itsBunch->getM() * 1e-6;
    return finalKineticEnergy;
}
856 857


858
double ParallelTTracker::getGlobalPhaseShift() {
859

860 861 862 863 864 865
    if(Options::autoPhase > 0) {
        double gPhaseSave = OpalData::getInstance()->getGlobalPhaseShift();
        OpalData::getInstance()->setGlobalPhaseShift(0.0);
        return gPhaseSave;
    } else
        return 0;
866

867
}
868

869 870
void ParallelTTracker::doAutoPhasing() {
    typedef std::vector<MaxPhasesT>::iterator iterator_t;
871

872
    if(Options::autoPhase == 0) return;
873

874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893
    if(OpalData::getInstance()->inRestartRun()) {
        itsDataSink_m->retriveCavityInformation(OpalData::getInstance()->getInputFn());
        iterator_t it = OpalData::getInstance()->getFirstMaxPhases();
        iterator_t end = OpalData::getInstance()->getLastMaxPhases();
        for(; it < end; ++ it)
            updateRFElement((*it).first, (*it).second);
    } else {
        if(OpalData::getInstance()->hasBunchAllocated()) {
            // we are in a followup track and the phase information is
            // already stored in the OPAL dictionary.
            iterator_t it = OpalData::getInstance()->getFirstMaxPhases();
            iterator_t end = OpalData::getInstance()->getLastMaxPhases();
            for(; it < end; ++ it) {
                updateRFElement((*it).first, (*it).second);
                INFOMSG("In follow-up track use saved phases for -> name: " << (*it).first
                        << " phi= " << (*it).second << " (rad)" << endl);
            }
        } else {
            int tag = 101;
            int Parent = 0;
894

895 896 897 898
            Vector_t stashedP = itsBunch->get_pmean();
            Vector_t initialR(0.0);
            // TODO: get correct thermal energy of particles for Astra model
            Vector_t initialP(0.0, 0.0, std::max(stashedP(2), 1e-6));
gsell's avatar
gsell committed
899

900 901 902 903
            if (itsBunch->getTotalNum() > 0) { // we are not emiting otherwise there wouldn't be any particles yet
                Vector_t stashedR = itsBunch->get_rmean();
                initialR(2) = stashedR(2);
            }
gsell's avatar
gsell committed
904

905 906 907 908 909 910 911 912 913 914 915 916
            itsBunch->stash();
            double zStop = itsOpalBeamline_m.calcBeamlineLenght();
            if(Ippl::myNode() == 0) {
                itsBunch->create(1);
                itsBunch->R[0] = initialR;
                itsBunch->P[0] = initialP;
                itsBunch->Bin[0] = 0;
                itsBunch->Q[0] = itsBunch->getChargePerParticle();
                itsBunch->PType[0] = 0;
                itsBunch->LastSection[0] = 0;
            }
            itsBunch->update();
gsell's avatar
gsell committed
917

918
            executeAutoPhase(Options::autoPhase, zStop);
gsell's avatar
gsell committed
919

920 921 922 923
            // the single particle may change the node during update
            Parent = (itsBunch->getLocalNum() == 1? Ippl::myNode(): -1);
            reduce(Parent, Parent, OpAddAssign());
            Parent += (Ippl::getNodes() - 1);
gsell's avatar
gsell committed
924

925 926 927 928
            if (Ippl::myNode() == Parent) {    // now send all max phases and names of the cavities to
                // all the other nodes for updating.
                Message *mess = new Message();
                putMessage(*mess, OpalData::getInstance()->getNumberOfMaxPhases());
gsell's avatar
gsell committed
929

930 931 932 933 934 935 936
                iterator_t it = OpalData::getInstance()->getFirstMaxPhases();
                iterator_t end = OpalData::getInstance()->getLastMaxPhases();
                for(; it < end; ++ it) {
                    putMessage(*mess, (*it).first);
                    putMessage(*mess, (*it).second);
                }
                Ippl::Comm->broadcast_all(mess, tag);
937

938
                itsBunch->destroy(1, 0);
939

940
                delete mess;
gsell's avatar
gsell committed
941
            } else {
942 943 944 945 946 947 948 949 950 951 952 953 954 955
                // receive max phases and names and update the structure
                int nData = 0;
                Message *mess = Ippl::Comm->receive_block(Parent, tag);
                getMessage(*mess, nData);
                for(int i = 0; i < nData; i++) {
                    std::string elName;
                    double maxPhi;
                    getMessage(*mess, elName);
                    getMessage(*mess, maxPhi);
                    updateRFElement(elName, maxPhi);
                    OpalData::getInstance()->setMaxPhase(elName, maxPhi);
                }

                delete mess;
gsell's avatar
gsell committed
956
            }
957 958 959
            itsBunch->update();
            itsBunch->pop();
            itsDataSink_m->storeCavityInformation();
gsell's avatar
gsell committed
960 961
        }
    }
962 963
    updateAllRFElements(OpalData::getInstance()->getGlobalPhaseShift());
    INFOMSG("finished autophasing" << endl);
gsell's avatar
gsell committed
964 965
}

966 967 968 969 970 971 972 973
void ParallelTTracker::execute() {
#ifdef HAVE_AMR_SOLVER
    executeAMRTracker();
#else
    if(timeIntegrator_m == 3) {
        executeAMTSTracker();
    } else {
        executeDefaultTracker();
gsell's avatar
gsell committed
974
    }
975
#endif
gsell's avatar
gsell committed
976 977
}

978 979 980 981 982 983
void ParallelTTracker::executeDefaultTracker() {
    Inform msg("ParallelTTracker ");
    const Vector_t vscaleFactor_m = Vector_t(scaleFactor_m);
    BorisPusher pusher(itsReference);
    secondaryFlg_m = false;
    dtCurrentTrack_m = itsBunch->getdT();
gsell's avatar
gsell committed
984

985 986 987 988
    // upper limit of particle number when we do field emission and secondary emission
    // simulation. Could be reset to another value in input file with MAXPARTSNUM.
    maxNparts_m = 100000000;
    nEmissionMode_m = true;
989

990
    prepareSections();
991

992 993 994 995 996 997 998
    if (OpalData::getInstance()->hasBunchAllocated()) {
        // delete last entry of sdds file and load balance file
        // if we are in a follow-up track
        itsDataSink_m->rewindLinesSDDS(1);
        itsDataSink_m->rewindLinesLBal(1);
    }

999 1000
    // do autophasing before tracking without a global phase shift!
    doAutoPhasing();
1001

1002
    numParticlesInSimulation_m = itsBunch->getTotalNum();
1003

1004
    OPALTimer::Timer myt1;
1005

1006
    setTime();
1007

1008
    double t = itsBunch->getT();
1009

1010
    unsigned long long step = itsBunch->getLocalTrackStep();
1011

1012
    msg << "Track start at: " << myt1.time() << ", t= " << t << "; zstop at: " << zStop_m.front() << " [m]" << endl;
1013

1014 1015
    gunSubTimeSteps_m = 10;
    prepareEmission();
1016

1017
    doSchottyRenormalization();
1018

1019 1020 1021
    msg << "Executing ParallelTTracker, initial DT " << itsBunch->getdT() << " [s];\n"
    << "max integration steps " << localTrackSteps_m.front() << ", next step= " << step << endl;
    msg << "Using default (Boris-Buneman) integrator" << endl;
1022

1023 1024 1025 1026
    if (Options::info)
      itsOpalBeamline_m.print(msg);
    else
      msg << "Silent track ... " << endl;
1027

1028
    setupSUV(!(OpalData::getInstance()->inRestartRun() || (OpalData::getInstance()->hasBunchAllocated() && !Options::scan)));
1029

1030 1031 1032
    // increase margin from 3.*c*dt to 10.*c*dt to prevent that fieldmaps are accessed
    // before they are allocated when increasing the timestep in the gun.
    switchElements(10.0);
1033

1034
    initializeBoundaryGeometry();
1035

1036
    setOptionalVariables();
gsell's avatar
gsell committed
1037

1038 1039 1040
    // there is no point to do repartitioning with one node
    if(Ippl::getNodes() == 1)
        repartFreq_m = 1000000;
gsell's avatar
gsell committed
1041

1042 1043
    wakeStatus_m = false;
    surfaceStatus_m = false;
1044

1045 1046 1047 1048
    while (localTrackSteps_m.size() > 0) {
        localTrackSteps_m.front() += step;
        dtCurrentTrack_m = dtAllTracks_m.front();
        changeDT();
1049

1050 1051 1052
        for(; step < localTrackSteps_m.front(); ++step) {
            bends_m = 0;
            numberOfFieldEmittedParticles_m = 0;
1053

1054
            itsOpalBeamline_m.resetStatus();
1055

1056 1057
            // we dump later, after one step.
            // dumpStats(step, true, true);
1058 1059


1060 1061
            timeIntegration1(pusher);
            timeIntegration1_bgf(pusher);
1062

1063
            itsBunch->calcBeamParameters();