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 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

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

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

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

gsell's avatar
gsell committed
94 95
class Beamline;
class PartData;
96

gsell's avatar
gsell committed
97 98 99
using Physics::pi;
using Physics::q_e;

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

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

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

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

216 217 218 219 220 221
/**
 * AAA
 *
 * @param none
 */
void ParallelCyclotronTracker::initializeBoundaryGeometry() {
222 223
    for(Component * component : myElements) {
        bgf_m = dynamic_cast<ElementBase *>(component)->getBoundaryGeometry();
224 225 226 227 228 229 230 231 232 233
        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;
    }
234
}
235

236 237 238 239 240 241
/**
 *
 *
 * @param fn Base file name
 */
void ParallelCyclotronTracker::bgf_main_collision_test() {
242
    if(!bgf_m) return;
243

244
    Inform msg("bgf_main_collision_test ");
245

246 247 248
    /**
     *Here we check if a particles is outside the domain, flag it for deletion
     */
249

250
    Vector_t intecoords = 0.0;
251

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

255
    int triId = 0;
frey_m's avatar
frey_m committed
256 257 258
    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);
259
        if(res >= 0) {
frey_m's avatar
frey_m committed
260
            itsBunch_m->Bin[i] = -1;
261
        }
262
    }
263 264 265
}


gsell's avatar
gsell committed
266 267 268 269 270
/**
 *
 *
 * @param fn Base file name
 */
271
void ParallelCyclotronTracker::openFiles(std::string SfileName) {
gsell's avatar
gsell committed
272

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

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

289
    SfileName2 = SfileName + std::string("-Angle2.dat");
gsell's avatar
gsell committed
290
    outfTheta2_m.precision(8);
291
    outfTheta2_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
292
    outfTheta2_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

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

299
    SfileName2 = SfileName + std::string("-afterEachTurn.dat");
gsell's avatar
gsell committed
300 301

    outfThetaEachTurn_m.precision(8);
302
    outfThetaEachTurn_m.setf(std::ios::scientific, std::ios::floatfield);
gsell's avatar
gsell committed
303 304

    outfThetaEachTurn_m.open(SfileName2.c_str());
305 306 307
    outfTheta2_m << "#  r [mm]      beta_r*gamma       "
                 << "theta [mm]      beta_theta*gamma        "
                 << "z [mm]          beta_z*gamma"  << std::endl;
gsell's avatar
gsell committed
308 309 310 311 312 313 314 315 316 317 318 319 320 321 322
}

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

323
/**
324 325 326
 *
 * @param ring
 */
kraus's avatar
kraus committed
327
void ParallelCyclotronTracker::visitRing(const Ring &ring) {
Daniel Winklehner's avatar
Daniel Winklehner committed
328

329
    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;
Daniel Winklehner's avatar
Daniel Winklehner committed
330 331

    if (opalRing_m != NULL)
332
        delete opalRing_m;
Daniel Winklehner's avatar
Daniel Winklehner committed
333

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

336
    myElements.push_back(opalRing_m);
Daniel Winklehner's avatar
Daniel Winklehner committed
337

frey_m's avatar
frey_m committed
338
    opalRing_m->initialise(itsBunch_m);
339 340 341 342

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

344
    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
kraus's avatar
kraus committed
345
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
346 347
                            "PHIINIT is out of [-180, 180)!");
    }
Daniel Winklehner's avatar
Daniel Winklehner committed
348 349

    referenceZ = 0.0;
350
    referencePz = 0.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
351 352 353 354

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

355 356
    if(referencePtot < 0.0)
        referencePt *= -1.0;
Daniel Winklehner's avatar
Daniel Winklehner committed
357

358 359
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
gsell's avatar
gsell committed
360

361
    double BcParameter[8];
Daniel Winklehner's avatar
Daniel Winklehner committed
362

363
    for(int i = 0; i < 8; i++)
Daniel Winklehner's avatar
Daniel Winklehner committed
364 365
        BcParameter[i] = 0.0;

366
    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);
367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382

    // 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
383 384 385 386 387 388 389 390

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

391
    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
gsell's avatar
gsell committed
392

393 394
    Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(elptr);
395

396
    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
397 398 399 400 401 402
    // useful information
    spiral_flag = elptr->getSpiralFlag();

    if(spiral_flag) {

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

    }

418
    // Fresh run (no restart):
419
    if(!OpalData::getInstance()->inRestartRun()) {
420

421 422
        // Get reference values from cyclotron element 
        // For now, these are still stored in mm. should be the only ones. -DW
423 424 425 426
        referenceR     = elptr->getRinit();
        referenceTheta = elptr->getPHIinit();
        referenceZ     = elptr->getZinit();
        referencePr    = elptr->getPRinit();
427 428
        referencePz    = elptr->getPZinit();

429
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
430 431
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
432 433 434 435 436
        }

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

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
437
        float insqrt = referencePtot * referencePtot - \
438
            referencePr * referencePr - referencePz * referencePz;
439 440 441 442 443 444 445 446 447

        if(insqrt < 0) {

            if(insqrt > -1.0e-10) {

	        referencePt = 0.0;

            } else {

448
	        throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
449
                                    "Pt imaginary!");
450 451 452 453 454 455 456
            }

        } else {

            referencePt = sqrt(insqrt);
        }

457
        if(referencePtot < 0.0)
Daniel Winklehner's avatar
Daniel Winklehner committed
458
            referencePt *= -1.0;
459 460
        // End calculate referencePt

461
        // Restart a run:
462 463
    } else {

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

480
        // Adjust some of the reference variables from the h5 file
481 482
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
483
        referencePtot = bega;
484
        if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
485
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
486 487
                                "PHIINIT is out of [-180, 180)!");
        }
488 489
    }

490 491
    sinRefTheta_m = sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = cos(referenceTheta * Physics::deg2rad);
492

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

gsell's avatar
gsell committed
508
    double sym = elptr->getSymmetry();
509
    *gmsg << "* " << sym << "-fold field symmerty " << endl;
gsell's avatar
gsell committed
510

511 512 513
    // 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
514

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

518
    std::string type = elptr->getCyclotronType();
519
    *gmsg << "* Type of cyclotron = " << type << " " << endl;
520

521 522
    double rmin = elptr->getMinR();
    double rmax = elptr->getMaxR();
523
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [mm] "<< endl;
524 525 526

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

    double h = elptr->getCyclHarm();
530
    *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
531
    *gmsg << "* Harmonic number h = " << h << " " << endl;
532

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

560
    // Read in cyclotron field maps (midplane + 3D fields if desired).
frey_m's avatar
frey_m committed
561
    elptr->initialise(itsBunch_m, fieldflag, elptr->getBScale());
gsell's avatar
gsell committed
562 563

    double BcParameter[8];
564

565
    for(int i = 0; i < 8; i++)
566
        BcParameter[i] = 0.0;
567

568 569
    BcParameter[0] = 0.001 * elptr->getRmin();
    BcParameter[1] = 0.001 * elptr->getRmax();
gsell's avatar
gsell committed
570

571
    // Store inner radius and outer radius of cyclotron field map in the list
572
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
gsell's avatar
gsell committed
573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589
}

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

590
    *gmsg << "* --------- Collimator -----------------------------" << endl;
gsell's avatar
gsell committed
591

592 593
    Collimator* elptr = dynamic_cast<Collimator *>(coll.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
594

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

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

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

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

607
    double zstart = elptr->getZStart();
adelmann's avatar
adelmann committed
608
    *gmsg << "* Zstart= " << zstart << " [mm]" << endl;
609 610

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

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

frey_m's avatar
frey_m committed
616
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
617 618 619 620

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

622 623 624 625 626 627
    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 ;

628
    buildupFieldList(BcParameter, ElementBase::COLLIMATOR, elptr);
gsell's avatar
gsell committed
629 630 631 632 633 634 635 636 637 638 639 640
}

/**
 *
 *
 * @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
641 642 643 644 645 646 647 648 649 650 651 652
/**
 *
 *
 * @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
653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682
/**
 *
 *
 * @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()));
}

683 684 685
void ParallelCyclotronTracker::visitOffset(const Offset & off) {
    if (opalRing_m == NULL)
        throw OpalException(
686 687
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
688 689
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
690
                                opalRing_m->getNextNormal());
691 692 693
    opalRing_m->appendElement(off);
}

gsell's avatar
gsell committed
694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715
/**
 *
 *
 * @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());
}

716

gsell's avatar
gsell committed
717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732
/**
 *
 *
 * @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) {
733
    *gmsg << "* -----------  Probe -------------------------------" << endl;
734 735
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
736

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

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

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

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

749
    double width = elptr->getWidth();
750
    *gmsg << "Width= " << width << " [mm]" << endl;
gsell's avatar
gsell committed
751 752 753


    // initialise, do nothing
frey_m's avatar
frey_m committed
754
    elptr->initialise(itsBunch_m);
gsell's avatar
gsell committed
755 756 757 758

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

760 761 762 763 764
    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
765 766

    // store probe parameters in the list
767
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
gsell's avatar
gsell committed
768 769 770 771 772 773 774 775 776 777 778 779
}

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

780 781 782 783 784 785
void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
    *gmsg << "Adding SBend3D" << endl;
    if (opalRing_m != NULL)
        opalRing_m->appendElement(bend);
    else
        throw OpalException("ParallelCyclotronTracker::visitSBend3D",
786
                            "Need to define a RINGDEFINITION to use SBend3D element");
787 788
}

789 790 791 792 793 794
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",
795
                            "Need to define a RINGDEFINITION to use VariableRFCavity element");
796 797
}

gsell's avatar
gsell committed
798 799 800 801 802 803 804 805
/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

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

807 808
    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
809 810 811 812 813 814 815

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

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

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

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

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

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

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

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

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

840 841 842
    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
843
    */
844

845 846 847 848
    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

849
    dvector_t  unityVec;
850 851 852 853
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
854

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

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

    if (elptr->getPhaseModelName() != "") {
870 871
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
872
    }
873 874
    else
        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
875

gsell's avatar
gsell committed
876
    // read cavity voltage profile data from file.
frey_m's avatar
frey_m committed
877
    elptr->initialise(itsBunch_m, freq_atd, ampl_atd, phase_atd);
gsell's avatar
gsell committed
878 879 880 881

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

883 884
    BcParameter[0] = 0.001 * rmin;
    BcParameter[1] = 0.001 * rmax;
885
    BcParameter[2] = 0.001 * pdis;
gsell's avatar
gsell committed
886 887
    BcParameter[3] = angle;

888
    buildupFieldList(BcParameter, ElementBase::RFCAVITY, elptr);
gsell's avatar
gsell committed
889 890 891 892 893 894 895 896
}

/**
 *
 *
 * @param rfq
 */
void ParallelCyclotronTracker::visitRFQuadrupole(const RFQuadrupole &rfq) {
897
    *gmsg << "In RFQuadrupole; L = " << rfq.getElementLength() << " however the element is missing " << endl;
gsell's avatar
gsell committed
898 899 900 901 902 903 904 905 906
    myElements.push_back(dynamic_cast<RFQuadrupole *>(rfq.clone()));
}

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

928
    *gmsg << endl << "* -----------------------------  Septum ------------------------------- *" << endl;
gsell's avatar
gsell committed
929

930 931
    Septum *elptr = dynamic_cast<Septum *>(sept.clone());
    myElements.push_back(elptr);
gsell's avatar
gsell committed
932

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

936
    dou