RFCavity.cpp 23 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
// ------------------------------------------------------------------------
// $RCSfile: RFCavity.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: RFCavity
//   Defines the abstract interface for an accelerating structure.
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:31 $
// $Author: andreas adelmann $
//
// ------------------------------------------------------------------------

#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/BeamlineVisitor.h"
23
#include "Algorithms/PartBunch.h"
24
#include "Algorithms/PartPusher.h"
25
#include "Fields/Fieldmap.h"
26
#include "Utilities/GeneralClassicException.h"
27
#include "Utilities/Util.h"
gsell's avatar
gsell committed
28 29

#include "gsl/gsl_interp.h"
30
#include "gsl/gsl_spline.h"
gsell's avatar
gsell committed
31 32
#include <iostream>
#include <fstream>
33 34 35 36
#ifdef OPAL_NOCPLUSPLUS11_NULLPTR
#define nullptr NULL
#endif

gsell's avatar
gsell committed
37 38 39 40 41 42 43 44 45
extern Inform *gmsg;

using namespace std;

// Class RFCavity
// ------------------------------------------------------------------------

RFCavity::RFCavity():
    Component(),
46 47 48
    phase_td_m(nullptr),
    amplitude_td_m(nullptr),
    frequency_td_m(nullptr),
gsell's avatar
gsell committed
49 50
    filename_m(""),
    scale_m(1.0),
51
    scaleError_m(0.0),
gsell's avatar
gsell committed
52
    phase_m(0.0),
53
    phaseError_m(0.0),
gsell's avatar
gsell committed
54
    frequency_m(0.0),
55 56 57
    fast_m(true),
    autophaseVeto_m(false),
    designEnergy_m(-1.0),
gsell's avatar
gsell committed
58 59
    startField_m(0.0),
    endField_m(0.0),
60
    length_m(0.0),
gsell's avatar
gsell committed
61 62 63 64 65 66 67 68 69
    type_m(SW),
    rmin_m(0.0),
    rmax_m(0.0),
    angle_m(0.0),
    sinAngle_m(0.0),
    cosAngle_m(0.0),
    pdis_m(0.0),
    gapwidth_m(0.0),
    phi0_m(0.0),
70 71 72
    RNormal_m(nullptr),
    VrNormal_m(nullptr),
    DvDr_m(nullptr),
73
    num_points_m(0)
74
{
gsell's avatar
gsell committed
75 76 77 78 79 80
    setElType(isRF);
}


RFCavity::RFCavity(const RFCavity &right):
    Component(right),
81 82 83 84 85 86
    phase_td_m(right.phase_td_m),
    phase_name_m(right.phase_name_m),
    amplitude_td_m(right.amplitude_td_m),
    amplitude_name_m(right.amplitude_name_m),
    frequency_td_m(right.frequency_td_m),
    frequency_name_m(right.frequency_name_m),
gsell's avatar
gsell committed
87 88
    filename_m(right.filename_m),
    scale_m(right.scale_m),
89
    scaleError_m(right.scaleError_m),
gsell's avatar
gsell committed
90
    phase_m(right.phase_m),
91
    phaseError_m(right.phaseError_m),
gsell's avatar
gsell committed
92
    frequency_m(right.frequency_m),
93 94 95
    fast_m(right.fast_m),
    autophaseVeto_m(right.autophaseVeto_m),
    designEnergy_m(right.designEnergy_m),
gsell's avatar
gsell committed
96 97
    startField_m(right.startField_m),
    endField_m(right.endField_m),
98
    length_m(right.length_m),
gsell's avatar
gsell committed
99 100 101 102 103 104 105 106 107
    type_m(right.type_m),
    rmin_m(right.rmin_m),
    rmax_m(right.rmax_m),
    angle_m(right.angle_m),
    sinAngle_m(right.sinAngle_m),
    cosAngle_m(right.cosAngle_m),
    pdis_m(right.pdis_m),
    gapwidth_m(right.gapwidth_m),
    phi0_m(right.phi0_m),
108 109 110
    RNormal_m(nullptr),
    VrNormal_m(nullptr),
    DvDr_m(nullptr),
111
    num_points_m(right.num_points_m)
112
{
gsell's avatar
gsell committed
113 114 115 116
    setElType(isRF);
}


117
RFCavity::RFCavity(const std::string &name):
gsell's avatar
gsell committed
118
    Component(name),
119 120 121
    phase_td_m(nullptr),
    amplitude_td_m(nullptr),
    frequency_td_m(nullptr),
gsell's avatar
gsell committed
122 123
    filename_m(""),
    scale_m(1.0),
124
    scaleError_m(0.0),
gsell's avatar
gsell committed
125
    phase_m(0.0),
126
    phaseError_m(0.0),
gsell's avatar
gsell committed
127
    frequency_m(0.0),
128 129 130
    fast_m(true),
    autophaseVeto_m(false),
    designEnergy_m(-1.0),
gsell's avatar
gsell committed
131 132
    startField_m(0.0),
    endField_m(0.0),
133
    length_m(0.0),
gsell's avatar
gsell committed
134 135 136 137 138 139 140 141 142
    type_m(SW),
    rmin_m(0.0),
    rmax_m(0.0),
    angle_m(0.0),
    sinAngle_m(0.0),
    cosAngle_m(0.0),
    pdis_m(0.0),
    gapwidth_m(0.0),
    phi0_m(0.0),
143 144 145
    RNormal_m(nullptr),
    VrNormal_m(nullptr),
    DvDr_m(nullptr),
146 147 148
    //     RNormal_m(std::nullptr_t(NULL)),
    //     VrNormal_m(std::nullptr_t(NULL)),
    //     DvDr_m(std::nullptr_t(NULL)),
149
    num_points_m(0)
150
{
gsell's avatar
gsell committed
151 152 153 154 155 156 157
    setElType(isRF);
}


RFCavity::~RFCavity() {
    // FIXME: in deleteFielmak, a map find makes problems
    //       Fieldmap::deleteFieldmap(filename_m);
158
    //~ if(RNormal_m) {
159 160 161
    //~ delete[] RNormal_m;
    //~ delete[] VrNormal_m;
    //~ delete[] DvDr_m;
162
    //~ }
gsell's avatar
gsell committed
163 164 165 166 167 168 169 170 171 172 173 174 175
}

void RFCavity::accept(BeamlineVisitor &visitor) const {
    visitor.visitRFCavity(*this);
}

/**
 * ENVELOPE COMPONENT for radial focussing of the beam
 * Calculates the transverse envelope component for the RF cavity
 * element and adds it to the K vector
*/
void RFCavity::addKR(int i, double t, Vector_t &K) {

176 177
    double pz = RefPartBunch_m->getZ(i) - startField_m;
    const Vector_t tmpR(RefPartBunch_m->getX(i), RefPartBunch_m->getY(i), pz);
gsell's avatar
gsell committed
178 179
    double k = -Physics::q_e / (2.0 * RefPartBunch_m->getGamma(i) * Physics::EMASS);

180 181 182
    Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
    fieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
    double Ez = tmpE(2);
gsell's avatar
gsell committed
183

184 185
    tmpE = Vector_t(0.0);
    fieldmap_m->getFieldDerivative(tmpR, tmpE, tmpB, DZ);
gsell's avatar
gsell committed
186

187 188 189
    double wtf = frequency_m * t + phase_m;
    double kj = k * scale_m * (tmpE(2) * cos(wtf) - RefPartBunch_m->getBeta(i) * frequency_m * Ez * sin(wtf) / Physics::c);
    K += Vector_t(kj, kj, 0.0);
gsell's avatar
gsell committed
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206
}


/**
 * ENVELOPE COMPONENT for transverse kick (only has an impact if x0, y0 != 0)
 * Calculates the transverse kick component for the RF cavity element and adds it to
 * the K vector. Only important for off track tracking, otherwise KT = 0.
*/
void RFCavity::addKT(int i, double t, Vector_t &K) {

    RefPartBunch_m->actT();

    //XXX: BET parameter, default is 1.
    //If cxy != 1, then cxy = true
    bool cxy = false; // default
    double kx = 0.0, ky = 0.0;
    if(cxy) {
207 208 209 210 211 212 213 214 215 216
        double pz = RefPartBunch_m->getZ(i) - startField_m;
        const Vector_t tmpA(RefPartBunch_m->getX(i), RefPartBunch_m->getY(i), pz);

        Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
        fieldmap_m->getFieldstrength(tmpA, tmpE, tmpB);

        double cwtf = cos(frequency_m * t + phase_m);
        double cf = -Physics::q_e / (RefPartBunch_m->getGamma(i) * Physics::m_e);
        kx += -cf * scale_m * tmpE(0) * cwtf;
        ky += -cf * scale_m * tmpE(1) * cwtf;
gsell's avatar
gsell committed
217 218
    }

219 220
    double dx = RefPartBunch_m->getX0(i);
    double dy = RefPartBunch_m->getY0(i);
gsell's avatar
gsell committed
221 222 223 224 225 226 227 228 229 230

    Vector_t KR(0.0, 0.0, 0.0);
    addKR(i, t, KR);
    //FIXME ?? different in bet src
    K += Vector_t(KR(1) * dx + kx, KR(1) * dy + ky, 0.0);
    //
    //K += Vector_t(kx - KR(1) * dx, ky - KR(1) * dy, 0.0);
}


231 232 233
bool RFCavity::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
    return apply(RefPartBunch_m->R[i], RefPartBunch_m->P[i], t, E, B);
}
234

235 236 237 238 239 240 241
bool RFCavity::apply(const Vector_t &R,
                     const Vector_t &P,
                     const double &t,
                     Vector_t &E,
                     Vector_t &B) {
    Vector_t tmpR(R(0), R(1), R(2) - startField_m);
    Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
gsell's avatar
gsell committed
242

243 244 245 246
    if (tmpR(2) >= 0.0 &&
        tmpR(2) < length_m) {
        bool outOfBounds = fieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
        if (outOfBounds) return true;
gsell's avatar
gsell committed
247

248 249
        E += (scale_m + scaleError_m) * cos(frequency_m * t + phase_m + phaseError_m) * tmpE;
        B -= (scale_m + scaleError_m) * sin(frequency_m * t + phase_m + phaseError_m) * tmpB;
gsell's avatar
gsell committed
250 251

    }
252
    return false;
gsell's avatar
gsell committed
253 254
}

255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270
bool RFCavity::applyToReferenceParticle(const Vector_t &R,
                                        const Vector_t &P,
                                        const double &t,
                                        Vector_t &E,
                                        Vector_t &B) {

    Vector_t tmpR(R(0), R(1), R(2) - startField_m);
    Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);

    if (tmpR(2) >= 0.0 &&
        tmpR(2) < length_m) {
        bool outOfBounds = fieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
        if (outOfBounds) return true;

        E += scale_m * cos(frequency_m * t + phase_m) * tmpE;
        B -= scale_m * sin(frequency_m * t + phase_m) * tmpB;
gsell's avatar
gsell committed
271 272

    }
273
    return false;
gsell's avatar
gsell committed
274 275
}

276
void RFCavity::initialise(PartBunch *bunch, double &startField, double &endField) {
gsell's avatar
gsell committed
277
    using Physics::two_pi;
278 279 280 281 282 283 284 285 286

    if (bunch == NULL) {
        startField = startField_m;
        endField = endField_m;

        return;
    }

    double rBegin = 0.0, rEnd = 0.0;
kraus's avatar
kraus committed
287
    Inform msg("RFCavity ", *gmsg);
288
    std::stringstream errormsg;
gsell's avatar
gsell committed
289 290
    RefPartBunch_m = bunch;

291
    fieldmap_m = Fieldmap::getFieldmap(filename_m, fast_m);
gsell's avatar
gsell committed
292

293 294 295 296 297 298 299 300 301 302 303 304 305 306
    fieldmap_m->getFieldDimensions(startField_m, endField_m, rBegin, rEnd);
    if(endField_m > startField_m) {
        msg << level2 << getName() << " using file ";
        fieldmap_m->getInfo(&msg);
        if(std::abs((frequency_m - fieldmap_m->getFrequency()) / frequency_m) > 0.01) {
            errormsg << "FREQUENCY IN INPUT FILE DIFFERENT THAN IN FIELD MAP '" << filename_m << "';\n"
                     << frequency_m / two_pi * 1e-6 << " MHz <> "
                     << fieldmap_m->getFrequency() / two_pi * 1e-6 << " MHz; TAKE ON THE LATTER";
            std::string errormsg_str = Fieldmap::typeset_msg(errormsg.str(), "warning");
            ERRORMSG(errormsg_str << "\n" << endl);
            if(Ippl::myNode() == 0) {
                std::ofstream omsg("errormsg.txt", std::ios_base::app);
                omsg << errormsg_str << std::endl;
                omsg.close();
gsell's avatar
gsell committed
307
            }
308
            frequency_m = fieldmap_m->getFrequency();
gsell's avatar
gsell committed
309
        }
310 311
        length_m = endField_m - startField_m;
        endField = startField + length_m;
gsell's avatar
gsell committed
312 313 314 315 316 317
    } else {
        endField = startField - 1e-3;
    }
}

// In current version ,this function reads in the cavity voltage profile data from file.
318 319 320 321
void RFCavity::initialise(PartBunch *bunch,
                          std::shared_ptr<AbstractTimeDependence> freq_atd,
                          std::shared_ptr<AbstractTimeDependence> ampl_atd,
                          std::shared_ptr<AbstractTimeDependence> phase_atd) {
gsell's avatar
gsell committed
322 323 324 325
    using Physics::pi;

    RefPartBunch_m = bunch;

326 327 328 329 330
    /// set the time dependent models
    setAmplitudeModel(ampl_atd);
    setPhaseModel(phase_atd);
    setFrequencyModel(freq_atd);

gsell's avatar
gsell committed
331 332
    ifstream in(filename_m.c_str());
    if(!in.good()) {
333 334
        throw GeneralClassicException("RFCavity::initialise",
                                      "failed to open file '" + filename_m + "', please check if it exists");
gsell's avatar
gsell committed
335
    }
adelmann's avatar
Cleanup  
adelmann committed
336
    *gmsg << "* Read cavity voltage profile data" << endl;
gsell's avatar
gsell committed
337 338 339

    in >> num_points_m;

340 341 342
    RNormal_m  = std::unique_ptr<double[]>(new double[num_points_m]);
    VrNormal_m = std::unique_ptr<double[]>(new double[num_points_m]);
    DvDr_m     = std::unique_ptr<double[]>(new double[num_points_m]);
gsell's avatar
gsell committed
343 344 345

    for(int i = 0; i < num_points_m; i++) {
        if(in.eof()) {
346 347
            throw GeneralClassicException("RFCavity::initialise",
                                          "not enough data in file '" + filename_m + "', please check the data format");
gsell's avatar
gsell committed
348 349 350 351 352 353 354 355 356
        }
        in >> RNormal_m[i] >> VrNormal_m[i] >> DvDr_m[i];

        VrNormal_m[i] *= RefPartBunch_m->getQ();
        DvDr_m[i]     *= RefPartBunch_m->getQ();
    }
    sinAngle_m = sin(angle_m / 180.0 * pi);
    cosAngle_m = cos(angle_m / 180.0 * pi);

adelmann's avatar
cleanup  
adelmann committed
357 358
    if (frequency_name_m != "")
      *gmsg << "* Timedependent frequency model " << frequency_name_m << endl;
359

360
    *gmsg << "* Cavity voltage data read successfully!" << endl;
gsell's avatar
gsell committed
361 362 363 364 365 366 367 368 369 370
}

void RFCavity::finalise()
{}

bool RFCavity::bends() const {
    return false;
}


kraus's avatar
kraus committed
371
void RFCavity::goOnline(const double &) {
372 373
    Fieldmap::readMap(filename_m);

gsell's avatar
gsell committed
374 375 376 377
    online_m = true;
}

void RFCavity::goOffline() {
378 379
    Fieldmap::freeMap(filename_m);

gsell's avatar
gsell committed
380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438
    online_m = false;
}

void  RFCavity::setRmin(double rmin) {
    rmin_m = rmin;
}

void  RFCavity::setRmax(double rmax) {
    rmax_m = rmax;
}

void  RFCavity::setAzimuth(double angle) {
    angle_m = angle;
}

void  RFCavity::setPerpenDistance(double pdis) {
    pdis_m = pdis;
}

void  RFCavity::setGapWidth(double gapwidth) {
    gapwidth_m = gapwidth;
}

void RFCavity::setPhi0(double phi0) {
    phi0_m = phi0;
}

double  RFCavity::getRmin() const {
    return rmin_m;
}

double  RFCavity::getRmax() const {
    return rmax_m;
}

double  RFCavity::getAzimuth() const {
    return angle_m;
}

double  RFCavity::getSinAzimuth() const {
    return sinAngle_m;
}

double  RFCavity::getCosAzimuth() const {
    return cosAngle_m;
}

double  RFCavity::getPerpenDistance() const {
    return pdis_m;
}

double  RFCavity::getGapWidth() const {
    return gapwidth_m;
}

double RFCavity::getPhi0() const {
    return phi0_m;
}

439
void RFCavity::setComponentType(std::string name) {
440
    name = Util::toUpper(name);
gsell's avatar
gsell committed
441 442 443 444
    if(name == "STANDING") {
        type_m = SW;
    } else if(name == "SINGLEGAP") {
        type_m = SGSW;
445 446 447 448 449 450 451 452 453
    } else if(name != "") {
        std::stringstream errormsg;
        errormsg << getName() << ": CAVITY TYPE " << name << " DOES NOT EXIST;";
        std::string errormsg_str = Fieldmap::typeset_msg(errormsg.str(), "warning");
        ERRORMSG(errormsg_str << "\n" << endl);
        if(Ippl::myNode() == 0) {
            ofstream omsg("errormsg.txt", ios_base::app);
            omsg << errormsg_str << endl;
            omsg.close();
gsell's avatar
gsell committed
454
        }
455 456
        throw GeneralClassicException("RFCavity::setComponentType", errormsg_str);
    } else {
gsell's avatar
gsell committed
457 458 459 460 461 462 463
        type_m = SW;
    }

}

string RFCavity::getComponentType()const {
    if(type_m == SGSW)
464
        return std::string("SINGLEGAP");
gsell's avatar
gsell committed
465
    else
466
        return std::string("STANDING");
gsell's avatar
gsell committed
467 468 469 470 471 472
}

double RFCavity::getCycFrequency()const {
    return  frequency_m;
}

473 474 475
/**
   \brief used in OPAL-cycl

476
   Is called from OPAL-cycl and can handle
477 478 479 480 481 482
   time dependent frequency, amplitude and phase

   At the moment (test) only the frequence is time
   dependent

 */
483
void RFCavity::getMomentaKick(const double normalRadius, double momentum[], const double t, const double dtCorrt, const int PID, const double restMass, const int chargenumber) {
gsell's avatar
gsell committed
484 485 486 487 488 489 490 491 492 493 494 495 496 497 498
    using Physics::two_pi;
    using Physics::pi;
    using Physics::c;
    double derivate;
    double Voltage;

    double momentum2  = momentum[0] * momentum[0] + momentum[1] * momentum[1] + momentum[2] * momentum[2];
    double betgam = sqrt(momentum2);

    double gamma = sqrt(1.0 + momentum2);
    double beta = betgam / gamma;

    Voltage = spline(normalRadius, &derivate) * scale_m * 1.0e6; // V

    double transit_factor = 0.0;
499
    double Ufactor = 1.0;
gsell's avatar
gsell committed
500

501
    double frequency = frequency_m * frequency_td_m->getValue(t);
502

gsell's avatar
gsell committed
503
    if(gapwidth_m > 0.0) {
504
        transit_factor = 0.5 * frequency * gapwidth_m * 1.0e-3 / (c * beta);
505
        Ufactor = sin(transit_factor) / transit_factor;
gsell's avatar
gsell committed
506 507
    }

508
    Voltage *= Ufactor;
gsell's avatar
gsell committed
509

510
    double nphase = (frequency * (t + dtCorrt) * 1.0e-9) - phi0_m / 180.0 * pi ; // rad/s, ns --> rad
snuverink_j's avatar
snuverink_j committed
511
    double dgam = Voltage * cos(nphase) / (restMass);
gsell's avatar
gsell committed
512 513 514 515 516 517 518 519 520 521

    double tempdegree = fmod(nphase * 360.0 / two_pi, 360.0);
    if(tempdegree > 270.0) tempdegree -= 360.0;

    gamma += dgam;

    double newmomentum2 = pow(gamma, 2) - 1.0;

    double pr = momentum[0] * cosAngle_m + momentum[1] * sinAngle_m;
    double ptheta = sqrt(newmomentum2 - pow(pr, 2));
522 523
    double px = pr * cosAngle_m - ptheta * sinAngle_m ; // x
    double py = pr * sinAngle_m + ptheta * cosAngle_m; // y
gsell's avatar
gsell committed
524

525
    double rotate = -derivate * (scale_m * 1.0e6) / ((rmax_m - rmin_m) / 1000.0) * sin(nphase) / (frequency * two_pi) / (betgam * restMass / c / chargenumber); // radian
gsell's avatar
gsell committed
526 527

    /// B field effects
528 529
    momentum[0] =  cos(rotate) * px + sin(rotate) * py;
    momentum[1] = -sin(rotate) * px + cos(rotate) * py;
gsell's avatar
gsell committed
530 531

    if(PID == 0) {
532
      *gmsg << "* Cavity " << getName() << " Phase= " << tempdegree << " [deg] transit time factor=  " << Ufactor
533 534
            << " dE= " << dgam *restMass * 1.0e-6 << " [MeV]"
            << " E_kin= " << (gamma - 1.0)*restMass * 1.0e-6 << " [MeV] Time dep freq = " << frequency_td_m->getValue(t) << endl;
gsell's avatar
gsell committed
535
    }
536

537
}
gsell's avatar
gsell committed
538 539 540 541 542 543 544

/* cubic spline subrutine */
double RFCavity::spline(double z, double *za) {
    double splint;

    // domain-test and handling of case "1-support-point"
    if(num_points_m < 1) {
545 546
        throw GeneralClassicException("RFCavity::spline",
                                      "no support points!");
gsell's avatar
gsell committed
547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605
    }
    if(num_points_m == 1) {
        splint = RNormal_m[0];
        *za = 0.0;
        return splint;
    }

    // search the two support-points
    int il, ih;
    il = 0;
    ih = num_points_m - 1;
    while((ih - il) > 1) {
        int i = (int)((il + ih) / 2.0);
        if(z < RNormal_m[i]) {
            ih = i;
        } else if(z > RNormal_m[i]) {
            il = i;
        } else if(z == RNormal_m[i]) {
            il = i;
            ih = i + 1;
            break;
        }
    }

    double x1 =  RNormal_m[il];
    double x2 =  RNormal_m[ih];
    double y1 =  VrNormal_m[il];
    double y2 =  VrNormal_m[ih];
    double y1a = DvDr_m[il];
    double y2a = DvDr_m[ih];
    //
    // determination of the requested function-values and its derivatives
    //
    double dx  = x2 - x1;
    double dy  = y2 - y1;
    double u   = (z - x1) / dx;
    double u2  = u * u;
    double u3  = u2 * u;
    double dy2 = -2.0 * dy;
    double ya2 = y2a + 2.0 * y1a;
    double dy3 = 3.0 * dy;
    double ya3 = y2a + y1a;
    double yb2 = dy2 + dx * ya3;
    double yb4 = dy3 - dx * ya2;
    splint = y1  + u * dx * y1a +       u2 * yb4 +        u3 * yb2;
    *za    =            y1a + 2.0 * u / dx * yb4 + 3.0 * u2 / dx * yb2;
    // if(m>=1) za=y1a+2.0*u/dx*yb4+3.0*u2/dx*yb2;
    // if(m>=2) za[1]=2.0/dx2*yb4+6.0*u/dx2*yb2;
    // if(m>=3) za[2]=6.0/dx3*yb2;

    return splint;
}

void RFCavity::getDimensions(double &zBegin, double &zEnd) const {
    zBegin = startField_m;
    zEnd = endField_m;
}


606 607
ElementBase::ElementType RFCavity::getType() const {
    return RFCAVITY;
gsell's avatar
gsell committed
608 609 610 611
}

double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const double &q, const double &mass) {
    vector<double> t, E, t2, E2;
612
    std::vector< double > F;
gsell's avatar
gsell committed
613
    std::vector< std::pair< double, double > > G;
614 615
    gsl_spline *onAxisInterpolants;
    gsl_interp_accel *onAxisAccel;
gsell's avatar
gsell committed
616 617 618 619

    unsigned int N;
    double A, B;
    double phi = 0.0, tmp_phi, dphi = 0.5 * Physics::pi / 180.;
620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 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
    double dz = 1.0, length = 0.0;
    fieldmap_m->getOnaxisEz(G);
    double begin = (G.front()).first;
    double end = (G.back()).first;
    std::unique_ptr<double[]> zvals(new double[G.size()]);
    std::unique_ptr<double[]> onAxisField(new double[G.size()]);

    for(size_t j = 0; j < G.size(); ++ j) {
        zvals[j] = G[j].first;
        onAxisField[j] = G[j].second;
    }
    onAxisInterpolants = gsl_spline_alloc(gsl_interp_cspline, G.size());
    onAxisAccel = gsl_interp_accel_alloc();
    gsl_spline_init(onAxisInterpolants, zvals.get(), onAxisField.get(), G.size());

    length = end - begin;
    dz = length / G.size();

    G.clear();

    N = (int)floor(length / dz + 1);
    dz = length / N;

    F.resize(N);
    double z = begin;
    for(size_t j = 0; j < N; ++ j, z += dz) {
        F[j] = gsl_spline_eval(onAxisInterpolants, z, onAxisAccel);
    }
    gsl_spline_free(onAxisInterpolants);
    gsl_interp_accel_free(onAxisAccel);

    t.resize(N, t0);
    t2.resize(N, t0);
    E.resize(N, E0);
    E2.resize(N, E0);

    z = begin + dz;
    for(unsigned int i = 1; i < N; ++ i, z += dz) {
        E[i] = E[i - 1] + dz * scale_m / mass;
        E2[i] = E[i];
    }

    for(int iter = 0; iter < 10; ++ iter) {
        A = B = 0.0;
        for(unsigned int i = 1; i < N; ++ i) {
            t[i] = t[i - 1] + getdT(i, E, dz, mass);
            t2[i] = t2[i - 1] + getdT(i, E2, dz, mass);
            A += scale_m * (1. + frequency_m * (t2[i] - t[i]) / dphi) * getdA(i, t, dz, frequency_m, F);
            B += scale_m * (1. + frequency_m * (t2[i] - t[i]) / dphi) * getdB(i, t, dz, frequency_m, F);
gsell's avatar
gsell committed
669
        }
670

671 672 673 674
        if(std::abs(B) > 0.0000001) {
            tmp_phi = atan(A / B);
        } else {
            tmp_phi = Physics::pi / 2;
gsell's avatar
gsell committed
675
        }
676 677
        if(q * (A * sin(tmp_phi) + B * cos(tmp_phi)) < 0) {
            tmp_phi += Physics::pi;
gsell's avatar
gsell committed
678 679
        }

680
        if(std::abs(phi - tmp_phi) < frequency_m * (t[N - 1] - t[0]) / (10 * N)) {
gsell's avatar
gsell committed
681
            for(unsigned int i = 1; i < N; ++ i) {
682 683
                E[i] = E[i - 1];
                E[i] += q * scale_m * getdE(i, t, dz, phi, frequency_m, F) ;
gsell's avatar
gsell committed
684
            }
685 686 687
            const int prevPrecision = Ippl::Info->precision(8);
            INFOMSG(level2 << "estimated phase= " << tmp_phi << " rad, "
                    << "Ekin= " << E[N - 1] << " MeV" << setprecision(prevPrecision) << endl);
gsell's avatar
gsell committed
688

689 690 691
            return tmp_phi;
        }
        phi = tmp_phi - floor(tmp_phi / Physics::two_pi + 0.5) * Physics::two_pi;
gsell's avatar
gsell committed
692 693


694 695 696 697 698
        for(unsigned int i = 1; i < N; ++ i) {
            E[i] = E[i - 1];
            E2[i] = E2[i - 1];
            E[i] += q * scale_m * getdE(i, t, dz, phi, frequency_m, F) ;
            E2[i] += q * scale_m * getdE(i, t2, dz, phi + dphi, frequency_m, F);
gsell's avatar
gsell committed
699

700 701
            t[i] = t[i - 1] + getdT(i, E, dz, mass);
            t2[i] = t2[i - 1] + getdT(i, E2, dz, mass);
702

703 704 705 706 707
            E[i] = E[i - 1];
            E2[i] = E2[i - 1];
            E[i] += q * scale_m * getdE(i, t, dz, phi, frequency_m, F) ;
            E2[i] += q * scale_m * getdE(i, t2, dz, phi + dphi, frequency_m, F);
        }
gsell's avatar
gsell committed
708

snuverink_j's avatar
snuverink_j committed
709
        double cosine_part = 0.0, sine_part = 0.0;
710 711 712 713
        double p0 = sqrt((E0 / mass + 1) * (E0 / mass + 1) - 1);
        cosine_part += scale_m * cos(frequency_m * t0) * F[0];
        sine_part += scale_m * sin(frequency_m * t0) * F[0];

snuverink_j's avatar
snuverink_j committed
714
        double totalEz0 = cos(phi) * cosine_part - sin(phi) * sine_part;
715 716 717 718 719 720 721 722

        if(p0 + q * totalEz0 * (t[1] - t[0]) * Physics::c / mass < 0) {
            // make totalEz0 = 0
            tmp_phi = atan(cosine_part / sine_part);
            if(abs(tmp_phi - phi) > Physics::pi) {
                phi = tmp_phi + Physics::pi;
            } else {
                phi = tmp_phi;
gsell's avatar
gsell committed
723 724
            }
        }
725
    }
gsell's avatar
gsell committed
726

727 728 729
    const int prevPrecision = Ippl::Info->precision(8);
    INFOMSG(level2 << "estimated phase= " << tmp_phi << " rad, "
            << "Ekin= " << E[N - 1] << " MeV" << setprecision(prevPrecision) << endl);
kraus's avatar
kraus committed
730

731
    return phi;
gsell's avatar
gsell committed
732 733 734
}

pair<double, double> RFCavity::trackOnAxisParticle(const double &p0,
735 736 737 738
                                                   const double &t0,
                                                   const double &dt,
                                                   const double &q,
                                                   const double &mass) {
739
    Vector_t p(0, 0, p0);
gsell's avatar
gsell committed
740
    double t = t0;
741
    BorisPusher integrator(*RefPartBunch_m->getReference());
742 743 744
    const double cdt = Physics::c * dt;
    const double zbegin = startField_m;
    const double zend = length_m + startField_m;
gsell's avatar
gsell committed
745

746
    Vector_t z(0.0, 0.0, zbegin);
747 748
    double dz = 0.5 * p(2) / sqrt(1.0 + dot(p, p)) * cdt;
    Vector_t Ef(0.0), Bf(0.0);
gsell's avatar
gsell committed
749

750 751 752 753 754 755 756
    while(z(2) + dz < zend && z(2) + dz > zbegin) {
        z /= cdt;
        integrator.push(z, p, dt);
        z *= cdt;

        if(z(2) >= zbegin && z(2) <= zend) {
            Ef = 0.0;
757
            Bf = 0.0;
758
            applyToReferenceParticle(z, p, t + 0.5 * dt, Ef, Bf);
gsell's avatar
gsell committed
759
        }
760 761 762 763 764 765
        integrator.kick(z, p, Ef, Bf, dt);

        dz = 0.5 * p(2) / sqrt(1.0 + dot(p, p)) * cdt;
        z /= cdt;
        integrator.push(z, p, dt);
        z *= cdt;
gsell's avatar
gsell committed
766 767 768
        t += dt;
    }

769 770
    const double beta = sqrt(1. - 1 / (dot(p, p) + 1.));
    const double tErr  = (z(2) - zend) / (Physics::c * beta);
gsell's avatar
gsell committed
771

772 773
    return pair<double, double>(p(2), t - tErr);
}
gsell's avatar
gsell committed
774

775 776 777 778 779 780
bool RFCavity::isInside(const Vector_t &r) const {
    if (isInsideTransverse(r)) {
        return fieldmap_m->isInside(r);
    }

    return false;
781
}