ParallelTTracker.cpp 121 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 62 63

#include "AbstractObjects/OpalData.h"

#include "BasicActions/Option.h"

#include "Distribution/Distribution.h"
64
#include "ValueDefinitions/RealVariable.h"
gsell's avatar
gsell committed
65 66
#include "Utilities/Timer.h"
#include "Utilities/OpalException.h"
67
#include "Solvers/SurfacePhysicsHandler.hh"
gsell's avatar
gsell committed
68 69 70 71 72 73 74 75 76 77
#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
78 79
                                   bool revTrack,
				   size_t N):
80 81 82 83 84 85 86 87 88 89 90 91
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),
92 93
//--------- Added by Xiaoying Pang 04/22/2014 ---------------
wakeFunction_m(NULL),
94 95 96 97 98 99 100 101 102 103 104
surfaceStatus_m(false),
secondaryFlg_m(false),
mpacflg_m(true),
nEmissionMode_m(false),
zStop_m(0.0),
scaleFactor_m(1.0),
vscaleFactor_m(scaleFactor_m),
recpGamma_m(1.0),
rescale_coeff_m(1.0),
dtTrack_m(0.0),
surfaceEmissionStop_m(-1),
105
specifiedNPart_m(N),
106 107 108 109 110 111 112 113 114 115 116 117
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()),
localTrackSteps_m(0),
maxNparts_m(0),
numberOfFieldEmittedParticles_m(std::numeric_limits<size_t>::max()),
bends_m(0),
numParticlesInSimulation_m(0),
kraus's avatar
kraus committed
118
space_orientation_m(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0),
119 120 121 122 123 124
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),
125 126
SeyNum_m(0.0),
sphys_m(NULL){
gsell's avatar
gsell committed
127 128 129 130 131 132 133 134 135
}

ParallelTTracker::ParallelTTracker(const Beamline &beamline,
                                   PartBunch &bunch,
                                   DataSink &ds,
                                   const PartData &reference,
                                   bool revBeam,
                                   bool revTrack,
                                   int maxSTEPS,
136
                                   double zstop,
adelmann's avatar
adelmann committed
137 138
                                   int timeIntegrator,
				   size_t N):
139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
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),
zStop_m(zstop),
scaleFactor_m(itsBunch->getdT() * Physics::c),
vscaleFactor_m(scaleFactor_m),
recpGamma_m(1.0),
rescale_coeff_m(1.0),
dtTrack_m(0.0),
surfaceEmissionStop_m(-1),
162
specifiedNPart_m(N),
163 164 165 166 167 168 169 170 171 172 173 174
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()),
localTrackSteps_m(maxSTEPS),
maxNparts_m(0),
numberOfFieldEmittedParticles_m(numeric_limits<size_t>::max()),
bends_m(0),
numParticlesInSimulation_m(0),
kraus's avatar
kraus committed
175
space_orientation_m(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0),
176 177 178 179 180 181 182 183
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) {
184

185
    //    itsBeamline = dynamic_cast<Beamline*>(beamline.clone());
186

gsell's avatar
gsell committed
187 188
}

189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
#ifdef HAVE_AMR_SOLVER
ParallelTTracker::ParallelTTracker(const Beamline &beamline,
                                   PartBunch &bunch,
                                   DataSink &ds,
                                   const PartData &reference,
                                   bool revBeam,
                                   bool revTrack,
                                   int maxSTEPS,
                                   double zstop,
                                   int timeIntegrator,
				   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),
zStop_m(zstop),
scaleFactor_m(itsBunch->getdT() * Physics::c),
vscaleFactor_m(scaleFactor_m),
recpGamma_m(1.0),
rescale_coeff_m(1.0),
dtTrack_m(0.0),
surfaceEmissionStop_m(-1),
224
specifiedNPart_m(N),
225 226 227 228 229 230 231 232 233 234 235 236
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()),
localTrackSteps_m(maxSTEPS),
maxNparts_m(0),
numberOfFieldEmittedParticles_m(numeric_limits<size_t>::max()),
bends_m(0),
numParticlesInSimulation_m(0),
kraus's avatar
kraus committed
237
space_orientation_m(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0),
238 239 240 241 242 243 244 245 246 247 248
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){
}
#endif
gsell's avatar
gsell committed
249 250

ParallelTTracker::~ParallelTTracker() {
251

gsell's avatar
gsell committed
252 253 254
}

void ParallelTTracker::applySchottkyCorrection(PartBunch &itsBunch, int ne, double t, double rescale_coeff) {
255

256
    Inform msg("ParallelTTracker ");
gsell's avatar
gsell committed
257 258
    const long ls = 0;
    /*
259 260
     Now I can calculate E_{rf} at each position
     of the newely generated particles and rescale Q
261

262 263 264 265
     Note:
     For now I only sample the field of the last emitted particles.
     Space charge is not yet included
     */
266 267


gsell's avatar
gsell committed
268 269
    double laser_erg = itsBunch.getLaserEnergy(); // 4.7322; energy of single photon of 262nm laser  [eV]
    double workFunction = itsBunch.getWorkFunctionRf(); // espace energy for copper (4.31)  [eV]
270

gsell's avatar
gsell committed
271
    const double schottky_coeff = 0.037947; // coeffecient for calculate schottky potenial from E field [eV/(MV^0.5)]
272

gsell's avatar
gsell committed
273 274
    if(ne == 0)
        return ;
275

gsell's avatar
gsell committed
276 277 278 279 280 281 282
    double Ez = 0;
    double obtain_erg = 0;
    double par_t = 0;
    for(int k = 0; k < ne; k++) {
        size_t n = itsBunch.getLocalNum() - k - 1;
        Vector_t externalE(0.0);
        Vector_t externalB(0.0);
283

gsell's avatar
gsell committed
284 285 286 287
        itsBunch.R[n] *= Vector_t(Physics::c * itsBunch.dt[n]);
        par_t = t + itsBunch.dt[n] / 2;
        itsOpalBeamline_m.getFieldAt(n, itsBunch.R[n], ls, par_t, externalE, externalB);
        Ez = externalE(2);
288

gsell's avatar
gsell committed
289 290 291 292 293 294
        // fabs(Ez): if the field of cathode surface is in the right direction, it will increase the
        // energy which electron obtain. If the field is in the wrong direction, this particle will
        // be back to the cathode surface and then be deleted automaticly by OPAL,  we don't add
        // another logical branch to handle this. So fabs is the simplest way to handle this
        obtain_erg = laser_erg - workFunction + schottky_coeff * sqrt(fabs(Ez) / 1E6);
        double schottkyScale = obtain_erg * obtain_erg * rescale_coeff;
295

gsell's avatar
gsell committed
296 297 298 299 300 301
        itsBunch.Q[n] *= schottkyScale;
        itsBunch.R[n] /= Vector_t(Physics::c * itsBunch.dt[n]);
    }
}

double ParallelTTracker::schottkyLoop(double rescale_coeff) {
302

303
    Inform msg("ParallelTTracker ");
304

gsell's avatar
gsell committed
305 306 307 308
    double recpgamma;
    double t = 0.0;
    double dt = itsBunch->getdT();
    Vector_t vscaleFactor = Vector_t(scaleFactor_m);
309

gsell's avatar
gsell committed
310 311
    unsigned long long step = 0;
    unsigned int emissionSteps = 0;
312

gsell's avatar
gsell committed
313 314 315 316
    Vector_t um, a, s;
    Vector_t externalE, externalB;
    BorisPusher pusher(itsReference);
    Vector_t rmin, rmax;
317

gsell's avatar
gsell committed
318
    bool global_EOL;
319

gsell's avatar
gsell committed
320 321
    bool hasSwitchedToTEmission = false;
    bool hasSwitchedBackToTTrack = false;
322

gsell's avatar
gsell committed
323
    size_t totalParticles_i = itsBunch->getTotalNum();
324

325 326 327
    msg << "*****************************************************************" << endl;
    msg << " Estimate Schottky correction                                    " << endl;
    msg << "*****************************************************************" << endl;
328

gsell's avatar
gsell committed
329 330 331 332 333 334 335 336
    double margin = 0.0;
    if(!mpacflg_m) {
        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));
        }
337 338

        if(!(itsBunch->WeHaveEnergyBins())) {
gsell's avatar
gsell committed
339 340 341 342 343
            IpplTimings::startTimer(BinRepartTimer_m);
            itsBunch->do_binaryRepart();
            IpplTimings::stopTimer(BinRepartTimer_m);
            Ippl::Comm->barrier();
        }
344

gsell's avatar
gsell committed
345 346 347 348 349 350 351 352 353
        // Check if there are any particles in simulation. If there are,
        // as in a restart, use the usual function to calculate beam
        // parameters. If not, calculate beam parameters of the initial
        // beam distribution.
        if(totalParticles_i == 0) {// fixme: maybe cause nonsense output if initialized momenta=0; Q: by Chuan.
            itsBunch->calcBeamParametersInitial();
        } else {
            itsBunch->calcBeamParameters();
        }
354

gsell's avatar
gsell committed
355 356
        RefPartR_suv_m = RefPartR_zxy_m = itsBunch->get_rmean();
        RefPartP_suv_m = RefPartP_zxy_m = itsBunch->get_pmean();
357

gsell's avatar
gsell committed
358 359 360
        if(!OpalData::getInstance()->hasBunchAllocated()) {
            updateSpaceOrientation(false);  // vec{p} = (0,0,p_z), vec{r} = (0,0,z)
        }
361

gsell's avatar
gsell committed
362 363 364 365 366
        RefPartR_suv_m = itsBunch->get_rmean();
        RefPartP_suv_m = itsBunch->get_pmean();
        /* Activate all elements which influence the particles when the simulation starts;
         *  mark all elements which are already past.
         */
367

gsell's avatar
gsell committed
368
        /*
369 370 371
         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.
         */
gsell's avatar
gsell committed
372
        itsBunch->get_bounds(rmin, rmax);
373
        margin = 10. * RefPartP_suv_m(2) * scaleFactor_m / sqrt(1.0 + pSqr(RefPartP_suv_m, RefPartP_suv_m));
gsell's avatar
gsell committed
374 375 376
        margin = 0.01 > margin ? 0.01 : margin;
        itsOpalBeamline_m.switchElements(rmin(2) - margin, rmax(2) + margin);
    }
377

gsell's avatar
gsell committed
378 379 380 381 382 383
    double minBinEmitted  = 10.0;
    RealVariable *ar = dynamic_cast<RealVariable *>(OpalData::getInstance()->find("MINBINEMITTED"));
    if(ar) {
        minBinEmitted = ar->getReal();  // the space charge solver crashes if we use less than ~10 particles.
        // This variable controls the number of particles to be emitted before we use
        // the space charge solver.
384
        msg << "MINBINEMITTED " << minBinEmitted << endl;
gsell's avatar
gsell committed
385
    }
386 387


gsell's avatar
gsell committed
388 389 390 391 392
    double minStepforReBin  = 10000.0;
    RealVariable *br = dynamic_cast<RealVariable *>(OpalData::getInstance()->find("MINSTEPFORREBIN"));
    if(br) {
        minStepforReBin = br->getReal();  // this variable controls the minimal number of steps of emission (using bins)
        // before we can merge the bins
393
        msg << "MINSTEPFORREBIN " << minStepforReBin << endl;
gsell's avatar
gsell committed
394
    }
395

gsell's avatar
gsell committed
396 397 398 399
    int repartFreq = 1000;
    RealVariable *rep = dynamic_cast<RealVariable *>(OpalData::getInstance()->find("REPARTFREQ"));
    if(rep) {
        repartFreq = static_cast<int>(rep->getReal());  // this variable controls the minimal number of steps until we repartition the particles
400
        msg << "REPARTFREQ " << repartFreq << endl;
gsell's avatar
gsell committed
401
    }
402

gsell's avatar
gsell committed
403 404 405
    // there is no point to do repartitioning with one node
    if(Ippl::getNodes() == 1)
        repartFreq = 1000000;
406

gsell's avatar
gsell committed
407
    size_t totalParticles_f = 0;
408

409
    for(; step < localTrackSteps_m; ++step) {
gsell's avatar
gsell committed
410
        global_EOL = true;  // check if any particle hasn't reached the end of the field from the last element
411

gsell's avatar
gsell committed
412
        itsOpalBeamline_m.resetStatus();
413

gsell's avatar
gsell committed
414
        IpplTimings::startTimer(timeIntegrationTimer1_m);
415

gsell's avatar
gsell committed
416
        // reset E and B to Vector_t(0.0) for every step
417

gsell's avatar
gsell committed
418 419
        itsBunch->Ef = Vector_t(0.0);
        itsBunch->Bf = Vector_t(0.0);
420

gsell's avatar
gsell committed
421 422
        Nimpact_m = 0; // Initial parallel plate benchmark variable.
        SeyNum_m = 0; // Initial parallel plate benchmark variable.
423

gsell's avatar
gsell committed
424 425 426 427 428 429 430 431 432 433
        for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) {
            //scale each particle with c*dt
            itsBunch->R[i] /= vscaleFactor;
            pusher.push(itsBunch->R[i], itsBunch->P[i], itsBunch->dt[i]);
            // update local coordinate system of particleInform &PartBunc
            itsBunch->X[i] /= vscaleFactor;
            pusher.push(itsBunch->X[i], TransformTo(itsBunch->P[i], itsOpalBeamline_m.getOrientation(itsBunch->LastSection[i])),
                        itsBunch->getdT());
            itsBunch->X[i] *= vscaleFactor;
        }
434

gsell's avatar
gsell committed
435 436 437
        if(totalParticles_i > minBinEmitted) {
            itsBunch->boundp();
        }
438

gsell's avatar
gsell committed
439
        IpplTimings::stopTimer(timeIntegrationTimer1_m);
440

gsell's avatar
gsell committed
441
        itsBunch->calcBeamParameters();
442 443


gsell's avatar
gsell committed
444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459
        /** \f[ Space Charge  \f]
         */
        if(itsBunch->hasFieldSolver() && totalParticles_i > minBinEmitted && fabs(itsBunch->getChargePerParticle()) > 0.0) {
            // Do repartition if we have enough particles.
            if(totalParticles_i > 1000 && (((step + 1) % repartFreq) == 0)) {
                INFOMSG("*****************************************************************" << endl);
                INFOMSG("do repartition because of repartFreq" << endl);
                INFOMSG("*****************************************************************" << endl);
                IpplTimings::startTimer(BinRepartTimer_m);
                itsBunch->do_binaryRepart();
                IpplTimings::stopTimer(BinRepartTimer_m);
                Ippl::Comm->barrier();
                INFOMSG("*****************************************************************" << endl);
                INFOMSG("do repartition done" << endl);
                INFOMSG("*****************************************************************" << endl);
            }
460

gsell's avatar
gsell committed
461
            // Calculate space charge.
462
            if(itsBunch->WeHaveEnergyBins()) {
gsell's avatar
gsell committed
463 464 465
                // When we have energy bins.
                itsBunch->calcGammas();
                ParticleAttrib<double> Q_back = itsBunch->Q;
466
                itsBunch->resetInterpolationCache();
gsell's avatar
gsell committed
467 468 469 470 471 472 473 474 475
                for(int binNumber = 0; binNumber <= itsBunch->getLastemittedBin() && binNumber < itsBunch->getNumBins(); ++binNumber) {
                    itsBunch->setBinCharge(binNumber);
                    itsBunch->computeSelfFields(binNumber);
                    itsBunch->Q = Q_back;
                }
            } else {
                // When we don't.
                itsBunch->computeSelfFields();
                /**
476 477 478 479
                 Need this maybe for the adaptive time integration scheme
                 pair<Vector_t,Vector_t> eExtrema = itsBunch->getEExtrema();
                 INFOMSG("maxE= " << eExtrema.first << " minE= " << eExtrema.second << endl);
                 */
gsell's avatar
gsell committed
480 481
            }
        }
482

gsell's avatar
gsell committed
483
        IpplTimings::startTimer(timeIntegrationTimer2_m);
484 485


gsell's avatar
gsell committed
486
        /*
487 488 489 490
         transport and emit particles
         that passed the cathode in the first
         half-step or that would pass it in the
         second half-step.
491

492 493
         to make IPPL and the field solver happy
         make sure that at least 10 particles are emitted
494

495 496
         also remember that node 0 has
         all the particles to be emitted
497

498 499 500
         this has to be done *after* the calculation of the
         space charges! thereby we neglect space charge effects
         in the very first step of a new-born particle.
501

502
         */
503 504 505

        if((itsBunch->WeHaveEnergyBins())) {

gsell's avatar
gsell committed
506 507
            // switch to TEmission
            if(!hasSwitchedToTEmission) {
508
                dt = itsBunch->GetEmissionDeltaT();
gsell's avatar
gsell committed
509 510 511
                itsBunch->setdT(dt);
                scaleFactor_m = dt * Physics::c;
                vscaleFactor = Vector_t(scaleFactor_m);
512
                msg << "Changing emission time step to: " << dt << endl;
gsell's avatar
gsell committed
513 514
                hasSwitchedToTEmission = true;
            }
515

gsell's avatar
gsell committed
516
            int ne = 0;
517 518 519 520 521 522 523 524 525
            Vector_t externalE = Vector_t(0.0);
            Vector_t externalB = Vector_t(0.0);
            itsOpalBeamline_m.getFieldAt(Vector_t(0.0),
                                         Vector_t(0.0),
                                         itsBunch->getT() + 0.5 * itsBunch->getdT(),
                                         externalE,
                                         externalB);
            ne += itsBunch->EmitParticles(externalE[2]);

gsell's avatar
gsell committed
526 527
            if(Options::schottkyCorrection && !hasSwitchedBackToTTrack)
                applySchottkyCorrection(*itsBunch, ne, t, rescale_coeff);
528

gsell's avatar
gsell committed
529 530
            reduce(ne, ne, OpAddAssign());
            totalParticles_i += ne;
531

gsell's avatar
gsell committed
532 533 534 535 536 537 538
            //emission has finished, reset to TTrack
            if(itsBunch->getNumBins() == itsBunch->getLastemittedBin() &&
               !hasSwitchedBackToTTrack) {
                //dt = dtTrack;
                itsBunch->setdT(dt);
                scaleFactor_m = dt * Physics::c;
                vscaleFactor = Vector_t(scaleFactor_m);
539
                msg << "Emission done. Switching back to track timestep: " << dt << endl;
gsell's avatar
gsell committed
540 541 542
                hasSwitchedBackToTTrack = true;
                break;
            }
543

gsell's avatar
gsell committed
544 545 546 547 548 549 550
        } else {
            //emission has finished, reset to TTrack
            if(!hasSwitchedBackToTTrack) {
                //dt = dtTrack;
                itsBunch->setdT(dt);
                scaleFactor_m = dt * Physics::c;
                vscaleFactor = Vector_t(scaleFactor_m);
551
                msg << "Emission done. Switching back to track timestep: " << dt << endl;
gsell's avatar
gsell committed
552 553
                hasSwitchedBackToTTrack = true;
            }
554

gsell's avatar
gsell committed
555
        }
556

gsell's avatar
gsell committed
557 558 559
        // push the reference particle by a half step
        recpgamma = 1.0 / sqrt(1.0 + dot(RefPartP_suv_m, RefPartP_suv_m));
        RefPartR_zxy_m += RefPartP_zxy_m * recpgamma / 2. * scaleFactor_m;
560

gsell's avatar
gsell committed
561 562 563 564 565 566 567
        //
        // get external fields for all particles
        //
        IpplTimings::startTimer(timeFieldEvaluation_m);
        for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) {
            //FIXME: rethink scaling!
            itsBunch->R[i] *= Vector_t(Physics::c * itsBunch->dt[i], Physics::c * itsBunch->dt[i], Physics::c * itsBunch->dt[i]);
568

gsell's avatar
gsell committed
569 570 571 572 573 574 575 576 577
            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;
            }
            const unsigned long rtv = itsOpalBeamline_m.getFieldAt(i, itsBunch->R[i], ls, t + itsBunch->dt[i] / 2., externalE, externalB);
578

gsell's avatar
gsell committed
579
            global_EOL = global_EOL && (rtv & BEAMLINE_EOL);
580

gsell's avatar
gsell committed
581 582 583 584 585
            // skip rest of the particle push if the
            // particle is out of bounds i.e. does not see
            // a E or B field
            if(rtv & BEAMLINE_OOB)
                itsBunch->Bin[i] = -1;
586 587


gsell's avatar
gsell committed
588 589
            itsBunch->Ef[i] += externalE;
            itsBunch->Bf[i] += externalB;
590

gsell's avatar
gsell committed
591
            itsBunch->R[i] /= Vector_t(Physics::c * itsBunch->dt[i], Physics::c * itsBunch->dt[i], Physics::c * itsBunch->dt[i]);
592

gsell's avatar
gsell committed
593
            // in case a particle is pushed behind the emission surface, delete the particle
594

gsell's avatar
gsell committed
595 596
            if(itsBunch->R[i](2) < 0)
                itsBunch->Bin[i] = -1;
597

gsell's avatar
gsell committed
598
        }
599

gsell's avatar
gsell committed
600
        IpplTimings::stopTimer(timeFieldEvaluation_m);
601

gsell's avatar
gsell committed
602 603
        //        if(itsBunch->getLocalNum() == 0)
        //    global_EOL = false;
604

gsell's avatar
gsell committed
605
        /**
606 607
         Delete marked particles.
         */
608

gsell's avatar
gsell committed
609 610 611 612 613
        bool globPartOutOfBounds = (min(itsBunch->Bin) < 0);
        size_t ne = 0;
        if(globPartOutOfBounds) {
            ne = itsBunch->boundp_destroyT();
        }
614

gsell's avatar
gsell committed
615 616
        totalParticles_f = totalParticles_i - ne;
        if(ne > 0)
617
            msg << "* Deleted in Shotky " << ne << " particles, remaining " << totalParticles_f << " particles" << endl; //benchmark output
618

gsell's avatar
gsell committed
619
        kickParticles(pusher);
620

gsell's avatar
gsell committed
621 622 623 624
        if(totalParticles_f > 0) {
            // none of the particles is in a bending element
            updateReferenceParticle();
        }
625

gsell's avatar
gsell committed
626 627
        itsBunch->RefPart_R = RefPartR_zxy_m;
        itsBunch->RefPart_P = RefPartP_zxy_m;
628

gsell's avatar
gsell committed
629 630 631
        // calculate the dimensions of the bunch and add a small margin to them; then decide which elements have to be triggered
        // when an element is triggered memory is allocated and the field map is read in
        itsBunch->get_bounds(rmin, rmax);
632

gsell's avatar
gsell committed
633 634 635 636
        // trigger the elements
        margin = 3. * RefPartP_suv_m(2) * recpgamma;
        margin = 0.01 > margin ? 0.01 : margin;
        itsOpalBeamline_m.switchElements((rmin(2) - margin)*scaleFactor_m, (rmax(2) + margin)*scaleFactor_m);
637

gsell's avatar
gsell committed
638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655
        // start normal particle loop part 2 for simulation without boundary geometry.
        for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) {
            /** \f[ \vec{x}_{n+1} = \vec{x}_{n+1/2} + \frac{1}{2}\vec{v}_{n+1/2}\quad (= \vec{x}_{n+1/2} + \frac{\Delta t}{2} \frac{\vec{\beta}_{n+1/2}\gamma_{n+1/2}}{\gamma_{n+1/2}}) \f]
             * \code
             * R[i] += 0.5 * P[i] * recpgamma;
             * \endcode
             */
            pusher.push(itsBunch->R[i], itsBunch->P[i], itsBunch->dt[i]);
            //and scale back to dimensions
            itsBunch->R[i] *= Vector_t(Physics::c * itsBunch->dt[i], Physics::c * itsBunch->dt[i], Physics::c * itsBunch->dt[i]);
            // update local coordinate system
            itsBunch->X[i] /= vscaleFactor;
            pusher.push(itsBunch->X[i], TransformTo(itsBunch->P[i], itsOpalBeamline_m.getOrientation(itsBunch->LastSection[i])), itsBunch->getdT());
            itsBunch->X[i] *= vscaleFactor;
            //reset time step if particle was emitted in the first half-step
            //the particle is now in sync with the simulation timestep
            itsBunch->dt[i] = itsBunch->getdT();
        }
656

gsell's avatar
gsell committed
657
        IpplTimings::stopTimer(timeIntegrationTimer2_m);
658

gsell's avatar
gsell committed
659 660
        if(totalParticles_f > minBinEmitted)
            itsBunch->boundp();
661

gsell's avatar
gsell committed
662
        totalParticles_i = itsBunch->getTotalNum();
663 664


gsell's avatar
gsell committed
665
        t += itsBunch->getdT(); //t after a full global timestep with dT "synchronization point" for simulation time
666

gsell's avatar
gsell committed
667
        itsBunch->setT(t);
668

gsell's avatar
gsell committed
669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696
        //IFF: cheap step dump regulation
        OPALTimer::Timer myt2;
        double sposRef = 0.0;
        if(totalParticles_f > 0) {
            sposRef = itsBunch->get_sPos();
            if(totalParticles_f <= minBinEmitted) {
                INFOMSG(myt2.time() << " Step " << step << "; only " << totalParticles_f << " particles emitted; t= " << t
                        << " [s] E=" << itsBunch->get_meanEnergy() << " [MeV] " << endl);
            } else if(std::isnan(sposRef) || std::isinf(sposRef)) {
                INFOMSG(myt2.time() << " Step " << step << "; there seems to be something wrong with the position of the bunch!" << endl);
            } else {
                INFOMSG(myt2.time() << " Step " << step << " at " << sposRef << " [m] t= " << t << " [s] E=" << itsBunch->get_meanEnergy() << " [MeV] " << endl);
                if(step % Options::psDumpFreq == 0 || step % Options::statDumpFreq == 0) {
                    size_t nLoc = itsBunch->getLocalNum();
                    reduce(nLoc, nLoc, OpMultipplyAssign());
                    if((nLoc == 0) || ((step + 1) % repartFreq == 0)) {
                        INFOMSG("*****************************************************************" << endl);
                        INFOMSG("do repartition because of zero particles or repartition frequency" << endl);
                        IpplTimings::startTimer(BinRepartTimer_m);
                        itsBunch->do_binaryRepart();
                        IpplTimings::stopTimer(BinRepartTimer_m);
                        Ippl::Comm->barrier();
                        INFOMSG("done" << endl);
                        INFOMSG("*****************************************************************" << endl);
                    }
                }
            }
            /**
697 698
             Stop simulation if beyond zStop_m
             */
699
            if(sposRef > zStop_m) {
700
                localTrackSteps_m = step;
gsell's avatar
gsell committed
701 702 703 704
            }
        } else {
            INFOMSG("Step " << step << " no emission yet "  << " t= " << t << " [s]" << endl);
        }
705

gsell's avatar
gsell committed
706 707 708 709 710 711 712 713 714 715 716 717
        if(step > emissionSteps) {
            reduce(&global_EOL, &global_EOL + 1, &global_EOL, OpBitwiseAndAssign());
            if(global_EOL) {
                break;
            }
        }
        // this seams to fix Ticket #12
        //  Ippl::Comm->barrier();
        itsBunch->get_bounds(rmin, rmax);
        // trigger the elements
        RefPartP_suv_m = itsBunch->get_pmean();
        recpgamma = 1. / sqrt(1.0 + dot(RefPartP_suv_m, RefPartP_suv_m));
718

gsell's avatar
gsell committed
719 720 721 722 723
        margin = 10. * RefPartP_suv_m(2) * recpgamma * scaleFactor_m;
        margin = 0.01 > margin ? 0.01 : margin;
    }
    OPALTimer::Timer myt3;
    OpalData::getInstance()->setLastStep(step);
724
    msg << "done executing Schottky loop " << myt3.time() << endl;
gsell's avatar
gsell committed
725 726 727 728
    return itsBunch->getCharge();
}


729
void ParallelTTracker::checkCavity(double s, Component *& comp, double & cavity_start_pos) {
730

731 732 733 734
    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())) {
735
            comp = (*fit).getElement().get();
736 737 738
            cavity_start_pos = (*fit).getStart();
            currently_ap_cavity_m = fit;
            return;
gsell's avatar
gsell committed
739
        }
740
    }
gsell's avatar
gsell committed
741 742
}

743

744
void ParallelTTracker::doOneStep(BorisPusher & pusher) {
745

746
    bends_m = 0;
747

gsell's avatar
gsell committed
748
    double t = itsBunch->getT();
749

750 751
    // 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.
752

753
    switchElements(10.0);
754

gsell's avatar
gsell committed
755
    itsOpalBeamline_m.resetStatus();
756

757
    timeIntegration1(pusher);
758

759
    itsBunch->calcBeamParameters();
760

761
    // reset E and B to Vector_t(0.0) for every step
gsell's avatar
gsell committed
762 763
    itsBunch->Ef = Vector_t(0.0);
    itsBunch->Bf = Vector_t(0.0);
764

765
    computeExternalFields();
766

767
    timeIntegration2(pusher);
768

769
    t += itsBunch->getdT();
gsell's avatar
gsell committed
770 771 772 773 774
    itsBunch->setT(t);
}


void ParallelTTracker::applyEntranceFringe(double angle, double curve,
775
                                           const BMultipoleField &field, double scale) {
gsell's avatar
gsell committed
776 777 778 779 780 781 782 783
}


void ParallelTTracker::applyExitFringe(double angle, double curve,
                                       const BMultipoleField &field, double scale) {
}


784
void ParallelTTracker::showCavities(Inform &msg) {
785

786
    msg << "Found the following cavities:" << endl;
787

gsell's avatar
gsell committed
788
    for(FieldList::iterator fit = cavities_m.begin(); fit != cavities_m.end(); ++ fit) {
789
        msg << (*fit).getElement()->getName()
790 791
        << " from " << (*fit).getStart() << " to "
        << (*fit).getEnd() << " (m) phi=";
gsell's avatar
gsell committed
792
        if((*fit).getElement()->getType() == "TravelingWave")
793
            msg << static_cast<TravelingWave *>((*fit).getElement().get())->getPhasem() / Physics::pi * 180.0 << endl;
gsell's avatar
gsell committed
794
        else
795
            msg << static_cast<RFCavity *>((*fit).getElement().get())->getPhasem() / Physics::pi * 180.0 << endl;
gsell's avatar
gsell committed
796
    }
797
    msg << endl << endl;
gsell's avatar
gsell committed
798 799 800
}


801
void ParallelTTracker::updateRFElement(std::string elName, double maxPhi) {
gsell's avatar
gsell committed
802
    /**
803 804 805 806
     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.
     */
gsell's avatar
gsell committed
807 808 809 810
    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") {
811
                phi  =  static_cast<TravelingWave *>((*fit).getElement().get())->getPhasem();
gsell's avatar
gsell committed
812
                phi += maxPhi;
813
                static_cast<TravelingWave *>((*fit).getElement().get())->updatePhasem(phi);
gsell's avatar
gsell committed
814
            } else {
815
                phi  = static_cast<RFCavity *>((*fit).getElement().get())->getPhasem();
gsell's avatar
gsell committed
816
                phi += maxPhi;
817
                static_cast<RFCavity *>((*fit).getElement().get())->updatePhasem(phi);
gsell's avatar
gsell committed
818 819 820 821 822 823 824 825
            }
        }
    }
}


void ParallelTTracker::updateAllRFElements(double phiShift) {
    /**
826 827 828
     All RF-Elements gets updated, where the phiShift is the
     global phase shift in units of seconds.
     */
829
    Inform msg("ParallelTTracker ");
gsell's avatar
gsell committed
830 831 832 833
    double phi = 0;
    double freq = 0.0;
    const double RADDEG = 1.0 / Physics::pi * 180.0;
    //   Inform m ("updateALLRFElements ",INFORM_ALL_NODES);
834
    msg << "\n-------------------------------------------------------------------------------------\n";
gsell's avatar
gsell committed
835 836
    for(FieldList::iterator fit = cavities_m.begin(); fit != cavities_m.end(); ++ fit) {
        if(fit != cavities_m.begin())
837
            msg << "\n";
gsell's avatar
gsell committed
838
        if((*fit).getElement()->getType() == "TravelingWave") {
839 840
            freq = static_cast<TravelingWave *>((*fit).getElement().get())->getFrequencym();
            phi = static_cast<TravelingWave *>((*fit).getElement().get())->getPhasem();
841
            msg << (*fit).getElement()->getName()
842 843
            << ": phi= phi_nom + phi_maxE + global phase shift= " << (phi*RADDEG)-(phiShift*freq*RADDEG) << " degree, "
            << "(global phase shift= " << -phiShift *freq *RADDEG << " degree)\n";
gsell's avatar
gsell committed
844
            phi -= (phiShift * freq);
845
            static_cast<TravelingWave *>((*fit).getElement().get())->updatePhasem(phi);
gsell's avatar
gsell committed
846
        } else {
847 848
            freq = static_cast<RFCavity *>((*fit).getElement().get())->getFrequencym();
            phi = static_cast<RFCavity *>((*fit).getElement().get())->getPhasem();
849
            msg << (*fit).getElement()->getName()
850 851
            << ": phi= phi_nom + phi_maxE + global phase shift= " << (phi*RADDEG)-(phiShift*freq*RADDEG) << " degree, "
            << "global phase shift= " << -phiShift *freq *RADDEG << " degree\n";
gsell's avatar
gsell committed
852
            phi -= (phiShift * freq);
853
            static_cast<RFCavity *>((*fit).getElement().get())->updatePhasem(phi);
gsell's avatar
gsell committed
854 855
        }
    }
856
    msg << "-------------------------------------------------------------------------------------\n"
857
	<< endl;
gsell's avatar
gsell committed
858 859 860 861 862
}


FieldList ParallelTTracker::executeAutoPhaseForSliceTracker() {
    Inform msg("executeAutoPhaseForSliceTracker ");
863

gsell's avatar
gsell committed
864
    double gPhaseSave;
865

gsell's avatar
gsell committed
866 867
    gPhaseSave = OpalData::getInstance()->getGlobalPhaseShift();
    OpalData::getInstance()->setGlobalPhaseShift(0.0);
868

gsell's avatar
gsell committed
869 870 871 872 873 874
    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);
875
        if(zbegin < zStop_m && zend >= zStop_m) {
876
            msg << "\033[0;31m"
877 878 879 880 881 882
            << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\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;
883
            static_cast<Monitor *>(it->getElement().get())->moveBy(-zend - 0.001);
gsell's avatar
gsell committed
884 885 886 887
            itsOpalBeamline_m.removeElement(it->getElement()->getName());
        }
    }
    itsOpalBeamline_m.prepareSections();
888

gsell's avatar
gsell committed
889
    cavities_m = itsOpalBeamline_m.getElementByType("RFCavity");
890
    currently_ap_cavity_m = cavities_m.end();
891 892
    FieldList travelingwaves = itsOpalBeamline_m.getElementByType("TravelingWave");
    cavities_m.merge(travelingwaves, OpalField::SortAsc);
893

gsell's avatar
gsell committed
894 895 896 897 898 899 900 901 902 903
    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;
904

gsell's avatar
gsell committed
905
    double zStop = itsOpalBeamline_m.calcBeamlineLenght();
906

gsell's avatar
gsell committed
907
    msg << "Preparation done zstop= " << zStop << endl;
908 909

    if(Ippl::myNode() == 0)
910 911 912
      itsBunch->create(1);
    itsBunch->update();

gsell's avatar
gsell committed
913
    if(Ippl::myNode() == 0) {
914 915 916 917 918 919
      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;
920

921 922
      RefPartP_suv_m = iniP;
      RefPartR_suv_m = iniR;
923 924 925 926 927
    }
    updateSpaceOrientation(false);
    executeAutoPhase(Options::autoPhase, zStop);
    if(Ippl::myNode() == 0) {
      RefPartP_suv_m = itsBunch->P[0];
928 929 930 931
      // need to rebuild for updateAllRFElements
      cavities_m = itsOpalBeamline_m.getElementByType("RFCavity");
      travelingwaves = itsOpalBeamline_m.getElementByType("TravelingWave");
      cavities_m.merge(travelingwaves, OpalField::SortAsc);
932

933 934 935 936
      // 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());
937

938
      for(std::vector<MaxPhasesT>::iterator it = OpalData::getInstance()->getFirstMaxPhases(); it < OpalData::getInstance()->getLastMaxPhases(); it++) {
939 940 941 942
	putMessage(*mess, (*it).first);
	putMessage(*mess, (*it).second);
      }
      Ippl::Comm->broadcast_all(mess, tag);
943 944

      delete mess;
gsell's avatar
gsell committed
945
    } else {
946 947 948 949 950
      // 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++) {
951
	std::string elName;
952 953 954 955 956 957
	double maxPhi;
	getMessage(*mess, elName);
	getMessage(*mess, maxPhi);
	updateRFElement(elName, maxPhi);
	OpalData::getInstance()->setMaxPhase(elName, maxPhi);
      }
958 959

      delete mess;
gsell's avatar
gsell committed
960
    }
961 962

    if(Ippl::myNode() == 0)
963
      itsBunch->destroy(1, 0);
964 965
    itsBunch->update();

gsell's avatar
gsell committed
966 967 968 969 970 971 972
    OpalData::getInstance()->setGlobalPhaseShift(gPhaseSave);
    return cavities_m;
}


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

gsell's avatar
gsell committed
974
    const double RADDEG = 180.0 / Physics::pi;
975

gsell's avatar
gsell committed
976
    Vector_t rmin, rmax;
977

978
    size_t maxStepsSave = localTrackSteps_m;
gsell's avatar
gsell committed
979
    size_t step = 0;
980

981 982
    double tSave = itsBunch->getT();
    double dTSave = itsBunch->getdT();
gsell's avatar
gsell committed
983
    double scaleFactorSave = scaleFactor_m;
984

985 986 987 988
    int dtfraction = 2;
    double newDt = itsBunch->getdT() / dtfraction;
    itsBunch->setdT(newDt);
    itsBunch->dt = newDt;
989

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

gsell's avatar
gsell committed
992
    scaleFactor_m = itsBunch->dt[0] * Physics::c;
993
    vscaleFactor_m = Vector_t(scaleFactor_m);