RFCavity.cpp 22.7 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/PartBunchBase.h"
24
#include "Steppers/BorisPusher.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"
31 32

#include <cmath>
gsell's avatar
gsell committed
33 34
#include <iostream>
#include <fstream>
35

gsell's avatar
gsell committed
36 37 38 39 40 41
extern Inform *gmsg;

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

RFCavity::RFCavity():
42 43
    RFCavity("")
{}
gsell's avatar
gsell committed
44 45 46 47


RFCavity::RFCavity(const RFCavity &right):
    Component(right),
48 49 50 51 52 53
    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
54 55
    filename_m(right.filename_m),
    scale_m(right.scale_m),
56
    scaleError_m(right.scaleError_m),
gsell's avatar
gsell committed
57
    phase_m(right.phase_m),
58
    phaseError_m(right.phaseError_m),
gsell's avatar
gsell committed
59
    frequency_m(right.frequency_m),
60 61 62
    fast_m(right.fast_m),
    autophaseVeto_m(right.autophaseVeto_m),
    designEnergy_m(right.designEnergy_m),
63
    fieldmap_m(right.fieldmap_m),
gsell's avatar
gsell committed
64
    startField_m(right.startField_m),
kraus's avatar
kraus committed
65
    endField_m(right.endField_m),
gsell's avatar
gsell committed
66 67 68 69 70 71 72 73 74
    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),
75 76 77
    RNormal_m(nullptr),
    VrNormal_m(nullptr),
    DvDr_m(nullptr),
78
    num_points_m(right.num_points_m)
snuverink_j's avatar
snuverink_j committed
79
{}
gsell's avatar
gsell committed
80 81


82
RFCavity::RFCavity(const std::string &name):
gsell's avatar
gsell committed
83
    Component(name),
84 85 86
    phase_td_m(nullptr),
    amplitude_td_m(nullptr),
    frequency_td_m(nullptr),
gsell's avatar
gsell committed
87 88
    filename_m(""),
    scale_m(1.0),
89
    scaleError_m(0.0),
gsell's avatar
gsell committed
90
    phase_m(0.0),
91
    phaseError_m(0.0),
gsell's avatar
gsell committed
92
    frequency_m(0.0),
93 94 95
    fast_m(true),
    autophaseVeto_m(false),
    designEnergy_m(-1.0),
96
    fieldmap_m(nullptr),
gsell's avatar
gsell committed
97
    startField_m(0.0),
kraus's avatar
kraus committed
98
    endField_m(0.0),
gsell's avatar
gsell committed
99 100 101 102 103 104 105 106 107
    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),
108 109 110
    RNormal_m(nullptr),
    VrNormal_m(nullptr),
    DvDr_m(nullptr),
111
    num_points_m(0)
snuverink_j's avatar
snuverink_j committed
112
{}
gsell's avatar
gsell committed
113 114 115 116 117 118 119 120 121


RFCavity::~RFCavity() {
}

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

122 123 124
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);
}
125

126
bool RFCavity::apply(const Vector_t &R,
127
                     const Vector_t &/*P*/,
128 129 130
                     const double &t,
                     Vector_t &E,
                     Vector_t &B) {
131
    if (R(2) >= startField_m &&
132
        R(2) < startField_m + getElementLength()) {
133
        Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
gsell's avatar
gsell committed
134

135
        bool outOfBounds = fieldmap_m->getFieldstrength(R, tmpE, tmpB);
136
        if (outOfBounds) return true;
gsell's avatar
gsell committed
137

138 139
        E += (scale_m + scaleError_m) * std::cos(frequency_m * t + phase_m + phaseError_m) * tmpE;
        B -= (scale_m + scaleError_m) * std::sin(frequency_m * t + phase_m + phaseError_m) * tmpB;
gsell's avatar
gsell committed
140 141

    }
142
    return false;
gsell's avatar
gsell committed
143 144
}

145
bool RFCavity::applyToReferenceParticle(const Vector_t &R,
146
                                        const Vector_t &/*P*/,
147 148 149 150
                                        const double &t,
                                        Vector_t &E,
                                        Vector_t &B) {

151
    if (R(2) >= startField_m &&
152
        R(2) < startField_m + getElementLength()) {
153
        Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
154

155
        bool outOfBounds = fieldmap_m->getFieldstrength(R, tmpE, tmpB);
156 157
        if (outOfBounds) return true;

158 159
        E += scale_m * std::cos(frequency_m * t + phase_m) * tmpE;
        B -= scale_m * std::sin(frequency_m * t + phase_m) * tmpB;
gsell's avatar
gsell committed
160 161

    }
162
    return false;
gsell's avatar
gsell committed
163 164
}

165
void RFCavity::initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField) {
166

kraus's avatar
kraus committed
167
    startField_m = endField_m = 0.0;
168 169
    if (bunch == NULL) {
        startField = startField_m;
kraus's avatar
kraus committed
170
        endField = endField_m;
171 172 173 174

        return;
    }

kraus's avatar
kraus committed
175
    Inform msg("RFCavity ", *gmsg);
176
    std::stringstream errormsg;
gsell's avatar
gsell committed
177 178
    RefPartBunch_m = bunch;

179
    fieldmap_m = Fieldmap::getFieldmap(filename_m, fast_m);
180
    fieldmap_m->getFieldDimensions(startField_m, endField);
181
    if (endField <= startField_m) {
182
        throw GeneralClassicException("RFCavity::initialise",
183
                                      "The length of the field map '" + filename_m + "' is zero or negative");
184 185 186 187
    }

    msg << level2 << getName() << " using file ";
    fieldmap_m->getInfo(&msg);
188
    if (std::abs((frequency_m - fieldmap_m->getFrequency()) / frequency_m) > 0.01) {
189
        errormsg << "FREQUENCY IN INPUT FILE DIFFERENT THAN IN FIELD MAP '" << filename_m << "';\n"
190 191
                 << frequency_m / Physics::two_pi * 1e-6 << " MHz <> "
                 << fieldmap_m->getFrequency() / Physics::two_pi * 1e-6 << " MHz; TAKE ON THE LATTER";
192 193
        std::string errormsg_str = Fieldmap::typeset_msg(errormsg.str(), "warning");
        ERRORMSG(errormsg_str << "\n" << endl);
194
        if (Ippl::myNode() == 0) {
195 196 197
            std::ofstream omsg("errormsg.txt", std::ios_base::app);
            omsg << errormsg_str << std::endl;
            omsg.close();
gsell's avatar
gsell committed
198
        }
199
        frequency_m = fieldmap_m->getFrequency();
gsell's avatar
gsell committed
200
    }
201
    setElementLength(endField - startField_m);
gsell's avatar
gsell committed
202 203 204
}

// In current version ,this function reads in the cavity voltage profile data from file.
205
void RFCavity::initialise(PartBunchBase<double, 3> *bunch,
206 207 208
                          std::shared_ptr<AbstractTimeDependence> freq_atd,
                          std::shared_ptr<AbstractTimeDependence> ampl_atd,
                          std::shared_ptr<AbstractTimeDependence> phase_atd) {
gsell's avatar
gsell committed
209 210 211

    RefPartBunch_m = bunch;

212 213 214 215 216
    /// set the time dependent models
    setAmplitudeModel(ampl_atd);
    setPhaseModel(phase_atd);
    setFrequencyModel(freq_atd);

217
    std::ifstream in(filename_m.c_str());
218
    if (!in.good()) {
219 220
        throw GeneralClassicException("RFCavity::initialise",
                                      "failed to open file '" + filename_m + "', please check if it exists");
gsell's avatar
gsell committed
221
    }
adelmann's avatar
Cleanup  
adelmann committed
222
    *gmsg << "* Read cavity voltage profile data" << endl;
gsell's avatar
gsell committed
223 224 225

    in >> num_points_m;

226 227 228
    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
229

230 231
    for (int i = 0; i < num_points_m; i++) {
        if (in.eof()) {
232 233
            throw GeneralClassicException("RFCavity::initialise",
                                          "not enough data in file '" + filename_m + "', please check the data format");
gsell's avatar
gsell committed
234 235 236 237 238 239
        }
        in >> RNormal_m[i] >> VrNormal_m[i] >> DvDr_m[i];

        VrNormal_m[i] *= RefPartBunch_m->getQ();
        DvDr_m[i]     *= RefPartBunch_m->getQ();
    }
240 241
    sinAngle_m = std::sin(angle_m * Physics::deg2rad);
    cosAngle_m = std::cos(angle_m * Physics::deg2rad);
gsell's avatar
gsell committed
242

adelmann's avatar
cleanup  
adelmann committed
243 244
    if (frequency_name_m != "")
      *gmsg << "* Timedependent frequency model " << frequency_name_m << endl;
245

246
    *gmsg << "* Cavity voltage data read successfully!" << endl;
gsell's avatar
gsell committed
247 248 249 250 251 252 253 254 255 256
}

void RFCavity::finalise()
{}

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


kraus's avatar
kraus committed
257
void RFCavity::goOnline(const double &) {
258 259
    Fieldmap::readMap(filename_m);

gsell's avatar
gsell committed
260 261 262 263
    online_m = true;
}

void RFCavity::goOffline() {
264 265
    Fieldmap::freeMap(filename_m);

gsell's avatar
gsell committed
266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324
    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;
}

325
void RFCavity::setComponentType(std::string name) {
326
    if (name == "STANDING") {
gsell's avatar
gsell committed
327
        type_m = SW;
328
    } else if (name == "SINGLEGAP") {
gsell's avatar
gsell committed
329
        type_m = SGSW;
330
    } else if (name != "") {
331 332 333 334
        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);
335
        if (Ippl::myNode() == 0) {
336 337
            std::ofstream omsg("errormsg.txt", std::ios_base::app);
            omsg << errormsg_str << std::endl;
338
            omsg.close();
gsell's avatar
gsell committed
339
        }
340 341
        throw GeneralClassicException("RFCavity::setComponentType", errormsg_str);
    } else {
gsell's avatar
gsell committed
342 343 344 345 346
        type_m = SW;
    }

}

347
std::string RFCavity::getComponentType()const {
348
    if (type_m == SGSW)
349
        return std::string("SINGLEGAP");
gsell's avatar
gsell committed
350
    else
351
        return std::string("STANDING");
gsell's avatar
gsell committed
352 353 354 355 356 357
}

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

358 359 360
/**
   \brief used in OPAL-cycl

361
   Is called from OPAL-cycl and can handle
362 363
   time dependent frequency, amplitude and phase

364
   At the moment (test) only the frequency is time
365 366 367
   dependent

 */
368
void RFCavity::getMomentaKick(const double normalRadius, double momentum[], const double t, const double dtCorrt, const int PID, const double restMass, const int chargenumber) {
369

gsell's avatar
gsell committed
370 371 372
    double derivate;

    double momentum2  = momentum[0] * momentum[0] + momentum[1] * momentum[1] + momentum[2] * momentum[2];
373
    double betgam = std::sqrt(momentum2);
gsell's avatar
gsell committed
374

375
    double gamma = std::sqrt(1.0 + momentum2);
gsell's avatar
gsell committed
376 377
    double beta = betgam / gamma;

378
    double Voltage = spline(normalRadius, &derivate) * scale_m * 1.0e6; // V
gsell's avatar
gsell committed
379

380
    double Ufactor = 1.0;
gsell's avatar
gsell committed
381

382
    double frequency = frequency_m * frequency_td_m->getValue(t);
383

384
    if (gapwidth_m > 0.0) {
385
    	double transit_factor = 0.5 * frequency * gapwidth_m * 1.0e-3 / (Physics::c * beta);
386
        Ufactor = std::sin(transit_factor) / transit_factor;
gsell's avatar
gsell committed
387 388
    }

389
    Voltage *= Ufactor;
390 391 392
    // rad/s, ns --> rad
    double nphase = (frequency * (t + dtCorrt) * 1.0e-9) - phi0_m / 180.0 * Physics::pi ;
    double dgam = Voltage * std::cos(nphase) / (restMass);
gsell's avatar
gsell committed
393

394
    double tempdegree = std::fmod(nphase * 360.0 / Physics::two_pi, 360.0);
395
    if (tempdegree > 270.0) tempdegree -= 360.0;
gsell's avatar
gsell committed
396 397 398

    gamma += dgam;

399
    double newmomentum2 = std::pow(gamma, 2) - 1.0;
gsell's avatar
gsell committed
400 401

    double pr = momentum[0] * cosAngle_m + momentum[1] * sinAngle_m;
402
    double ptheta = std::sqrt(newmomentum2 - std::pow(pr, 2));
403 404
    double px = pr * cosAngle_m - ptheta * sinAngle_m ; // x
    double py = pr * sinAngle_m + ptheta * cosAngle_m; // y
gsell's avatar
gsell committed
405

406
    double rotate = -derivate * (scale_m * 1.0e6) / ((rmax_m - rmin_m) / 1000.0) * std::sin(nphase) / (frequency * Physics::two_pi) / (betgam * restMass / Physics::c / chargenumber); // radian
gsell's avatar
gsell committed
407 408

    /// B field effects
409 410
    momentum[0] =  std::cos(rotate) * px + std::sin(rotate) * py;
    momentum[1] = -std::sin(rotate) * px + std::cos(rotate) * py;
gsell's avatar
gsell committed
411

412
    if (PID == 0) {
413

kraus's avatar
kraus committed
414
        Inform  m("OPAL", *gmsg, Ippl::myNode());
415

kraus's avatar
kraus committed
416 417 418
        m << "* Cavity " << getName() << " Phase= " << tempdegree << " [deg] transit time factor=  " << Ufactor
          << " 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
419
    }
420

421
}
gsell's avatar
gsell committed
422 423 424 425 426 427

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

    // domain-test and handling of case "1-support-point"
428
    if (num_points_m < 1) {
429 430
        throw GeneralClassicException("RFCavity::spline",
                                      "no support points!");
gsell's avatar
gsell committed
431
    }
432
    if (num_points_m == 1) {
gsell's avatar
gsell committed
433 434 435 436 437 438 439 440 441
        splint = RNormal_m[0];
        *za = 0.0;
        return splint;
    }

    // search the two support-points
    int il, ih;
    il = 0;
    ih = num_points_m - 1;
442
    while ((ih - il) > 1) {
gsell's avatar
gsell committed
443
        int i = (int)((il + ih) / 2.0);
444
        if (z < RNormal_m[i]) {
gsell's avatar
gsell committed
445
            ih = i;
446
        } else if (z > RNormal_m[i]) {
gsell's avatar
gsell committed
447
            il = i;
448
        } else if (z == RNormal_m[i]) {
gsell's avatar
gsell committed
449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485
            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;
kraus's avatar
kraus committed
486
    zEnd = endField_m;
gsell's avatar
gsell committed
487 488 489
}


490 491
ElementBase::ElementType RFCavity::getType() const {
    return RFCAVITY;
gsell's avatar
gsell committed
492 493
}

494
double RFCavity::getAutoPhaseEstimateFallback(double E0, double t0, double q, double mass) {
495

496
    const double dt = 1e-13;
497
    const double p0 = Util::getBetaGamma(E0, mass);
498
    const double origPhase =getPhasem();
499
    double dphi = Physics::pi / 18;
500 501 502

    double phi = 0.0;
    setPhasem(phi);
503
    std::pair<double, double> ret = trackOnAxisParticle(p0, t0, dt, q, mass);
504
    double phimax = 0.0;
505
    double Emax = Util::getKineticEnergy(Vector_t(0.0, 0.0, ret.first), mass);
506 507
    phi += dphi;

508 509 510 511
    for (unsigned int j = 0; j < 2; ++ j) {
        for (unsigned int i = 0; i < 36; ++ i, phi += dphi) {
            setPhasem(phi);
            ret = trackOnAxisParticle(p0, t0, dt, q, mass);
512
            double Ekin = Util::getKineticEnergy(Vector_t(0.0, 0.0, ret.first), mass);
513 514 515 516
            if (Ekin > Emax) {
                Emax = Ekin;
                phimax = phi;
            }
517
        }
518 519 520

        phi = phimax - dphi;
        dphi = dphi / 17.5;
521 522
    }

523
    phimax = phimax - std::round(phimax / Physics::two_pi) * Physics::two_pi;
524
    phimax = std::fmod(phimax, Physics::two_pi);
525

526
    const int prevPrecision = Ippl::Info->precision(8);
kraus's avatar
kraus committed
527
    INFOMSG(level2
528
            << "estimated phase= " << phimax << " rad = "
529
            << phimax * Physics::rad2deg << " deg \n"
530
            << "Ekin= " << Emax << " MeV" << std::setprecision(prevPrecision) << "\n" << endl);
531

532
    setPhasem(origPhase);
533 534 535
    return phimax;
}

gsell's avatar
gsell committed
536
double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const double &q, const double &mass) {
537
    std::vector<double> t, E, t2, E2;
538
    std::vector<double> F;
gsell's avatar
gsell committed
539
    std::vector< std::pair< double, double > > G;
540 541
    gsl_spline *onAxisInterpolants;
    gsl_interp_accel *onAxisAccel;
gsell's avatar
gsell committed
542 543

    double phi = 0.0, tmp_phi, dphi = 0.5 * Physics::pi / 180.;
544 545
    double dz = 1.0, length = 0.0;
    fieldmap_m->getOnaxisEz(G);
546
    if (G.size() == 0) return 0.0;
547
    double begin = (G.front()).first;
548 549
    double end   = (G.back()).first;
    std::unique_ptr<double[]> zvals(      new double[G.size()]);
550 551
    std::unique_ptr<double[]> onAxisField(new double[G.size()]);

552
    for (size_t j = 0; j < G.size(); ++ j) {
553 554 555 556 557 558 559 560 561 562 563 564
        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();

565
    unsigned int N = (int)floor(length / dz + 1);
566 567 568 569
    dz = length / N;

    F.resize(N);
    double z = begin;
570
    for (size_t j = 0; j < N; ++ j, z += dz) {
571 572 573 574 575 576 577 578 579 580 581
        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;
582
    for (unsigned int i = 1; i < N; ++ i, z += dz) {
583 584 585 586
        E[i] = E[i - 1] + dz * scale_m / mass;
        E2[i] = E[i];
    }

587
    for (int iter = 0; iter < 10; ++ iter) {
588 589
        double A = 0.0;
        double B = 0.0;
590
        for (unsigned int i = 1; i < N; ++ i) {
591 592 593 594
            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
595
        }
596

597
        if (std::abs(B) > 0.0000001) {
598 599 600
            tmp_phi = atan(A / B);
        } else {
            tmp_phi = Physics::pi / 2;
gsell's avatar
gsell committed
601
        }
602
        if (q * (A * sin(tmp_phi) + B * cos(tmp_phi)) < 0) {
603
            tmp_phi += Physics::pi;
gsell's avatar
gsell committed
604 605
        }

606 607
        if (std::abs (phi - tmp_phi) < frequency_m * (t[N - 1] - t[0]) / (10 * N)) {
            for (unsigned int i = 1; i < N; ++ i) {
608 609
                E[i] = E[i - 1];
                E[i] += q * scale_m * getdE(i, t, dz, phi, frequency_m, F) ;
gsell's avatar
gsell committed
610
            }
611
            const int prevPrecision = Ippl::Info->precision(8);
kraus's avatar
kraus committed
612
            INFOMSG(level2 << "estimated phase= " << tmp_phi << " rad = "
613
                    << tmp_phi * Physics::rad2deg << " deg \n"
614
                    << "Ekin= " << E[N - 1] << " MeV" << std::setprecision(prevPrecision) << "\n" << endl);
gsell's avatar
gsell committed
615

616 617
            return tmp_phi;
        }
618
        phi = tmp_phi - std::round(tmp_phi / Physics::two_pi) * Physics::two_pi;
gsell's avatar
gsell committed
619

620
        for (unsigned int i = 1; i < N; ++ i) {
621 622 623 624
            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);
625 626
            double a = E[i], b = E2[i];
            if (std::isnan(a) || std::isnan(b)) {
627
                return getAutoPhaseEstimateFallback(E0, t0, q, mass);
628
            }
629 630
            t[i] = t[i - 1] + getdT(i, E, dz, mass);
            t2[i] = t2[i - 1] + getdT(i, E2, dz, mass);
631

632
            E[i]  = E [i - 1];
633
            E2[i] = E2[i - 1];
634
            E[i]  += q * scale_m * getdE(i, t, dz, phi, frequency_m, F) ;
635 636
            E2[i] += q * scale_m * getdE(i, t2, dz, phi + dphi, frequency_m, F);
        }
gsell's avatar
gsell committed
637

snuverink_j's avatar
snuverink_j committed
638
        double cosine_part = 0.0, sine_part = 0.0;
639
        double p0 = Util::getBetaGamma(E0, mass);
640 641
        cosine_part += scale_m * std::cos(frequency_m * t0) * F[0];
        sine_part += scale_m * std::sin(frequency_m * t0) * F[0];
642

643
        double totalEz0 = std::cos(phi) * cosine_part - sin(phi) * sine_part;
644

645
        if (p0 + q * totalEz0 * (t[1] - t[0]) * Physics::c / mass < 0) {
646
            // make totalEz0 = 0
647
            tmp_phi = std::atan(cosine_part / sine_part);
648
            if (std::abs (tmp_phi - phi) > Physics::pi) {
649 650 651
                phi = tmp_phi + Physics::pi;
            } else {
                phi = tmp_phi;
gsell's avatar
gsell committed
652 653
            }
        }
654
    }
gsell's avatar
gsell committed
655

656
    const int prevPrecision = Ippl::Info->precision(8);
kraus's avatar
kraus committed
657
    INFOMSG(level2
658
            << "estimated phase= " << tmp_phi << " rad = "
659
            << tmp_phi * Physics::rad2deg << " deg \n"
660
            << "Ekin= " << E[N - 1] << " MeV" << std::setprecision(prevPrecision) << "\n" << endl);
kraus's avatar
kraus committed
661

662
    return phi;
gsell's avatar
gsell committed
663 664
}

665
std::pair<double, double> RFCavity::trackOnAxisParticle(const double &p0,
666 667 668 669 670
                                                        const double &t0,
                                                        const double &dt,
                                                        const double &/*q*/,
                                                        const double &mass,
                                                        std::ofstream *out) {
671
    Vector_t p(0, 0, p0);
gsell's avatar
gsell committed
672
    double t = t0;
673
    BorisPusher integrator(*RefPartBunch_m->getReference());
674 675
    const double cdt = Physics::c * dt;
    const double zbegin = startField_m;
676
    const double zend = getElementLength() + startField_m;
gsell's avatar
gsell committed
677

678
    Vector_t z(0.0, 0.0, zbegin);
679
    double dz = 0.5 * p(2) / Util::getGamma(p) * cdt;
680
    Vector_t Ef(0.0), Bf(0.0);
gsell's avatar
gsell committed
681

682
    if (out) *out << std::setw(18) << z[2]
683
                  << std::setw(18) << Util::getKineticEnergy(p, mass)
kraus's avatar
kraus committed
684
                  << std::endl;
685
    while (z(2) + dz < zend && z(2) + dz > zbegin) {
686 687 688 689
        z /= cdt;
        integrator.push(z, p, dt);
        z *= cdt;

690 691
        Ef = 0.0;
        Bf = 0.0;
692
        if (z(2) >= zbegin && z(2) <= zend) {
693
            applyToReferenceParticle(z, p, t + 0.5 * dt, Ef, Bf);
gsell's avatar
gsell committed
694
        }
695 696
        integrator.kick(z, p, Ef, Bf, dt);

697
        dz = 0.5 * p(2) / std::sqrt(1.0 + dot(p, p)) * cdt;
698 699 700
        z /= cdt;
        integrator.push(z, p, dt);
        z *= cdt;
gsell's avatar
gsell committed
701
        t += dt;
702

kraus's avatar
kraus committed
703
        if (out) *out << std::setw(18) << z[2]
704
                      << std::setw(18) << Util::getKineticEnergy(p, mass)
kraus's avatar
kraus committed
705
                      << std::endl;
gsell's avatar
gsell committed
706 707
    }

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

711
    return std::pair<double, double>(p(2), t - tErr);
712
}
gsell's avatar
gsell committed
713

714 715 716 717 718 719
bool RFCavity::isInside(const Vector_t &r) const {
    if (isInsideTransverse(r)) {
        return fieldmap_m->isInside(r);
    }

    return false;
720 721 722 723 724
}

double RFCavity::getElementLength() const {
    double length = ElementBase::getElementLength();
    if (length < 1e-10 && fieldmap_m != NULL) {
gsell's avatar
gsell committed
725
        double start, end;
726
        fieldmap_m->getFieldDimensions(start, end);
727 728 729 730 731 732 733 734
        length = end - start;
    }

    return length;
}

void RFCavity::getElementDimensions(double &begin,
                                    double &end) const {
735
    fieldmap_m->getFieldDimensions(begin, end);
gsell's avatar
gsell committed
736 737 738 739 740 741 742 743
}

// vi: set et ts=4 sw=4 sts=4:
// Local Variables:
// mode:c++
// c-basic-offset: 4
// indent-tabs-mode: nil
// require-final-newline: nil
744
// End: