ParallelCyclotronTracker.cpp 125 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14
// ------------------------------------------------------------------------
// $RCSfile: ParallelCyclotronTracker.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1 $initialLocalNum_m
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: ParallelCyclotronTracker
//   The class for tracking particles with 3D space charge in Cyclotrons and FFAG's
//
// ------------------------------------------------------------------------
//
// $Date: 2007/10/17 04:00:08 $
15
// $Author: adelmann, yang, winklehner $
gsell's avatar
gsell committed
16 17
//
// ------------------------------------------------------------------------
kraus's avatar
kraus committed
18

19
#include <Algorithms/Ctunes.hpp>
kraus's avatar
kraus committed
20
#include "Algorithms/ParallelCyclotronTracker.h"
21 22
#include "Algorithms/PolynomialTimeDependence.h"
#include "Elements/OpalPolynomialTimeDependence.h"
23

gsell's avatar
gsell committed
24 25 26 27
#include <cfloat>
#include <iostream>
#include <fstream>
#include <vector>
28
#include "AbstractObjects/OpalData.h"
gsell's avatar
gsell committed
29 30 31 32

#include "AbsBeamline/Collimator.h"
#include "AbsBeamline/Corrector.h"
#include "AbsBeamline/Cyclotron.h"
adelmann's avatar
adelmann committed
33
#include "AbsBeamline/Degrader.h"
gsell's avatar
gsell committed
34 35 36 37
#include "AbsBeamline/Diagnostic.h"
#include "AbsBeamline/Drift.h"
#include "AbsBeamline/ElementBase.h"
#include "AbsBeamline/Lambertson.h"
38
#include "AbsBeamline/Offset.h"
gsell's avatar
gsell committed
39 40 41
#include "AbsBeamline/Marker.h"
#include "AbsBeamline/Monitor.h"
#include "AbsBeamline/Multipole.h"
ext-rogers_c's avatar
ext-rogers_c committed
42
#include "AbsBeamline/MultipoleT.h"
gsell's avatar
gsell committed
43 44 45 46 47
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/RFQuadrupole.h"
#include "AbsBeamline/SBend.h"
48
#include "AbsBeamline/SBend3D.h"
49
#include "AbsBeamline/ScalingFFAGMagnet.h"
gsell's avatar
gsell committed
50 51 52 53 54
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/CyclotronValley.h"
#include "AbsBeamline/Stripper.h"
55
#include "AbsBeamline/VariableRFCavity.h"
56

57 58
#include "AbstractObjects/Element.h"

59
#include "Beamlines/FlaggedBeamline.h"
60
#include "Elements/OpalBeamline.h"
kraus's avatar
kraus committed
61
#include "AbsBeamline/Ring.h"
gsell's avatar
gsell committed
62 63 64

#include "BeamlineGeometry/Euclid3D.h"
#include "BeamlineGeometry/PlanarArcGeometry.h"
Jianjun Yang's avatar
Jianjun Yang committed
65
#include "BeamlineGeometry/RBendGeometry.h"
ext-rogers_c's avatar
ext-rogers_c committed
66
#include "BeamlineGeometry/StraightGeometry.h"
gsell's avatar
gsell committed
67 68 69 70 71 72 73 74 75 76 77
#include "Beamlines/Beamline.h"

#include "Fields/BMultipoleField.h"
#include "FixedAlgebra/FTps.h"
#include "FixedAlgebra/FTpsMath.h"
#include "FixedAlgebra/FVps.h"

#include "Physics/Physics.h"

#include "Utilities/OpalException.h"

78
#include "BasicActions/DumpFields.h"
79
#include "BasicActions/DumpEMFields.h"
80

81
#include "Structure/H5PartWrapperForPC.h"
82
#include "Structure/BoundaryGeometry.h"
83
#include "Utilities/Options.h"
gsell's avatar
gsell committed
84 85 86 87 88 89 90

#include "Ctunes.h"
#include <cassert>

#include <hdf5.h>
#include "H5hut.h"

frey_m's avatar
frey_m committed
91 92
//FIXME Remove headers and dynamic_cast in readOneBunchFromFile
#include "Algorithms/PartBunch.h"
93
#ifdef ENABLE_AMR
frey_m's avatar
frey_m committed
94 95
    #include "Algorithms/AmrPartBunch.h"
#endif
frey_m's avatar
frey_m committed
96

gsell's avatar
gsell committed
97 98
class Beamline;
class PartData;
99

gsell's avatar
gsell committed
100 101 102
using Physics::pi;
using Physics::q_e;

103

104
const double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns
105

106 107 108
Vector_t const ParallelCyclotronTracker::xaxis = Vector_t(1.0, 0.0, 0.0);
Vector_t const ParallelCyclotronTracker::yaxis = Vector_t(0.0, 1.0, 0.0);
Vector_t const ParallelCyclotronTracker::zaxis = Vector_t(0.0, 0.0, 1.0);
gsell's avatar
gsell committed
109

snuverink_j's avatar
snuverink_j committed
110
//#define PSdim 6
gsell's avatar
gsell committed
111 112 113 114 115 116 117 118 119 120 121 122 123 124

extern Inform *gmsg;

// typedef FVector<double, PSdim> Vector;

/**
 * Constructor ParallelCyclotronTracker
 *
 * @param beamline
 * @param reference
 * @param revBeam
 * @param revTrack
 */
ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
125 126
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack):
gsell's avatar
gsell committed
127
    Tracker(beamline, reference, revBeam, revTrack),
128 129
    itsDataSink(nullptr),
    bgf_m(nullptr),
130
    lastDumpedStep_m(0),
131
    eta_m(0.01),
gsell's avatar
gsell committed
132 133
    myNode_m(Ippl::myNode()),
    initialLocalNum_m(0),
134
    initialTotalNum_m(0),
135
    onebunch_m(OpalData::getInstance()->getInputBasename() + "-onebunch.h5"),
136
    opalRing_m(nullptr),
137 138
    itsStepper_mp(nullptr),
    mode_m(MODE::UNDEFINED),
139
    stepper_m(stepper::INTEGRATOR::UNDEFINED) {
gsell's avatar
gsell committed
140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
}

/**
 * Constructor ParallelCyclotronTracker
 *
 * @param beamline
 * @param bunch
 * @param ds
 * @param reference
 * @param revBeam
 * @param revTrack
 * @param maxSTEPS
 * @param timeIntegrator
 */
ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
frey_m's avatar
frey_m committed
156
                                                   PartBunchBase<double, 3> *bunch,
gsell's avatar
gsell committed
157 158 159
                                                   DataSink &ds,
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack,
160
                                                   int maxSTEPS, int timeIntegrator):
161
    Tracker(beamline, bunch, reference, revBeam, revTrack),
162
    bgf_m(nullptr),
gsell's avatar
gsell committed
163
    maxSteps_m(maxSTEPS),
164
    lastDumpedStep_m(0),
165
    eta_m(0.01),
gsell's avatar
gsell committed
166
    myNode_m(Ippl::myNode()),
frey_m's avatar
frey_m committed
167 168
    initialLocalNum_m(bunch->getLocalNum()),
    initialTotalNum_m(bunch->getTotalNum()),
169
    onebunch_m(OpalData::getInstance()->getInputBasename() + "-onebunch.h5"),
170 171
    opalRing_m(nullptr),
    itsStepper_mp(nullptr) {
gsell's avatar
gsell committed
172 173
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;
frey_m's avatar
frey_m committed
174
    //  scaleFactor_m = itsBunch_m->getdT() * c;
gsell's avatar
gsell committed
175
    scaleFactor_m = 1;
176
    multiBunchMode_m = MB_MODE::NONE;
gsell's avatar
gsell committed
177 178 179 180 181

    IntegrationTimer_m = IpplTimings::getTimer("Integration");
    TransformTimer_m   = IpplTimings::getTimer("Frametransform");
    DumpTimer_m        = IpplTimings::getTimer("Dump");
    BinRepartTimer_m   = IpplTimings::getTimer("Binaryrepart");
182

183 184 185 186 187 188 189 190 191
    // FIXME Change track command
    if ( initialTotalNum_m == 1 ) {
        mode_m = MODE::SINGLE;
    } else if ( initialTotalNum_m == 2 ) {
        mode_m = MODE::SEO;
    } else if ( initialTotalNum_m > 2 ) {
        mode_m = MODE::BUNCH;
    } else
        mode_m = MODE::UNDEFINED;
192

193 194 195 196 197 198 199 200
    if ( timeIntegrator == 0 ) {
        stepper_m = stepper::INTEGRATOR::RK4;
    } else if ( timeIntegrator == 1) {
        stepper_m = stepper::INTEGRATOR::LF2;
    } else if ( timeIntegrator == 2) {
        stepper_m = stepper::INTEGRATOR::MTS;
    } else
        stepper_m = stepper::INTEGRATOR::UNDEFINED;
gsell's avatar
gsell committed
201 202 203 204 205 206 207
}

/**
 * Destructor ParallelCyclotronTracker
 *
 */
ParallelCyclotronTracker::~ParallelCyclotronTracker() {
208 209
    for(Component* component : myElements) {
        delete(component);
gsell's avatar
gsell committed
210
    }
211 212
    for(auto fd : FieldDimensions) {
        delete(fd);
gsell's avatar
gsell committed
213 214
    }
    delete itsBeamline;
215
    // delete opalRing_m;
gsell's avatar
gsell committed
216 217
}

218 219 220 221 222 223
/**
 * AAA
 *
 * @param none
 */
void ParallelCyclotronTracker::initializeBoundaryGeometry() {
224 225
    for(Component * component : myElements) {
        bgf_m = dynamic_cast<ElementBase *>(component)->getBoundaryGeometry();
226 227 228 229 230 231 232 233 234 235
        if(!bgf_m)
            continue;
        else
            break;
    }
    if (bgf_m) {
        itsDataSink->writeGeomToVtk(*bgf_m, std::string("data/testGeometry-00000.vtk"));
        OpalData::getInstance()->setGlobalGeometry(bgf_m);
        *gmsg << "* Boundary geometry initialized " << endl;
    }
236
}
237

238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259

/// set the working sub-mode for multi-bunch mode: "FORCE" or "AUTO"
inline
void ParallelCyclotronTracker::setMultiBunchMode(const std::string& mbmode)
{
    if ( mbmode.compare("FORCE") == 0 ) {
        *gmsg << "FORCE mode: The multi bunches will be injected consecutively "
              << "after each revolution, until get \"TURNS\" bunches." << endl;
        multiBunchMode_m = MB_MODE::FORCE;
    } else if ( mbmode.compare("AUTO") == 0 ) {
        *gmsg << "AUTO mode: The multi bunches will be injected only when "
                      << "the distance between two neighboring bunches " << endl
                      << "is below the limitation. The control parameter is set to "
                      << CoeffDBunches_m << endl;
        multiBunchMode_m = MB_MODE::AUTO;
    } else if ( mbmode.compare("NONE") == 0 )
        multiBunchMode_m = MB_MODE::NONE;
    else
        throw OpalException("ParallelCyclotronTracker::setMultiBunchMode()",
                            "MBMODE name \"" + mbmode + "\" unknown.");
}

260 261 262 263 264 265
/**
 *
 *
 * @param fn Base file name
 */
void ParallelCyclotronTracker::bgf_main_collision_test() {
266
    if(!bgf_m) return;
267

268
    Inform msg("bgf_main_collision_test ");
269

270 271 272
    /**
     *Here we check if a particles is outside the domain, flag it for deletion
     */
273

274
    Vector_t intecoords = 0.0;
275

276
    // This has to match the dT in the rk4 pusher
frey_m's avatar
frey_m committed
277
    double dtime = itsBunch_m->getdT() * getHarmonicNumber();
278

279
    int triId = 0;
frey_m's avatar
frey_m committed
280 281 282
    for(size_t i = 0; i < itsBunch_m->getLocalNum(); i++) {
        int res = bgf_m->partInside(itsBunch_m->R[i], itsBunch_m->P[i], dtime, itsBunch_m->PType[i], itsBunch_m->Q[i], intecoords, triId);
        //int res = bgf_m->partInside(itsBunch_m->R[i]*1.0e-3, itsBunch_m->P[i], dtime, itsBunch_m->PType[i], itsBunch_m->Q[i], intecoords, triId);
283
        if(res >= 0) {
frey_m's avatar
frey_m committed
284
            itsBunch_m->Bin[i] = -1;
285
        }
286
    }
287 288 289
}


gsell's avatar
gsell committed
290 291 292 293 294
/**
 *
 *
 * @param fn Base file name
 */
295
void ParallelCyclotronTracker::openFiles(std::string SfileName) {
gsell's avatar
gsell committed
296

297
    std::string  SfileName2 = SfileName + std::string("-Angle0.dat");
gsell's avatar
gsell committed
298
    outfTheta0_m.precision(8);
299
    outfTheta0_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
300
    outfTheta0_m.open(SfileName2.c_str());
301 302 303
    outfTheta0_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma" << std::endl;
gsell's avatar
gsell committed
304

305
    SfileName2 = SfileName + std::string("-Angle1.dat");
gsell's avatar
gsell committed
306
    outfTheta1_m.precision(8);
307
    outfTheta1_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
308
    outfTheta1_m.open(SfileName2.c_str());
309 310 311
    outfTheta1_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
312

313
    SfileName2 = SfileName + std::string("-Angle2.dat");
gsell's avatar
gsell committed
314
    outfTheta2_m.precision(8);
315
    outfTheta2_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
316
    outfTheta2_m.open(SfileName2.c_str());
317 318 319
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
320 321 322

    // for single Particle Mode, output after each turn, to define matched initial phase ellipse.

323
    SfileName2 = SfileName + std::string("-afterEachTurn.dat");
gsell's avatar
gsell committed
324 325

    outfThetaEachTurn_m.precision(8);
326
    outfThetaEachTurn_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
327 328

    outfThetaEachTurn_m.open(SfileName2.c_str());
329 330 331
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
332 333 334 335 336 337 338 339 340 341 342 343 344 345 346
}

/**
 * Close all files related to
 * special output in the Cyclotron
 * mode.
 */
void ParallelCyclotronTracker::closeFiles() {

    outfTheta0_m.close();
    outfTheta1_m.close();
    outfTheta2_m.close();
    outfThetaEachTurn_m.close();
}

347
/**
348 349 350
 *
 * @param ring
 */
kraus's avatar
kraus committed
351
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
352

353
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
354

snuverink_j's avatar
snuverink_j committed
355
    delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
356

kraus's avatar
kraus committed
357
    opalRing_m = dynamic_cast<Ring*>(ring.clone());
Daniel Winklehner's avatar
Daniel Winklehner committed
358

359
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
360

frey_m's avatar
frey_m committed
361
    opalRing_m->initialise(itsBunch_m);
362 363 364 365

    referenceR = opalRing_m->getBeamRInit();
    referencePr = opalRing_m->getBeamPRInit();
    referenceTheta = opalRing_m->getBeamPhiInit();
Daniel Winklehner's avatar
Daniel Winklehner committed
366

367
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
368
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
369 370
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
371 372

    referenceZ = 0.0;
373
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
374 375 376 377

    referencePtot = itsReference.getGamma() * itsReference.getBeta();
    referencePt = sqrt(referencePtot * referencePtot - referencePr * referencePr);

378 379
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
380

381 382
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
383

384
    double BcParameter[8];
Daniel Winklehner's avatar
Daniel Winklehner committed
385

386
    for(int i = 0; i < 8; i++)
Daniel Winklehner's avatar
Daniel Winklehner committed
387 388
        BcParameter[i] = 0.0;

389
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405

    // Finally print some diagnostic
    *gmsg << "* Initial beam radius = " << referenceR << " [mm] " << endl;
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
    *gmsg << "* Total reference momentum   = " << referencePtot * 1000.0
          << " [MCU]" << endl;
    *gmsg << "* Reference azimuthal momentum  = " << referencePt * 1000.0
          << " [MCU]" << endl;
    *gmsg << "* Reference radial momentum     = " << referencePr * 1000.0
          << " [MCU]" << endl;
    *gmsg << "* " << opalRing_m->getSymmetry() << " fold field symmetry "
          << endl;
    *gmsg << "* Harmonic number h= " << opalRing_m->getHarmonicNumber() << " "
          << endl;
}
gsell's avatar
gsell committed
406 407 408 409 410 411 412 413

/**
 *
 *
 * @param cycl
 */
void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {

414
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
415

416 417
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
418

419
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
420 421 422 423 424 425
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
426
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
427 428
        *gmsg         << "*     (Use BANDRF type cyclotron and use RFMAPFN to load both magnetic" << endl;
        *gmsg         << "*     and electric fields, setting SUPERPOSE to an array of TRUE values.)" << endl;
snuverink_j's avatar
snuverink_j committed
429
        *gmsg         << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
430 431
        *gmsg         << "*     FFT does not give the correct results (boundaty conditions are missing)." << endl;
        *gmsg         << "* 3.) The whole geometry will be meshed and used for the fieldsolve." << endl;
432 433
        *gmsg         << "*     There will be no transformations of the bunch into a local frame und consequently," << endl;
        *gmsg         << "*     the problem will be treated non-relativistically!" << endl;
434 435 436 437 438 439 440
        *gmsg         << "*     (This is not an issue for spiral inflectors as they are typically < 100 keV/amu.)" << endl;
        *gmsg << endl << "* Note: For now, multi-bunch mode (MBM) needs to be de-activated for spiral inflector" << endl;
        *gmsg         << "* and space charge needs to be solved every time-step. numBunch_m and scSolveFreq are reset." << endl;
        numBunch_m = 1;

    }

441
    // Fresh run (no restart):
442
    if(!OpalData::getInstance()->inRestartRun()) {
443

444
        // Get reference values from cyclotron element
445
        // For now, these are still stored in mm. should be the only ones. -DW
446 447 448 449
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
450 451
        referencePz    = elptr->getPZinit();

452
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
453 454
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
455 456 457 458 459
        }

        referencePtot =  itsReference.getGamma() * itsReference.getBeta();

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
460
        float insqrt = referencePtot * referencePtot - \
461
            referencePr * referencePr - referencePz * referencePz;
462 463 464 465 466 467 468 469 470

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

	        referencePt = 0.0;

            } else {

471
	        throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
472
                                    "Pt imaginary!");
473 474 475 476 477 478 479
            }

        } else {

            referencePt = sqrt(insqrt);
        }

480
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
481
            referencePt *= -1.0;
482 483
        // End calculate referencePt

484
        // Restart a run:
485 486
    } else {

487
        // If the user wants to save the restarted run in local frame,
488
        // make sure the previous h5 file was local too
489
      if (Options::psDumpFrame != Options::GLOBAL) {
winklehner_d's avatar
-DW  
winklehner_d committed
490
        if (!previousH5Local) {
491 492
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a local restart from a global h5 file!");
winklehner_d's avatar
-DW  
winklehner_d committed
493
        }
494 495
            // Else, if the user wants to save the restarted run in global frame,
            // make sure the previous h5 file was global too
winklehner_d's avatar
-DW  
winklehner_d committed
496 497
      } else {
        if (previousH5Local) {
498 499
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
500
            }
501
        }
502

503
        // Adjust some of the reference variables from the h5 file
504 505
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
506
        referencePtot = bega;
507
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
508
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
509 510
                                "PHIINIT is out of [-180, 180)!");
        }
511 512
    }

513 514
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
515

516
    *gmsg << endl;
adelmann's avatar
adelmann committed
517
    *gmsg << "* Bunch global starting position:" << endl;
518 519
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
520
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
521
    *gmsg << endl;
adelmann's avatar
adelmann committed
522
    *gmsg << "* Bunch global starting momenta:" << endl;
gsell's avatar
gsell committed
523 524
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
525
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
526 527 528
    *gmsg << "* Reference azimuthal momentum (Pt) = " << referencePt * 1000.0 << " [MCU]" << endl;
    *gmsg << "* Reference radial momentum (Pr) = " << referencePr * 1000.0 << " [MCU]" << endl;
    *gmsg << "* Reference axial momentum (Pz) = " << referencePz * 1000.0 << " [MCU]" << endl;
529
    *gmsg << endl;
adelmann's avatar
adelmann committed
530

gsell's avatar
gsell committed
531
    double sym = elptr->getSymmetry();
snuverink_j's avatar
snuverink_j committed
532
    *gmsg << "* " << sym << "-fold field symmetry " << endl;
gsell's avatar
gsell committed
533

534 535 536
    // ckr: this just returned the default value as defined in Component.h
    // double rff = elptr->getRfFrequ();
    // *gmsg << "* Rf frequency= " << rff << " [MHz]" << endl;
gsell's avatar
gsell committed
537

538
    std::string fmfn = elptr->getFieldMapFN();
539
    *gmsg << "* Field map file name = " << fmfn << " " << endl;
gsell's avatar
gsell committed
540

541
    std::string type = elptr->getCyclotronType();
542
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
543

544 545
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
adelmann's avatar
adelmann committed
546
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
547 548 549

    double zmin = elptr->getMinZ();
    double zmax = elptr->getMaxZ();
adelmann's avatar
adelmann committed
550
    *gmsg << "* Vertical aperture = " << zmin << " ... " << zmax<<" [m]"<< endl;
gsell's avatar
gsell committed
551 552

    double h = elptr->getCyclHarm();
553
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
554
    *gmsg << "* Harmonic number h = " << h << " " << endl;
555

gsell's avatar
gsell committed
556 557 558 559 560 561 562 563 564
    /**
     * To ease the initialise() function, set a integral parameter fieldflag internally.
     * Its value is  by the option "TYPE" of the element  "CYCLOTRON"
     * fieldflag = 1, readin PSI format measured field file (default)
     * fieldflag = 2, readin carbon cyclotron field file created by Jianjun Yang, TYPE=CARBONCYCL
     * fieldflag = 3, readin ANSYS format file for CYCIAE-100 created by Jianjun Yang, TYPE=CYCIAE
     * fieldflag = 4, readin AVFEQ format file for Riken cyclotrons
     * fieldflag = 5, readin FFAG format file for MSU/FNAL FFAG
     * fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron
565
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
566 567
     */
    int  fieldflag;
568
    if(type == std::string("CARBONCYCL")) {
gsell's avatar
gsell committed
569
        fieldflag = 2;
570
    } else if(type == std::string("CYCIAE")) {
gsell's avatar
gsell committed
571
        fieldflag = 3;
572
    } else if(type == std::string("AVFEQ")) {
gsell's avatar
gsell committed
573
        fieldflag = 4;
574
    } else if(type == std::string("FFAG")) {
gsell's avatar
gsell committed
575
        fieldflag = 5;
576
    } else if(type == std::string("BANDRF")) {
gsell's avatar
gsell committed
577
        fieldflag = 6;
578 579
    } else if(type == std::string("SYNCHROCYCLOTRON")) {
	fieldflag = 7;
580
    } else //(type == "RING")
gsell's avatar
gsell committed
581 582
        fieldflag = 1;

583
    // Read in cyclotron field maps (midplane + 3D fields if desired).
frey_m's avatar
frey_m committed
584
    elptr->initialise(itsBunch_m, fieldflag, elptr->getBScale());
gsell's avatar
gsell committed
585 586

    double BcParameter[8];
587

588
    for(int i = 0; i < 8; i++)
589
        BcParameter[i] = 0.0;
590

591 592
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
593

594
    // Store inner radius and outer radius of cyclotron field map in the list
595
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612
}

/**
 * Not implemented and most probable never used
 *
 */
void ParallelCyclotronTracker::visitBeamBeam(const BeamBeam &) {
    *gmsg << "In BeamBeam tracker is missing " << endl;
}

/**
 *
 *
 * @param coll
 */
void ParallelCyclotronTracker::visitCollimator(const Collimator &coll) {

613
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
614

615 616
    Collimator* elptr = dynamic_cast<Collimator *>(coll.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
617

618
    double xstart = elptr->getXStart();
adelmann's avatar
adelmann committed
619
    *gmsg << "* Xstart= " << xstart << " [mm]" << endl;
gsell's avatar
gsell committed
620

621
    double xend = elptr->getXEnd();
adelmann's avatar
adelmann committed
622
    *gmsg << "* Xend= " << xend << " [mm]" << endl;
gsell's avatar
gsell committed
623

624
    double ystart = elptr->getYStart();
adelmann's avatar
adelmann committed
625
    *gmsg << "* Ystart= " << ystart << " [mm]" << endl;
gsell's avatar
gsell committed
626

627
    double yend = elptr->getYEnd();
628
    *gmsg << "* Yend= " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
629

630
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
631
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
632 633

    double zend = elptr->getZEnd();
634
    *gmsg << "* Zend= " << zend << " [mm]" << endl;
635

636
    double width = elptr->getWidth();
adelmann's avatar
adelmann committed
637
    *gmsg << "* Width= " << width << " [mm]" << endl;
gsell's avatar
gsell committed
638

frey_m's avatar
frey_m committed
639
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
640 641 642 643

    double BcParameter[8];
    for(int i = 0; i < 8; i++)
        BcParameter[i] = 0.0;
644

645 646 647 648 649 650
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
    BcParameter[4] = 0.001 * width ;

651
    buildupFieldList(BcParameter, ElementBase::COLLIMATOR, elptr);
gsell's avatar
gsell committed
652 653 654 655 656 657 658 659 660 661 662 663
}

/**
 *
 *
 * @param corr
 */
void ParallelCyclotronTracker::visitCorrector(const Corrector &corr) {
    *gmsg << "In Corrector; L= " << corr.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Corrector *>(corr.clone()));
}

adelmann's avatar
adelmann committed
664 665 666 667 668 669 670 671 672 673 674 675
/**
 *
 *
 * @param degrader
 */
void ParallelCyclotronTracker::visitDegrader(const Degrader &deg) {
    *gmsg << "In Degrader; L= " << deg.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Degrader *>(deg.clone()));

}


gsell's avatar
gsell committed
676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705
/**
 *
 *
 * @param diag
 */
void ParallelCyclotronTracker::visitDiagnostic(const Diagnostic &diag) {
    *gmsg << "In Diagnostic; L= " << diag.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Diagnostic *>(diag.clone()));
}

/**
 *
 *
 * @param drift
 */
void ParallelCyclotronTracker::visitDrift(const Drift &drift) {
    *gmsg << "In drift L= " << drift.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Drift *>(drift.clone()));
}

/**
 *
 *
 * @param lamb
 */
void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) {
    *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Lambertson *>(lamb.clone()));
}

706 707 708
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
709 710
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
711 712
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
713
                                opalRing_m->getNextNormal());
714 715 716
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738
/**
 *
 *
 * @param marker
 */
void ParallelCyclotronTracker::visitMarker(const Marker &marker) {
    //   *gmsg << "In Marker; L= " << marker.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Marker *>(marker.clone()));
    // Do nothing.
}

/**
 *
 *
 * @param corr
 */
void ParallelCyclotronTracker::visitMonitor(const Monitor &corr) {
    //   *gmsg << "In Monitor; L= " << corr.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Monitor *>(corr.clone()));
    //   applyDrift(flip_s * corr.getElementLength());
}

739

gsell's avatar
gsell committed
740 741 742 743 744 745 746 747 748 749
/**
 *
 *
 * @param mult
 */
void ParallelCyclotronTracker::visitMultipole(const Multipole &mult) {
    *gmsg << "In Multipole; L= " << mult.getElementLength() << " however the element is missing " << endl;
    myElements.push_back(dynamic_cast<Multipole *>(mult.clone()));
}

ext-rogers_c's avatar
ext-rogers_c committed
750 751 752 753 754 755 756
/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
757
    if (opalRing_m != NULL) {
ext-rogers_c's avatar
ext-rogers_c committed
758
        opalRing_m->appendElement(multT);
759
    } else {
ext-rogers_c's avatar
ext-rogers_c committed
760 761
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
762
    }
ext-rogers_c's avatar
ext-rogers_c committed
763 764 765
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

gsell's avatar
gsell committed
766 767 768 769 770 771
/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
772
    *gmsg << "* -----------  Probe -------------------------------" << endl;
773 774
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
775

776
    double xstart = elptr->getXstart();
777
    *gmsg << "XStart= " << xstart << " [mm]" << endl;
gsell's avatar
gsell committed
778

779
    double xend = elptr->getXend();
780
    *gmsg << "XEnd= " << xend << " [mm]" << endl;
gsell's avatar
gsell committed
781

782
    double ystart = elptr->getYstart();
783
    *gmsg << "YStart= " << ystart << " [mm]" << endl;
gsell's avatar
gsell committed
784

785
    double yend = elptr->getYend();
786
    *gmsg << "YEnd= " << yend << " [mm]" << endl;
gsell's avatar
gsell committed
787

788
    double width = elptr->getWidth();
789
    *gmsg << "Width= " << width << " [mm]" << endl;
gsell's avatar
gsell committed
790 791 792


    // initialise, do nothing
frey_m's avatar
frey_m committed
793
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
794 795 796 797

    double BcParameter[8];
    for(int i = 0; i < 8; i++)
        BcParameter[i] = 0.0;
798

799 800 801 802 803
    BcParameter[0] = 0.001 * xstart ;
    BcParameter[1] = 0.001 * xend;
    BcParameter[2] = 0.001 * ystart ;
    BcParameter[3] = 0.001 * yend;
    BcParameter[4] = 0.001 * width ;
gsell's avatar
gsell committed
804 805

    // store probe parameters in the list
806
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
807 808 809 810 811 812 813 814 815 816 817 818
}

/**
 *
 *
 * @param bend
 */
void ParallelCyclotronTracker::visitRBend(const RBend &bend) {
    *gmsg << "In RBend; L= " << bend.getElementLength() << " however the element is missing " << endl;
    myElements.push_back(dynamic_cast<RBend *>(bend.clone()));
}

819 820 821 822 823 824
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
825
                            "Need to define a RINGDEFINITION to use SBend3D element");
826 827
}

828 829
void ParallelCyclotronTracker::visitScalingFFAGMagnet(const ScalingFFAGMagnet &bend) {
    *gmsg << "Adding ScalingFFAGMagnet" << endl;
ext-rogers_c's avatar
ext-rogers_c committed
830
    if (opalRing_m != NULL) {
831
        opalRing_m->appendElement(bend);
ext-rogers_c's avatar
ext-rogers_c committed
832
    } else {
833 834
        throw OpalException("ParallelCyclotronTracker::visitScalingFFAGMagnet",
                            "Need to define a RINGDEFINITION to use ScalingFFAGMagnet element");
ext-rogers_c's avatar
ext-rogers_c committed
835
    }
836 837
}

838 839 840 841 842 843
void ParallelCyclotronTracker::visitVariableRFCavity(const VariableRFCavity &cav) {
    *gmsg << "Adding Variable RF Cavity" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(cav);
    else
        throw OpalException("ParallelCyclotronTracker::visitVariableRFCavity",
844
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
845 846
}

gsell's avatar
gsell committed
847 848 849 850 851 852 853 854
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

    *gmsg << "* --------- RFCavity ------------------------------" << endl;
855

856 857
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
858 859 860 861 862 863 864

    if((elptr->getComponentType() != "SINGLEGAP") && (elptr->getComponentType() != "DOUBLEGAP")) {
        *gmsg << (elptr->getComponentType()) << endl;
        throw OpalException("ParallelCyclotronTracker::visitRFCavity",
                            "The ParallelCyclotronTracker can only play with cyclotron type RF system currently ...");
    }

865
    double rmin = elptr->getRmin();
gsell's avatar
gsell committed
866 867
    *gmsg << "* Minimal radius of cavity= " << rmin << " [mm]" << endl;

868
    double rmax = elptr->getRmax();
gsell's avatar
gsell committed
869 870
    *gmsg << "* Maximal radius of cavity= " << rmax << " [mm]" << endl;

871
    double rff = elptr->getCycFrequency();
gsell's avatar
gsell committed
872 873
    *gmsg << "* RF frequency (2*pi*f)= " << rff << " [rad/s]" << endl;

874
    std::string fmfn = elptr->getFieldMapFN();
gsell's avatar
gsell committed
875
    *gmsg << "* RF Field map file name= " << fmfn << endl;
876

877
    double angle = elptr->getAzimuth();
gsell's avatar
gsell committed
878 879
    *gmsg << "* Cavity azimuth position= " << angle << " [deg] " << endl;

880
    double gap = elptr->getGapWidth();
gsell's avatar
gsell committed
881 882
    *gmsg << "* Cavity gap width= " << gap << " [mm] " << endl;

883
    double pdis = elptr->getPerpenDistance();
gsell's avatar
gsell committed
884 885
    *gmsg << "* Cavity Shift distance= " << pdis << " [mm] " << endl;

886
    double phi0 = elptr->getPhi0();
gsell's avatar
gsell committed
887 888
    *gmsg << "* Initial RF phase (t=0)= " << phi0 << " [deg] " << endl;

889 890 891
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
892
    */
893

894 895 896 897
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

898
    dvector_t  unityVec;
899 900 901 902
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
903

904
    if (elptr->getFrequencyModelName() != "") {
905 906
        freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName());
        *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl;
907
    }
908 909
    else
        freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
910 911

    if (elptr->getAmplitudeModelName() != "") {
912 913
        ampl_atd = AbstractTimeDependence::getTimeDependence(