ParallelCyclotronTracker.cpp 124 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 42 43 44 45 46
#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/RFQuadrupole.h"
#include "AbsBeamline/SBend.h"
47
#include "AbsBeamline/SBend3D.h"
gsell's avatar
gsell committed
48 49 50 51 52
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/CyclotronValley.h"
#include "AbsBeamline/Stripper.h"
53
#include "AbsBeamline/VariableRFCavity.h"
54

55 56
#include "AbstractObjects/Element.h"

57
#include "Elements/OpalBeamline.h"
kraus's avatar
kraus committed
58
#include "AbsBeamline/Ring.h"
gsell's avatar
gsell committed
59 60 61

#include "BeamlineGeometry/Euclid3D.h"
#include "BeamlineGeometry/PlanarArcGeometry.h"
Jianjun Yang's avatar
Jianjun Yang committed
62
#include "BeamlineGeometry/RBendGeometry.h"
gsell's avatar
gsell committed
63 64 65 66 67 68 69 70 71 72 73 74
#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/NumToStr.h"
#include "Utilities/OpalException.h"

75
#include "BasicActions/DumpFields.h"
76
#include "BasicActions/DumpEMFields.h"
77

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

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

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

class Beamline;
class PartData;
90

gsell's avatar
gsell committed
91 92 93
using Physics::pi;
using Physics::q_e;

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

96 97 98
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
99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114

#define PSdim 6

extern Inform *gmsg;

// typedef FVector<double, PSdim> Vector;

/**
 * Constructor ParallelCyclotronTracker
 *
 * @param beamline
 * @param reference
 * @param revBeam
 * @param revTrack
 */
ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
115 116
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack):
gsell's avatar
gsell committed
117
    Tracker(beamline, reference, revBeam, revTrack),
118
    lastDumpedStep_m(0),
119
    eta_m(0.01),
gsell's avatar
gsell committed
120 121
    myNode_m(Ippl::myNode()),
    initialLocalNum_m(0),
122
    initialTotalNum_m(0),
123 124 125 126
    opalRing_m(NULL),
    itsStepper_mp(nullptr),
    mode_m(MODE::UNDEFINED),
    stepper_m(stepper::INTEGRATOR::UNDEFINED) {
gsell's avatar
gsell committed
127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
    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,
                                                   PartBunch &bunch,
                                                   DataSink &ds,
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack,
147
                                                   int maxSTEPS, int timeIntegrator):
gsell's avatar
gsell committed
148 149
    Tracker(beamline, reference, revBeam, revTrack),
    maxSteps_m(maxSTEPS),
150
    lastDumpedStep_m(0),
151
    eta_m(0.01),
gsell's avatar
gsell committed
152 153
    myNode_m(Ippl::myNode()),
    initialLocalNum_m(bunch.getLocalNum()),
154
    initialTotalNum_m(bunch.getTotalNum()),
155 156
    opalRing_m(NULL),
    itsStepper_mp(nullptr) {
gsell's avatar
gsell committed
157 158 159 160 161 162 163 164 165 166 167
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsBunch = &bunch;
    itsDataSink = &ds;
    //  scaleFactor_m = itsBunch->getdT() * c;
    scaleFactor_m = 1;
    multiBunchMode_m = 0;

    IntegrationTimer_m = IpplTimings::getTimer("Integration");
    TransformTimer_m   = IpplTimings::getTimer("Frametransform");
    DumpTimer_m        = IpplTimings::getTimer("Dump");
    BinRepartTimer_m   = IpplTimings::getTimer("Binaryrepart");
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186
    
    // 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;
    
    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
187 188 189 190 191 192 193
}

/**
 * Destructor ParallelCyclotronTracker
 *
 */
ParallelCyclotronTracker::~ParallelCyclotronTracker() {
194 195
    for(Component* component : myElements) {
        delete(component);
gsell's avatar
gsell committed
196
    }
197 198
    for(auto fd : FieldDimensions) {
        delete(fd);
gsell's avatar
gsell committed
199 200
    }
    delete itsBeamline;
201
    // delete opalRing_m;
gsell's avatar
gsell committed
202 203
}

204 205 206 207 208 209
/**
 * AAA
 *
 * @param none
 */
void ParallelCyclotronTracker::initializeBoundaryGeometry() {
210 211
    for(Component * component : myElements) {
        bgf_m = dynamic_cast<ElementBase *>(component)->getBoundaryGeometry();
212 213 214 215 216 217 218 219 220 221
        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;
    }
222
}
223

224 225 226 227 228 229
/**
 *
 *
 * @param fn Base file name
 */
void ParallelCyclotronTracker::bgf_main_collision_test() {
230
    if(!bgf_m) return;
231

232
    Inform msg("bgf_main_collision_test ");
233

234 235 236
    /**
     *Here we check if a particles is outside the domain, flag it for deletion
     */
237

238
    Vector_t intecoords = 0.0;
239

240 241
    // This has to match the dT in the rk4 pusher
    double dtime = itsBunch->getdT() * getHarmonicNumber();
242

243 244
    int triId = 0;
    for(size_t i = 0; i < itsBunch->getLocalNum(); i++) {
245 246
        int res = bgf_m->partInside(itsBunch->R[i], itsBunch->P[i], dtime, itsBunch->PType[i], itsBunch->Q[i], intecoords, triId);
        //int res = bgf_m->partInside(itsBunch->R[i]*1.0e-3, itsBunch->P[i], dtime, itsBunch->PType[i], itsBunch->Q[i], intecoords, triId);
247 248 249
        if(res >= 0) {
            itsBunch->Bin[i] = -1;
        }
250
    }
251 252 253
}


gsell's avatar
gsell committed
254 255 256 257 258
/**
 *
 *
 * @param fn Base file name
 */
259
void ParallelCyclotronTracker::openFiles(std::string SfileName) {
gsell's avatar
gsell committed
260

261
    std::string  SfileName2 = SfileName + std::string("-Angle0.dat");
gsell's avatar
gsell committed
262
    outfTheta0_m.precision(8);
263
    outfTheta0_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
264
    outfTheta0_m.open(SfileName2.c_str());
265 266 267
    outfTheta0_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma" << std::endl;
gsell's avatar
gsell committed
268

269
    SfileName2 = SfileName + std::string("-Angle1.dat");
gsell's avatar
gsell committed
270
    outfTheta1_m.precision(8);
271
    outfTheta1_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
272
    outfTheta1_m.open(SfileName2.c_str());
273 274 275
    outfTheta1_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
276

277
    SfileName2 = SfileName + std::string("-Angle2.dat");
gsell's avatar
gsell committed
278
    outfTheta2_m.precision(8);
279
    outfTheta2_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
280
    outfTheta2_m.open(SfileName2.c_str());
281 282 283
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
284 285 286

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

287
    SfileName2 = SfileName + std::string("-afterEachTurn.dat");
gsell's avatar
gsell committed
288 289

    outfThetaEachTurn_m.precision(8);
290
    outfThetaEachTurn_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
291 292

    outfThetaEachTurn_m.open(SfileName2.c_str());
293 294 295
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
296 297 298 299 300 301 302 303 304 305 306 307 308 309 310
}

/**
 * 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();
}

311
/**
312 313 314
 *
 * @param ring
 */
kraus's avatar
kraus committed
315
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
316

317
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
318 319

    if (opalRing_m != NULL)
320
        delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
321

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

324
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
325

326 327 328 329 330
    opalRing_m->initialise(itsBunch);

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

332
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
333
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
334 335
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
336 337

    referenceZ = 0.0;
338
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
339 340 341 342

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

343 344
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
345

346 347
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
348

349
    double BcParameter[8];
Daniel Winklehner's avatar
Daniel Winklehner committed
350

351
    for(int i = 0; i < 8; i++)
Daniel Winklehner's avatar
Daniel Winklehner committed
352 353
        BcParameter[i] = 0.0;

354
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370

    // 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
371 372 373 374 375 376 377 378

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

379
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
380

381 382
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
383

384
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
385 386 387 388 389 390
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
391
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
392 393
        *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;
394
        *gmsg         << "* 2.) For high currentst is strongly recommended to use the SAAMG fieldsolver," << endl;
395 396
        *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;
397 398
        *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;
399 400 401 402 403 404 405
        *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;

    }

406
    // Fresh run (no restart):
407
    if(!OpalData::getInstance()->inRestartRun()) {
408

409 410
        // Get reference values from cyclotron element 
        // For now, these are still stored in mm. should be the only ones. -DW
411 412 413 414
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
415 416
        referencePz    = elptr->getPZinit();

417
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
418 419
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
420 421 422 423 424
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
425
        float insqrt = referencePtot * referencePtot - \
426
            referencePr * referencePr - referencePz * referencePz;
427 428 429 430 431 432 433 434 435

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

	        referencePt = 0.0;

            } else {

436
	        throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
437
                                    "Pt imaginary!");
438 439 440 441 442 443 444
            }

        } else {

            referencePt = sqrt(insqrt);
        }

445
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
446
            referencePt *= -1.0;
447 448
        // End calculate referencePt

449
        // Restart a run:
450 451
    } else {

452
        // If the user wants to save the restarted run in local frame,
453
        // make sure the previous h5 file was local too
winklehner_d's avatar
-DW  
winklehner_d committed
454 455
      if (Options::psDumpLocalFrame != Options::GLOBAL) {
        if (!previousH5Local) {
456 457
                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
458
        }
459 460
            // 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
461 462
      } else {
        if (previousH5Local) {
463 464
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
465
            }
466
        }
467

468
        // Adjust some of the reference variables from the h5 file
469 470
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
471
        referencePtot = bega;
472
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
473
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
474 475
                                "PHIINIT is out of [-180, 180)!");
        }
476 477
    }

478 479
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
480

481
    *gmsg << endl;
adelmann's avatar
adelmann committed
482
    *gmsg << "* Bunch global starting position:" << endl;
483 484
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
485
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
486
    *gmsg << endl;
adelmann's avatar
adelmann committed
487
    *gmsg << "* Bunch global starting momenta:" << endl;
gsell's avatar
gsell committed
488 489
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
490
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
491 492 493
    *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;
494
    *gmsg << endl;
adelmann's avatar
adelmann committed
495

gsell's avatar
gsell committed
496
    double sym = elptr->getSymmetry();
497
    *gmsg << "* " << sym << "-fold field symmerty " << endl;
gsell's avatar
gsell committed
498

499 500 501
    // 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
502

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

506
    std::string type = elptr->getCyclotronType();
507
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
508

509 510
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
511
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [mm] "<< endl;
512 513 514

    double zmin = elptr->getMinZ();
    double zmax = elptr->getMaxZ();
515
    *gmsg << "* Vertical aperture = " << zmin << " ... " << zmax<<" [mm]"<< endl;
gsell's avatar
gsell committed
516 517

    double h = elptr->getCyclHarm();
518
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
519
    *gmsg << "* Harmonic number h = " << h << " " << endl;
520

gsell's avatar
gsell committed
521 522 523 524 525 526 527 528 529
    /**
     * 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
530
     * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
gsell's avatar
gsell committed
531 532
     */
    int  fieldflag;
533
    if(type == std::string("CARBONCYCL")) {
gsell's avatar
gsell committed
534
        fieldflag = 2;
535
    } else if(type == std::string("CYCIAE")) {
gsell's avatar
gsell committed
536
        fieldflag = 3;
537
    } else if(type == std::string("AVFEQ")) {
gsell's avatar
gsell committed
538
        fieldflag = 4;
539
    } else if(type == std::string("FFAG")) {
gsell's avatar
gsell committed
540
        fieldflag = 5;
541
    } else if(type == std::string("BANDRF")) {
gsell's avatar
gsell committed
542
        fieldflag = 6;
543 544
    } else if(type == std::string("SYNCHROCYCLOTRON")) {
	fieldflag = 7;
545
    } else //(type == "RING")
gsell's avatar
gsell committed
546 547
        fieldflag = 1;

548
    // Read in cyclotron field maps (midplane + 3D fields if desired).
549
    elptr->initialise(itsBunch, fieldflag, elptr->getBScale());
gsell's avatar
gsell committed
550 551

    double BcParameter[8];
552

553
    for(int i = 0; i < 8; i++)
554
        BcParameter[i] = 0.0;
555

556 557
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
558

559
    // Store inner radius and outer radius of cyclotron field map in the list
560
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577
}

/**
 * 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) {

578
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
579

580 581
    Collimator* elptr = dynamic_cast<Collimator *>(coll.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
582

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

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

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

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

595
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
596
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
597 598

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

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

604
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
605 606 607 608

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

610 611 612 613 614 615
    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 ;

616
    buildupFieldList(BcParameter, ElementBase::COLLIMATOR, elptr);
gsell's avatar
gsell committed
617 618 619 620 621 622 623 624 625 626 627 628
}

/**
 *
 *
 * @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
629 630 631 632 633 634 635 636 637 638 639 640
/**
 *
 *
 * @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
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
/**
 *
 *
 * @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()));
}

671 672 673
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
674 675
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
676 677
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
678
                                opalRing_m->getNextNormal());
679 680 681
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703
/**
 *
 *
 * @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());
}

704

gsell's avatar
gsell committed
705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720
/**
 *
 *
 * @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()));
}

/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
721
    *gmsg << "* -----------  Probe -------------------------------" << endl;
722 723
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
724

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

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

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

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

737
    double width = elptr->getWidth();
738
    *gmsg << "Width= " << width << " [mm]" << endl;
gsell's avatar
gsell committed
739 740 741


    // initialise, do nothing
742
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
743 744 745 746

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

748 749 750 751 752
    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
753 754

    // store probe parameters in the list
755
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
756 757 758 759 760 761 762 763 764 765 766 767
}

/**
 *
 *
 * @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()));
}

768 769 770 771 772 773
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
774
                            "Need to define a RINGDEFINITION to use SBend3D element");
775 776
}

777 778 779 780 781 782
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",
783
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
784 785
}

gsell's avatar
gsell committed
786 787 788 789 790 791 792 793
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

795 796
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
797 798 799 800 801 802 803

    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 ...");
    }

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

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

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

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

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

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

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

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

828 829 830
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
831
    */
832

833 834 835 836
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

837
    dvector_t  unityVec;
838 839 840 841
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
842

843
    if (elptr->getFrequencyModelName() != "") {
844 845
        freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName());
        *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl;
846
    }
847 848
    else
        freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
849 850

    if (elptr->getAmplitudeModelName() != "") {
851 852
        ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName());
        *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl;
853
    }
854 855
    else
        ampl_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
856 857

    if (elptr->getPhaseModelName() != "") {
858 859
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
860
    }
861 862
    else
        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
863

gsell's avatar
gsell committed
864
    // read cavity voltage profile data from file.
865
    elptr->initialise(itsBunch, freq_atd, ampl_atd, phase_atd);
gsell's avatar
gsell committed
866 867 868 869

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

871 872
    BcParameter[0] = 0.001 * rmin;
    BcParameter[1] = 0.001 * rmax;
873
    BcParameter[2] = 0.001 * pdis;
gsell's avatar
gsell committed
874 875
    BcParameter[3] = angle;

876
    buildupFieldList(BcParameter, ElementBase::RFCAVITY, elptr);
gsell's avatar
gsell committed
877 878 879 880 881 882 883 884
}

/**
 *
 *
 * @param rfq
 */
void ParallelCyclotronTracker::visitRFQuadrupole(const RFQuadrupole &rfq) {
885
    *gmsg << "In RFQuadrupole; L = " << rfq.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
886 887 888 889 890 891 892 893 894
    myElements.push_back(dynamic_cast<RFQuadrupole *>(rfq.clone()));
}

/**
 *
 *
 * @param bend
 */
void ParallelCyclotronTracker::visitSBend(const SBend &bend) {
895
    *gmsg << "In SBend; L = " << bend.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914
    myElements.push_back(dynamic_cast<SBend *>(bend.clone()));
}

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

/**
 *
 *
 * @param sept
 */
void ParallelCyclotronTracker::visitSeptum(const Septum &sept) {
915

916
    *gmsg << endl << "* -----------------------------  Septum ------------------------------- *" << endl;
gsell's avatar
gsell committed
917

918 919
    Septum *elptr = dynamic_cast<Septum *>(sept.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
920

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

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

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

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

933
    double width = elptr->getWidth();
934
    *gmsg << "Width = " << width << " [mm]" << endl;
gsell's avatar
gsell committed
935 936 937


    // initialise, do nothing
938
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
939 940 941 942

    double BcParameter[8];
    for(