ParallelCyclotronTracker.cpp 142 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 94 95 96
using Physics::m_p; // GeV
using Physics::PMASS;
using Physics::PCHARGE;
using Physics::pi;
using Physics::q_e;

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

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

#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,
118 119
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack):
gsell's avatar
gsell committed
120
    Tracker(beamline, reference, revBeam, revTrack),
121
    lastDumpedStep_m(0),
122
    eta_m(0.01),
gsell's avatar
gsell committed
123 124
    myNode_m(Ippl::myNode()),
    initialLocalNum_m(0),
125
    initialTotalNum_m(0),
126 127 128 129
    opalRing_m(NULL),
    itsStepper_mp(nullptr),
    mode_m(MODE::UNDEFINED),
    stepper_m(stepper::INTEGRATOR::UNDEFINED) {
gsell's avatar
gsell committed
130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
    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,
150
                                                   int maxSTEPS, int timeIntegrator):
gsell's avatar
gsell committed
151 152
    Tracker(beamline, reference, revBeam, revTrack),
    maxSteps_m(maxSTEPS),
153
    lastDumpedStep_m(0),
154
    eta_m(0.01),
gsell's avatar
gsell committed
155 156
    myNode_m(Ippl::myNode()),
    initialLocalNum_m(bunch.getLocalNum()),
157
    initialTotalNum_m(bunch.getTotalNum()),
158 159
    opalRing_m(NULL),
    itsStepper_mp(nullptr) {
gsell's avatar
gsell committed
160 161 162 163 164 165 166 167 168 169 170
    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");
171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189
    
    // 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
190 191 192 193 194 195 196
}

/**
 * Destructor ParallelCyclotronTracker
 *
 */
ParallelCyclotronTracker::~ParallelCyclotronTracker() {
197
    for(std::list<Component *>::iterator compindex = myElements.begin(); compindex != myElements.end(); compindex++) {
gsell's avatar
gsell committed
198 199 200 201 202 203
        delete(*compindex);
    }
    for(beamline_list::iterator fdindex = FieldDimensions.begin(); fdindex != FieldDimensions.end(); fdindex++) {
        delete(*fdindex);
    }
    delete itsBeamline;
204 205 206
    if (opalRing_m != NULL) {
        // delete opalRing_m;
    }
gsell's avatar
gsell committed
207 208
}

209 210 211 212 213 214
/**
 * AAA
 *
 * @param none
 */
void ParallelCyclotronTracker::initializeBoundaryGeometry() {
215 216 217 218 219 220 221 222 223 224 225 226
    for(std::list<Component *>::iterator compindex = myElements.begin(); compindex != myElements.end(); compindex++) {
        bgf_m = dynamic_cast<ElementBase *>(*compindex)->getBoundaryGeometry();
        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;
    }
227
}
228

229 230 231 232 233 234
/**
 *
 *
 * @param fn Base file name
 */
void ParallelCyclotronTracker::bgf_main_collision_test() {
235
    if(!bgf_m) return;
236

237
    Inform msg("bgf_main_collision_test ");
238

239 240 241
    /**
     *Here we check if a particles is outside the domain, flag it for deletion
     */
242

243
    Vector_t intecoords = 0.0;
244

245 246
    // This has to match the dT in the rk4 pusher
    double dtime = itsBunch->getdT() * getHarmonicNumber();
247

248 249
    int triId = 0;
    for(size_t i = 0; i < itsBunch->getLocalNum(); i++) {
250 251
        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);
252 253 254
        if(res >= 0) {
            itsBunch->Bin[i] = -1;
        }
255
    }
256 257 258
}


gsell's avatar
gsell committed
259 260 261 262 263
/**
 *
 *
 * @param fn Base file name
 */
264
void ParallelCyclotronTracker::openFiles(std::string SfileName) {
gsell's avatar
gsell committed
265

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

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

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

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

292
    SfileName2 = SfileName + std::string("-afterEachTurn.dat");
gsell's avatar
gsell committed
293 294

    outfThetaEachTurn_m.precision(8);
295
    outfThetaEachTurn_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
296 297

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

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

316
/**
317 318 319
 *
 * @param ring
 */
kraus's avatar
kraus committed
320
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
321

322
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
323 324

    if (opalRing_m != NULL)
325
        delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
326

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

329
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
330

331 332 333 334 335
    opalRing_m->initialise(itsBunch);

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

337
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
338
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
339 340
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
341 342

    referenceZ = 0.0;
343
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
344 345 346 347

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

348 349
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
350

351 352
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
353

354
    double BcParameter[8];
Daniel Winklehner's avatar
Daniel Winklehner committed
355

356
    for(int i = 0; i < 8; i++)
Daniel Winklehner's avatar
Daniel Winklehner committed
357 358
        BcParameter[i] = 0.0;

359
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375

    // 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
376 377 378 379 380 381 382 383

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

384
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
385

386 387
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
388

389
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
390 391 392 393 394 395
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

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

    }

411
    // Fresh run (no restart):
412
    if(!OpalData::getInstance()->inRestartRun()) {
413

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

422
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
423 424
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
425 426 427 428 429
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
430
        float insqrt = referencePtot * referencePtot - \
431
            referencePr * referencePr - referencePz * referencePz;
432 433 434 435 436 437 438 439 440

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

	        referencePt = 0.0;

            } else {

441
	        throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
442
                                    "Pt imaginary!");
443 444 445 446 447 448 449
            }

        } else {

            referencePt = sqrt(insqrt);
        }

450
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
451
            referencePt *= -1.0;
452 453
        // End calculate referencePt

454
        // Restart a run:
455 456
    } else {

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

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

483 484
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
485

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

gsell's avatar
gsell committed
501
    double sym = elptr->getSymmetry();
502
    *gmsg << "* " << sym << "-fold field symmerty " << endl;
gsell's avatar
gsell committed
503

504 505 506
    // 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
507

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

511
    std::string type = elptr->getCyclotronType();
512
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
513

514 515
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
516
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [mm] "<< endl;
517 518 519

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

    double h = elptr->getCyclHarm();
523
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
524
    *gmsg << "* Harmonic number h = " << h << " " << endl;
525

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

553
    // Read in cyclotron field maps (midplane + 3D fields if desired).
554
    elptr->initialise(itsBunch, fieldflag, elptr->getBScale());
gsell's avatar
gsell committed
555 556

    double BcParameter[8];
557

558
    for(int i = 0; i < 8; i++)
559
        BcParameter[i] = 0.0;
560

561 562
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
563

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

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

583
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
584

585 586
    Collimator* elptr = dynamic_cast<Collimator *>(coll.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
587

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

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

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

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

600
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
601
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
602 603

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

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

609
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
610 611 612 613

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

615 616 617 618 619 620
    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 ;

621
    buildupFieldList(BcParameter, ElementBase::COLLIMATOR, elptr);
gsell's avatar
gsell committed
622 623 624 625 626 627 628 629 630 631 632 633
}

/**
 *
 *
 * @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
634 635 636 637 638 639 640 641 642 643 644 645
/**
 *
 *
 * @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
646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675
/**
 *
 *
 * @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()));
}

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

gsell's avatar
gsell committed
687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708
/**
 *
 *
 * @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());
}

709

gsell's avatar
gsell committed
710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725
/**
 *
 *
 * @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) {
726
    *gmsg << "* -----------  Probe -------------------------------" << endl;
727 728
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
729

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

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

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

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

742
    double width = elptr->getWidth();
743
    *gmsg << "Width= " << width << " [mm]" << endl;
gsell's avatar
gsell committed
744 745 746


    // initialise, do nothing
747
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
748 749 750 751

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

753 754 755 756 757
    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
758 759

    // store probe parameters in the list
760
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
761 762 763 764 765 766 767 768 769 770 771 772
}

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

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

782 783 784 785 786 787
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",
788
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
789 790
}

gsell's avatar
gsell committed
791 792 793 794 795 796 797 798
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

800 801
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
802 803 804 805 806 807 808

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

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

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

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

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

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

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

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

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

833 834 835
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
836
    */
837

838 839 840 841
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

842 843 844 845 846
    std::vector<double>  unityVec;
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
847

848
    if (elptr->getFrequencyModelName() != "") {
849 850
        freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName());
        *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl;
851
    }
852 853
    else
        freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
854 855

    if (elptr->getAmplitudeModelName() != "") {
856 857
        ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName());
        *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl;
858
    }
859 860
    else
        ampl_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
861 862

    if (elptr->getPhaseModelName() != "") {
863 864
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
865
    }
866 867
    else
        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
868

gsell's avatar
gsell committed
869
    // read cavity voltage profile data from file.
870
    elptr->initialise(itsBunch, freq_atd, ampl_atd, phase_atd);
gsell's avatar
gsell committed
871 872 873 874

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

876 877
    BcParameter[0] = 0.001 * rmin;
    BcParameter[1] = 0.001 * rmax;
878
    BcParameter[2] = 0.001 * pdis;
gsell's avatar
gsell committed
879 880
    BcParameter[3] = angle;

881
    buildupFieldList(BcParameter, ElementBase::RFCAVITY, elptr);
gsell's avatar
gsell committed
882 883 884 885 886 887 888 889
}

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

/**
 *
 *
 * @param bend
 */
void ParallelCyclotronTracker::visitSBend(const SBend &bend) {
900
    *gmsg << "In SBend; L = " << bend.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919
    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) {
920

921
    *gmsg << endl << "* -----------------------------  Septum ------------------------------- *" << endl;
gsell's avatar
gsell committed
922

923 924
    Septum *elptr = dynamic_cast<Septum *>(sept.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
925

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

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

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

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

938
    double width = elptr->getWidth();
939
    *gmsg << "Width = " << width << " [mm]" << endl;
gsell's avatar
gsell committed
940 941 942


    // initialise, do nothing
943
    elptr->initialise(itsBunch);
gsell's avatar
gsell committed
944 945 946 947

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