RFCavity.cpp 34.5 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 "Fields/Fieldmap.h"
25
#include "Utilities/GeneralClassicException.h"
gsell's avatar
gsell committed
26 27 28 29 30 31

#include "gsl/gsl_interp.h"
#include "gsl/gsl_spline.h"
#include <iostream>
#include <fstream>

32 33 34 35
#ifdef OPAL_NOCPLUSPLUS11_NULLPTR
#define nullptr NULL
#endif

gsell's avatar
gsell committed
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
extern Inform *gmsg;

using namespace std;

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

RFCavity::RFCavity():
    Component(),
    filename_m(""),
    scale_m(1.0),
    phase_m(0.0),
    frequency_m(0.0),
    ElementEdge_m(0.0),
    startField_m(0.0),
    endField_m(0.0),
    type_m(SW),
    fast_m(false),
    autophaseVeto_m(false),
    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),
63 64 65
    RNormal_m(nullptr),
    VrNormal_m(nullptr),
    DvDr_m(nullptr),
66 67 68 69 70 71 72 73
    num_points_m(0),
    phase_td_m(nullptr),
    amplitude_td_m(nullptr),
    frequency_td_m(nullptr),
    phase_name_m(""),
    amplitude_name_m(""),
    frequency_name_m("")
{
gsell's avatar
gsell committed
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
    setElType(isRF);
}


RFCavity::RFCavity(const RFCavity &right):
    Component(right),
    filename_m(right.filename_m),
    scale_m(right.scale_m),
    phase_m(right.phase_m),
    frequency_m(right.frequency_m),
    ElementEdge_m(right.ElementEdge_m),
    startField_m(right.startField_m),
    endField_m(right.endField_m),
    type_m(right.type_m),
    fast_m(right.fast_m),
    autophaseVeto_m(right.autophaseVeto_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),
98 99 100
    RNormal_m(nullptr),
    VrNormal_m(nullptr),
    DvDr_m(nullptr),
101 102 103 104 105 106 107 108
    num_points_m(right.num_points_m),
    phase_td_m(right.phase_td_m),
    amplitude_td_m(right.amplitude_td_m),
    frequency_td_m(right.frequency_td_m),
    phase_name_m(right.phase_name_m),
    amplitude_name_m(right.amplitude_name_m),
    frequency_name_m(right.frequency_name_m) {

gsell's avatar
gsell committed
109
    setElType(isRF);
110

gsell's avatar
gsell committed
111
    std::vector<string>::const_iterator fname_it;
112
    for(fname_it = right.multiFilenames_m.begin(); fname_it != right.multiFilenames_m.end(); ++ fname_it) {
gsell's avatar
gsell committed
113 114
        multiFilenames_m.push_back(*fname_it);
    }
115 116
    std::vector<Fieldmap *>::const_iterator fmap_it;
    for(fmap_it = right.multiFieldmaps_m.begin(); fmap_it != right.multiFieldmaps_m.end(); ++ fmap_it) {
gsell's avatar
gsell committed
117 118 119
        multiFieldmaps_m.push_back(*fmap_it);
    }
    std::vector<double>::const_iterator scale_it;
120
    for(scale_it = right.multiScales_m.begin(); scale_it != right.multiScales_m.end(); ++ scale_it) {
gsell's avatar
gsell committed
121 122 123
        multiScales_m.push_back(*scale_it);
    }
    std::vector<double>::const_iterator phase_it;
124
    for(phase_it = right.multiPhases_m.begin(); phase_it != right.multiPhases_m.end(); ++ phase_it) {
gsell's avatar
gsell committed
125 126 127
        multiPhases_m.push_back(*phase_it);
    }
    std::vector<double>::const_iterator freq_it;
128
    for(freq_it = right.multiFrequencies_m.begin(); freq_it != right.multiFrequencies_m.end(); ++ freq_it) {
gsell's avatar
gsell committed
129 130 131 132 133
        multiFrequencies_m.push_back(*freq_it);
    }
}


134
RFCavity::RFCavity(const std::string &name):
gsell's avatar
gsell committed
135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153
    Component(name),
    filename_m(""),
    scale_m(1.0),
    phase_m(0.0),
    frequency_m(0.0),
    ElementEdge_m(0.0),
    startField_m(0.0),
    endField_m(0.0),
    type_m(SW),
    fast_m(false),
    autophaseVeto_m(false),
    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),
154 155 156
    RNormal_m(nullptr),
    VrNormal_m(nullptr),
    DvDr_m(nullptr),
157 158 159
    phase_td_m(nullptr),
    amplitude_td_m(nullptr),
    frequency_td_m(nullptr),
160 161 162
    //     RNormal_m(std::nullptr_t(NULL)),
    //     VrNormal_m(std::nullptr_t(NULL)),
    //     DvDr_m(std::nullptr_t(NULL)),
gsell's avatar
gsell committed
163 164 165 166 167 168 169 170
    num_points_m(0) {
    setElType(isRF);
}


RFCavity::~RFCavity() {
    // FIXME: in deleteFielmak, a map find makes problems
    //       Fieldmap::deleteFieldmap(filename_m);
171
    //~ if(RNormal_m) {
172 173 174
    //~ delete[] RNormal_m;
    //~ delete[] VrNormal_m;
    //~ delete[] DvDr_m;
175
    //~ }
gsell's avatar
gsell committed
176 177 178
}


179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211
std::shared_ptr<AbstractTimeDependence> RFCavity::getAmplitudeModel() const {
  return amplitude_td_m;
}

std::shared_ptr<AbstractTimeDependence> RFCavity::getPhaseModel() const {
  return phase_td_m;
}

std::shared_ptr<AbstractTimeDependence> RFCavity::getFrequencyModel() const {
  return frequency_td_m;
}

void RFCavity::setAmplitudeModel(std::shared_ptr<AbstractTimeDependence> amplitude_td) {
  amplitude_td_m = amplitude_td;
}

void RFCavity::setPhaseModel(std::shared_ptr<AbstractTimeDependence> phase_td) {
  phase_td_m = phase_td;
}

void RFCavity::setFrequencyModel(std::shared_ptr<AbstractTimeDependence> frequency_td) {
  INFOMSG("frequency_td " << frequency_td << endl);
  frequency_td_m = frequency_td;
}


void RFCavity::setFrequencyModelName(std::string name) {
  frequency_name_m=name;
}




gsell's avatar
gsell committed
212 213 214 215 216
void RFCavity::accept(BeamlineVisitor &visitor) const {
    visitor.visitRFCavity(*this);
}

void RFCavity::dropFieldmaps() {
217
    std::vector<Fieldmap *>::iterator fmap_it;
gsell's avatar
gsell committed
218 219 220 221 222 223 224 225 226 227
    for(fmap_it = multiFieldmaps_m.begin(); fmap_it != multiFieldmaps_m.end(); ++ fmap_it) {
        *fmap_it = NULL;
    }
    multiFilenames_m.clear();
    multiFieldmaps_m.clear();
    multiScales_m.clear();
    multiPhases_m.clear();
    multiFrequencies_m.clear();
}

228
void RFCavity::setFieldMapFN(std::string fn) {
gsell's avatar
gsell committed
229 230 231 232 233
    multiFilenames_m.push_back(fn);
    filename_m = multiFilenames_m[0];
}

string RFCavity::getFieldMapFN() const {
234
    std::string filename("No_fieldmap");
235
    if(numFieldmaps() > 0) {
gsell's avatar
gsell committed
236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265
        filename = multiFilenames_m[0];
    }
    return filename;
}

void RFCavity::setAmplitudem(double vPeak) {
    multiScales_m.push_back(vPeak);
    scale_m = multiScales_m[0];
}

void RFCavity::setFrequency(double freq) {
    frequency_m = freq;
}

void RFCavity::setFrequencym(double freq) {
    multiFrequencies_m.push_back(freq);
    frequency_m = multiFrequencies_m[0];
}


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

void RFCavity::setPhasem(double phase) {
    multiPhases_m.push_back(phase);
    phase_m = multiPhases_m[0];
}

void RFCavity::updatePhasem(double phase) {
266
    if(multiPhases_m.size() == 0) {
gsell's avatar
gsell committed
267 268 269 270
        multiPhases_m.push_back(phase);
    } else {
        double diff = phase - multiPhases_m[0];
        multiPhases_m[0] = phase;
271
        for(size_t i = 1; i < numFieldmaps(); ++ i) {
gsell's avatar
gsell committed
272 273 274 275 276 277 278 279 280 281
            multiPhases_m[i] += diff;
        }
    }
    phase_m = multiPhases_m[0];
}

double RFCavity::getPhasem() const {
    return phase_m;
}

282
void RFCavity::setCavityType(std::string type) {
gsell's avatar
gsell committed
283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300

}

string RFCavity::getCavityType() const {
    return "SW";
}

void RFCavity::setFast(bool fast) {
    fast_m = fast;
}


bool RFCavity::getFast() const {
    return fast_m;
}



301
void RFCavity::setAutophaseVeto(bool veto) {
gsell's avatar
gsell committed
302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321
    autophaseVeto_m = veto;
}


bool RFCavity::getAutophaseVeto() const {
    return autophaseVeto_m;
}


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

    double pz = RefPartBunch_m->getZ(i) - startField_m - ds_m;
    const Vector_t tmpR(RefPartBunch_m->getX(i) - dx_m, RefPartBunch_m->getY(i) - dy_m, pz);
    double k = -Physics::q_e / (2.0 * RefPartBunch_m->getGamma(i) * Physics::EMASS);

322 323
    for(size_t j = 0; j < numFieldmaps(); ++ j) {
        Fieldmap *fieldmap = multiFieldmaps_m[j];
gsell's avatar
gsell committed
324 325 326 327 328 329 330 331 332
        double frequency = multiFrequencies_m[j];
        double scale = multiScales_m[j];
        double phase = multiPhases_m[j];

        Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
        fieldmap->getFieldstrength(tmpR, tmpE, tmpB);
        double Ez = tmpE(2);

        tmpE = Vector_t(0.0);
333
        fieldmap->getFieldDerivative(tmpR, tmpE, tmpB, DZ);
gsell's avatar
gsell committed
334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355

        double wtf = frequency * t + phase;
        double kj = k * scale * (tmpE(2) * cos(wtf) - RefPartBunch_m->getBeta(i) * frequency * Ez * sin(wtf) / Physics::c);
        K += Vector_t(kj, kj, 0.0);
    }
}


/**
 * 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) {
356 357
        for(size_t j = 0; j < numFieldmaps(); ++ j) {
            Fieldmap *fieldmap = multiFieldmaps_m[j];
gsell's avatar
gsell committed
358 359 360 361 362 363 364 365 366
            double scale = multiScales_m[j];
            double frequency = multiFrequencies_m[j];
            double phase = multiPhases_m[j];

            double pz = RefPartBunch_m->getZ(i) - startField_m - ds_m;
            const Vector_t tmpA(RefPartBunch_m->getX(i) - dx_m, RefPartBunch_m->getY(i) - dy_m, pz);

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

gsell's avatar
gsell committed
368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386
            double cwtf = cos(frequency * t + phase);
            double cf = -Physics::q_e / (RefPartBunch_m->getGamma(i) * Physics::m_e);
            kx += -cf * scale * tmpE(0) * cwtf;
            ky += -cf * scale * tmpE(1) * cwtf;
        }
    }

    double dx = RefPartBunch_m->getX0(i) - dx_m;
    double dy = RefPartBunch_m->getY0(i) - dy_m;

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


387
bool RFCavity::apply(const size_t &i, const double &t, double E[], double B[]) {
gsell's avatar
gsell committed
388 389
    Vector_t Ev(0, 0, 0), Bv(0, 0, 0);
    Vector_t Rt(RefPartBunch_m->getX(i), RefPartBunch_m->getY(i), RefPartBunch_m->getZ(i));
390

gsell's avatar
gsell committed
391 392 393 394 395 396 397 398 399 400 401 402
    if(apply(Rt, Vector_t(0.0), t, Ev, Bv)) return true;

    E[0] = Ev(0);
    E[1] = Ev(1);
    E[2] = Ev(2);
    B[0] = Bv(0);
    B[1] = Bv(1);
    B[2] = Bv(2);

    return false;
}

403
bool RFCavity::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
gsell's avatar
gsell committed
404 405
    bool out_of_bounds = true;
    const Vector_t tmpR(RefPartBunch_m->getX(i) - dx_m, RefPartBunch_m->getY(i) - dy_m , RefPartBunch_m->getZ(i) - startField_m - ds_m);
406

407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422
    if (tmpR(2) >= 0.0) {
        // inside the cavity
        for(size_t j = 0; j < numFieldmaps(); ++ j) {
            Fieldmap *fieldmap = multiFieldmaps_m[j];
            Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
            const std::pair<double, double> & start_end = multi_start_end_field_m[j];

            if(tmpR(2) > start_end.first &&
               tmpR(2) < start_end.second &&
               !fieldmap->getFieldstrength(tmpR, tmpE, tmpB)) {
                const double phase = multiFrequencies_m[j] * t + multiPhases_m[j];
                const double &scale = multiScales_m[j];
                E += scale * cos(phase) * tmpE;
                B -= scale * sin(phase) * tmpB;
                out_of_bounds = false;
            }
gsell's avatar
gsell committed
423 424
        }
    }
425
    else {
426
        /*
427 428 429 430 431
           some of the bunch is still outside of the cavity
           so let them drift in
        */
        out_of_bounds = false;
    }
gsell's avatar
gsell committed
432 433 434 435 436 437 438
    return out_of_bounds;
}

bool RFCavity::apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B) {
    bool out_of_bounds = true;
    const Vector_t tmpR(R(0) - dx_m, R(1) - dy_m, R(2) - startField_m - ds_m);

439 440
    for(size_t i = 0; i < numFieldmaps(); ++ i) {
        Fieldmap *fieldmap = multiFieldmaps_m[i];
gsell's avatar
gsell committed
441 442 443
        Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
        const std::pair<double, double> & start_end = multi_start_end_field_m[i];

444
        if(tmpR(2) > start_end.first &&
gsell's avatar
gsell committed
445 446 447
           tmpR(2) < start_end.second &&
           !fieldmap->getFieldstrength(tmpR, tmpE, tmpB)) {
            const double phase = multiFrequencies_m[i] * t + multiPhases_m[i];
448
            const double &scale = multiScales_m[i];
gsell's avatar
gsell committed
449 450 451 452 453 454 455 456 457 458 459
            E += scale * cos(phase) * tmpE;
            B -= scale * sin(phase) * tmpB;
            out_of_bounds = false;
        }
    }
    return out_of_bounds;
}

void RFCavity::initialise(PartBunch *bunch, double &startField, double &endField, const double &scaleFactor) {
    using Physics::two_pi;
    double zBegin = 0.0, zEnd = 0.0, rBegin = 0.0, rEnd = 0.0;
kraus's avatar
kraus committed
460
    Inform msg("RFCavity ", *gmsg);
461
    std::stringstream errormsg;
gsell's avatar
gsell committed
462 463
    RefPartBunch_m = bunch;

464
    for(size_t i = 0; i < numFieldmaps(); ++ i) {
465
        std::string fname = multiFilenames_m[i];
466 467
        Fieldmap *fmap = Fieldmap::getFieldmap(fname, fast_m);
        if(fmap != NULL) {
gsell's avatar
gsell committed
468 469
            multiFieldmaps_m.push_back(fmap);
        } else {
470
            for(size_t j = i; j < numFieldmaps() - 1; ++ j) {
gsell's avatar
gsell committed
471 472 473 474 475 476
                multiFieldmaps_m[j] = multiFieldmaps_m[j + 1];
            }
            multiFieldmaps_m.pop_back();
            -- i;
        }
    }
477 478
    if(multiScales_m.size() > 0) {
        while(multiScales_m.size() < numFieldmaps()) {
gsell's avatar
gsell committed
479 480 481
            multiScales_m.push_back(multiScales_m[0]);
        }
    } else {
482
        while(multiScales_m.size() < numFieldmaps()) {
gsell's avatar
gsell committed
483 484 485
            multiScales_m.push_back(1.0);
        }
    }
486 487
    if(multiFrequencies_m.size() > 0) {
        while(multiFrequencies_m.size() < numFieldmaps()) {
gsell's avatar
gsell committed
488 489 490
            multiFrequencies_m.push_back(multiFrequencies_m[0]);
        }
    } else {
491
        while(multiFrequencies_m.size() < numFieldmaps()) {
gsell's avatar
gsell committed
492 493 494
            multiFrequencies_m.push_back(0.0);
        }
    }
495 496
    if(multiPhases_m.size() > 0) {
        while(multiPhases_m.size() < numFieldmaps()) {
gsell's avatar
gsell committed
497 498 499
            multiPhases_m.push_back(multiPhases_m[0]);
        }
    } else {
500
        while(multiPhases_m.size() < multiFilenames_m.size()) {
gsell's avatar
gsell committed
501 502 503 504
            multiPhases_m.push_back(0.0);
        }
    }

505
    if(numFieldmaps() > 0) {
gsell's avatar
gsell committed
506
        double overall_zBegin = 999999999.99, overall_zEnd = -999999999.99;
507 508 509
        for(size_t i = 0; i < numFieldmaps(); ++ i) {
            Fieldmap *fmap = multiFieldmaps_m[i];
            double &frequency = multiFrequencies_m[i];
510
            const std::string &filename = multiFilenames_m[i];
gsell's avatar
gsell committed
511 512 513

            fmap->getFieldDimensions(zBegin, zEnd, rBegin, rEnd);
            if(zEnd > zBegin) {
kraus's avatar
kraus committed
514 515
                msg << level2 << getName() << " using file ";
                fmap->getInfo(&msg);
gsell's avatar
gsell committed
516 517 518 519
                if(fabs((frequency - fmap->getFrequency()) / frequency) > 0.01) {
                    errormsg << "FREQUENCY IN INPUT FILE DIFFERENT THAN IN FIELD MAP '" << filename << "';\n"
                             << frequency / two_pi * 1e-6 << " MHz <> "
                             << fmap->getFrequency() / two_pi * 1e-6 << " MHz; TAKE ON THE LATTER";
520
                    std::string errormsg_str = Fieldmap::typeset_msg(errormsg.str(), "warning");
kraus's avatar
kraus committed
521
                    ERRORMSG(errormsg_str << "\n" << endl);
gsell's avatar
gsell committed
522 523 524 525 526 527 528 529 530 531 532 533
                    if(Ippl::myNode() == 0) {
                        std::ofstream omsg("errormsg.txt", std::ios_base::app);
                        omsg << errormsg_str << std::endl;
                        omsg.close();
                    }
                    frequency = fmap->getFrequency();
                }
                multi_start_end_field_m.push_back(std::pair<double, double>(zBegin, zEnd));

                overall_zBegin = std::min(overall_zBegin, zBegin);
                overall_zEnd = std::max(overall_zEnd, zEnd);
            } else {
534
                for(size_t j = i; j < numFieldmaps() - 1; ++ j) {
gsell's avatar
gsell committed
535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551
                    multiFilenames_m[j] = multiFilenames_m[j + 1];
                    multiFieldmaps_m[j] = multiFieldmaps_m[j + 1];
                    multiScales_m[j] = multiScales_m[j + 1];
                    multiPhases_m[j] = multiPhases_m[j + 1];
                    multiFrequencies_m[j] = multiFrequencies_m[j + 1];
                }
                multiFilenames_m.pop_back();
                multiFieldmaps_m.back() = NULL;
                multiFieldmaps_m.pop_back();
                multiScales_m.pop_back();
                multiPhases_m.pop_back();
                multiFrequencies_m.pop_back();
                -- i;
            }
        }
        zBegin = overall_zBegin;
        zEnd = overall_zEnd;
552
        for(size_t i = 0; i < numFieldmaps(); ++ i) {
gsell's avatar
gsell committed
553 554 555 556 557
            multi_start_end_field_m[i].first -= zBegin;
            multi_start_end_field_m[i].second -= zBegin;
        }
    }

558
    if(zEnd > zBegin) {
gsell's avatar
gsell committed
559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574
        ElementEdge_m = startField;
        startField_m = startField = ElementEdge_m + zBegin;
        endField_m = endField = ElementEdge_m + zEnd;
    } else {
        endField = startField - 1e-3;
    }
}

// In current version ,this function reads in the cavity voltage profile data from file.
void RFCavity::initialise(PartBunch *bunch, const double &scaleFactor) {
    using Physics::pi;

    RefPartBunch_m = bunch;

    ifstream in(filename_m.c_str());
    if(!in.good()) {
575 576
        throw GeneralClassicException("RFCavity::initialise",
                                      "failed to open file '" + filename_m + "', please check if it exists");
gsell's avatar
gsell committed
577
    }
adelmann's avatar
Cleanup  
adelmann committed
578
    *gmsg << "* Read cavity voltage profile data" << endl;
gsell's avatar
gsell committed
579 580 581

    in >> num_points_m;

582 583 584
    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
585 586 587

    for(int i = 0; i < num_points_m; i++) {
        if(in.eof()) {
588 589
            throw GeneralClassicException("RFCavity::initialise",
                                          "not enough data in file '" + filename_m + "', please check the data format");
gsell's avatar
gsell committed
590 591 592 593 594 595 596 597 598
        }
        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);

599 600 601
    if (frequency_td_m)
      *gmsg << "* timedependent frequency model " << frequency_name_m << endl;
    
602
    *gmsg << "* Cavity voltage data read successfully!" << endl;
gsell's avatar
gsell committed
603 604 605 606 607 608 609 610 611 612
}

void RFCavity::finalise()
{}

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


kraus's avatar
kraus committed
613
void RFCavity::goOnline(const double &) {
gsell's avatar
gsell committed
614
    std::vector<string>::iterator fmap_it;
615
    for(fmap_it = multiFilenames_m.begin(); fmap_it != multiFilenames_m.end(); ++ fmap_it) {
gsell's avatar
gsell committed
616 617 618 619 620 621 622
        Fieldmap::readMap(*fmap_it);
    }
    online_m = true;
}

void RFCavity::goOffline() {
    std::vector<string>::iterator fmap_it;
623
    for(fmap_it = multiFilenames_m.begin(); fmap_it != multiFilenames_m.end(); ++ fmap_it) {
gsell's avatar
gsell committed
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 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684
        Fieldmap::freeMap(*fmap_it);
    }
    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;
}

685
void RFCavity::setComponentType(std::string name) {
gsell's avatar
gsell committed
686 687 688 689 690 691
    if(name == "STANDING") {
        type_m = SW;
    } else if(name == "SINGLEGAP") {
        type_m = SGSW;
    } else {
        if(name != "") {
692
            std::stringstream errormsg;
gsell's avatar
gsell committed
693 694
            errormsg << "CAVITY TYPE " << name << " DOES NOT EXIST; \n"
                     << "CHANGING TO REGULAR STANDING WAVE";
695
            std::string errormsg_str = Fieldmap::typeset_msg(errormsg.str(), "warning");
kraus's avatar
kraus committed
696
            ERRORMSG(errormsg_str << "\n" << endl);
gsell's avatar
gsell committed
697 698 699 700 701 702 703 704 705 706 707 708 709
            if(Ippl::myNode() == 0) {
                ofstream omsg("errormsg.txt", ios_base::app);
                omsg << errormsg_str << endl;
                omsg.close();
            }
        }
        type_m = SW;
    }

}

string RFCavity::getComponentType()const {
    if(type_m == SGSW)
710
        return std::string("SINGLEGAP");
gsell's avatar
gsell committed
711
    else
712
        return std::string("STANDING");
gsell's avatar
gsell committed
713 714 715 716 717 718
}

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

719
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
720 721 722 723 724 725 726 727 728 729 730 731 732 733 734
    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;
735
    double Ufactor = 1.0;
gsell's avatar
gsell committed
736 737 738

    if(gapwidth_m > 0.0) {
        transit_factor = 0.5 * frequency_m * gapwidth_m * 1.0e-3 / (c * beta);
739
        Ufactor = sin(transit_factor) / transit_factor;
gsell's avatar
gsell committed
740 741
    }

742
    Voltage *= Ufactor;
gsell's avatar
gsell committed
743 744 745 746 747 748 749 750 751 752 753 754 755 756 757

    double dgam = 0.0;
    double nphase = (frequency_m * (t + dtCorrt) * 1.0e-9) - phi0_m / 180.0 * pi ; // rad/s, ns --> rad

    dgam = Voltage * cos(nphase) / (restMass);

    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));
758 759
    double px = pr * cosAngle_m - ptheta * sinAngle_m ; // x
    double py = pr * sinAngle_m + ptheta * cosAngle_m; // y
gsell's avatar
gsell committed
760

761
    double rotate = -derivate * (scale_m * 1.0e6) / ((rmax_m - rmin_m) / 1000.0) * sin(nphase) / (frequency_m * two_pi) / (betgam * restMass / c / chargenumber); // radian
gsell's avatar
gsell committed
762 763

    /// B field effects
764 765
    momentum[0] =  cos(rotate) * px + sin(rotate) * py;
    momentum[1] = -sin(rotate) * px + cos(rotate) * py;
gsell's avatar
gsell committed
766 767

    if(PID == 0) {
768
      *gmsg << "* Cavity " << getName() << " Phase= " << tempdegree << " [deg] transit time factor=  " << Ufactor
adelmann's avatar
Cleanup  
adelmann committed
769 770
	      << " dE= " << dgam *restMass * 1.0e-6 << " [MeV]"
	      << " E_kin= " << (gamma - 1.0)*restMass * 1.0e-6 << " [MeV]" << endl;
gsell's avatar
gsell committed
771 772
    }

773
}
gsell's avatar
gsell committed
774 775 776 777 778 779 780

/* 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) {
781 782
        throw GeneralClassicException("RFCavity::spline",
                                      "no support points!");
gsell's avatar
gsell committed
783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841
    }
    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;
}


842 843
ElementBase::ElementType RFCavity::getType() const {
    return RFCAVITY;
gsell's avatar
gsell committed
844 845 846 847 848 849 850 851
}

double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const double &q, const double &mass) {
    vector<double> t, E, t2, E2;
    std::vector< std::vector< double > > F(numFieldmaps());
    std::vector< std::pair< double, double > > G;
    std::vector< double > begin(numFieldmaps());
    std::vector< double > end(numFieldmaps());
852 853
    std::vector<gsl_spline *> onAxisInterpolants(numFieldmaps());
    std::vector<gsl_interp_accel *> onAxisAccel(numFieldmaps());
gsell's avatar
gsell committed
854 855 856 857 858 859

    unsigned int N;
    double A, B;
    double phi = 0.0, tmp_phi, dphi = 0.5 * Physics::pi / 180.;
    double min_dz = 1.0, max_length = 0.0, min_begin = 99999999.99;
    double max_frequency = 0.0;
860
    for(size_t i = 0; i < numFieldmaps(); ++ i) {
gsell's avatar
gsell committed
861 862 863 864
        multiFieldmaps_m[i]->getOnaxisEz(G);
        begin[i] = (G.front()).first;
        end[i] = (G.back()).first;
        max_frequency = std::max(max_frequency, multiFrequencies_m[i]);
865 866
        std::unique_ptr<double[]> zvals(new double[G.size()]);
        std::unique_ptr<double[]> onAxisField(new double[G.size()]);
gsell's avatar
gsell committed
867

868
        for(size_t j = 0; j < G.size(); ++ j) {
gsell's avatar
gsell committed
869 870 871 872 873
            zvals[j] = G[j].first;
            onAxisField[j] = G[j].second;
        }
        onAxisInterpolants[i] = gsl_spline_alloc(gsl_interp_cspline, G.size());
        onAxisAccel[i] = gsl_interp_accel_alloc();
874
        gsl_spline_init(onAxisInterpolants[i], zvals.get(), onAxisField.get(), G.size());
gsell's avatar
gsell committed
875 876 877 878 879 880 881

        double length = end[i] - begin[i];
        min_dz = std::min(min_dz, length / G.size());
        max_length = std::max(max_length, length);
        min_begin = std::min(min_begin, begin[i]);

        G.clear();
882 883
        //~ delete[] zvals;
        //~ delete[] onAxisField;
gsell's avatar
gsell committed
884
    }
885

gsell's avatar
gsell committed
886 887 888
    N = (int)floor(max_length / min_dz + 1);
    min_dz = max_length / N;

889
    for(size_t i = 0; i < numFieldmaps(); ++ i) {
gsell's avatar
gsell committed
890 891
        F[i].resize(N);
        double z = min_begin;
892 893
        for(size_t j = 0; j < N; ++ j, z += min_dz) {
            if(z >= begin[i] && z <= end[i]) {
gsell's avatar
gsell committed
894 895 896 897 898 899 900 901
                F[i][j] = gsl_spline_eval(onAxisInterpolants[i], z, onAxisAccel[i]);
            } else {
                F[i][j] = 0.0;
            }
        }
        gsl_spline_free(onAxisInterpolants[i]);
        gsl_interp_accel_free(onAxisAccel[i]);
    }
902 903 904

    //~ delete[] onAxisInterpolants;
    //~ delete[] onAxisAccel;
gsell's avatar
gsell committed
905 906 907 908 909 910 911 912 913 914

    if(N > 0) {

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

        double z = min_begin + min_dz;
        for(unsigned int i = 1; i < N; ++ i, z += min_dz) {
915
            E[i] = E[i - 1] + min_dz * scale_m / mass;
gsell's avatar
gsell committed
916 917 918 919 920 921
            E2[i] = E[i];
        }

        for(int iter = 0; iter < 10; ++ iter) {
            A = B = 0.0;
            for(unsigned int i = 1; i < N; ++ i) {
922 923 924 925 926
                t[i] = t[i - 1] + getdT(i, E, min_dz, mass);
                t2[i] = t2[i - 1] + getdT(i, E2, min_dz, mass);
                for(size_t j = 0; j < numFieldmaps(); ++ j) {
                    const double &frequency = multiFrequencies_m[j];
                    const double &scale = multiScales_m[j];
gsell's avatar
gsell committed
927 928 929 930 931 932 933 934 935 936 937 938 939 940 941
                    const double Dphi = multiPhases_m[j] - multiPhases_m[0];
                    A += scale * (1. + frequency * (t2[i] - t[i]) / dphi) * getdA(i, t, min_dz, Dphi, frequency, F[j]);
                    B += scale * (1. + frequency * (t2[i] - t[i]) / dphi) * getdB(i, t, min_dz, Dphi, frequency, F[j]);
                }
            }

            if(fabs(B) > 0.0000001) {
                tmp_phi = atan(A / B);
            } else {
                tmp_phi = Physics::pi / 2;
            }
            if(q * (A * sin(tmp_phi) + B * cos(tmp_phi)) < 0) {
                tmp_phi += Physics::pi;
            }

942
            if(fabs(phi - tmp_phi) < max_frequency * (t[N - 1] - t[0]) / (10 * N)) {
gsell's avatar
gsell committed
943
                for(unsigned int i = 1; i < N; ++ i) {
944 945 946 947
                    E[i] = E[i - 1];
                    for(size_t j = 0; j < numFieldmaps(); ++ j) {
                        const double &frequency = multiFrequencies_m[j];
                        const double &scale = multiScales_m[j];
gsell's avatar
gsell committed
948 949 950 951
                        const double Dphi = multiPhases_m[j] - multiPhases_m[0];
                        E[i] += q * scale * getdE(i, t, min_dz, phi + Dphi, frequency, F[j]) ;
                    }
                }
kraus's avatar
kraus committed
952
                const int prevPrecision = Ippl::Info->precision(8);
953
                INFOMSG(level2 << "estimated phase= " << tmp_phi << " rad, "
kraus's avatar
kraus committed
954 955
                        << "Ekin= " << E[N - 1] << " MeV" << setprecision(prevPrecision) << endl);

gsell's avatar
gsell committed
956 957 958 959 960 961
                return tmp_phi;
            }
            phi = tmp_phi - floor(tmp_phi / Physics::two_pi + 0.5) * Physics::two_pi;


            for(unsigned int i = 1; i < N; ++ i) {
962 963 964 965 966
                E[i] = E[i - 1];
                E2[i] = E2[i - 1];
                for(size_t j = 0; j < numFieldmaps(); ++ j) {
                    const double &frequency = multiFrequencies_m[j];
                    const double &scale = multiScales_m[j];
gsell's avatar
gsell committed
967 968 969 970
                    const double Dphi = multiPhases_m[j] - multiPhases_m[0];
                    E[i] += q * scale * getdE(i, t, min_dz, phi + Dphi, frequency, F[j]) ;
                    E2[i] += q * scale * getdE(i, t2, min_dz, phi + dphi + Dphi, frequency, F[j]);
                }
971 972 973 974 975 976 977 978
                t[i] = t[i - 1] + getdT(i, E, min_dz, mass);
                t2[i] = t2[i - 1] + getdT(i, E2, min_dz, mass);

                E[i] = E[i - 1];
                E2[i] = E2[i - 1];
                for(size_t j = 0; j < numFieldmaps(); ++ j) {
                    const double &frequency = multiFrequencies_m[j];
                    const double &scale = multiScales_m[j];
gsell's avatar
gsell committed
979 980 981 982 983 984 985 986
                    const double Dphi = multiPhases_m[j] - multiPhases_m[0];
                    E[i] += q * scale * getdE(i, t, min_dz, phi + Dphi, frequency, F[j]) ;
                    E2[i] += q * scale * getdE(i, t2, min_dz, phi + dphi + Dphi, frequency, F[j]);
                }
            }

            double totalEz0 = 0.0, cosine_part = 0.0, sine_part = 0.0;
            double p0 = sqrt((E0 / mass + 1) * (E0 / mass + 1) - 1);
987 988 989
            for(size_t j = 0; j < numFieldmaps(); ++ j) {
                const double &frequency = multiFrequencies_m[j];
                const double &scale = multiScales_m[j];
gsell's avatar
gsell committed
990 991 992 993 994 995
                const double Dphi = multiPhases_m[j] - multiPhases_m[0];
                cosine_part += scale * cos(frequency * t0 + Dphi) * F[j][0];
                sine_part += scale * sin(frequency * t0 + Dphi) * F[j][0];
            }
            totalEz0 = cos(phi) * cosine_part - sin(phi) * sine_part;

996
            if(p0 + q * totalEz0 * (t[1] - t[0]) * Physics::c / mass < 0) {
gsell's avatar
gsell committed
997 998
                // make totalEz0 = 0
                tmp_phi = atan(cosine_part / sine_part);
999
                if(abs(tmp_phi - phi) > Physics::pi) {
gsell's avatar
gsell committed
1000 1001 1002 1003 1004 1005 1006
                    phi = tmp_phi + Physics::pi;
                } else {
                    phi = tmp_phi;
                }
            }
        }

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

gsell's avatar
gsell committed
1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026
        return phi;
    } else {
        return 0.0;
    }
}

pair<double, double> RFCavity::trackOnAxisParticle(const double &p0,
        const double &t0,
        const double &dt,
        const double &q,
        const double &mass) {
    double p = p0;
    double t = t0;
    double cdt = Physics::c * dt;
    vector<pair<double, double> > F;
    vector<size_t> length_fieldmaps(numFieldmaps());
1027
    for(size_t i = 0; i < numFieldmaps(); ++ i) {
gsell's avatar
gsell committed
1028 1029 1030 1031 1032 1033
        vector<pair<double, double> > G;
        multiFieldmaps_m[i]->getOnaxisEz(G);
        length_fieldmaps[i] = G.size();
        F.insert(F.end(), G.begin(), G.end());
    }

1034 1035
    std::unique_ptr<double[]> zvals(new double[F.size()]);
    std::unique_ptr<double[]> onAxisField(new double[F.size()]);
gsell's avatar
gsell committed
1036 1037 1038 1039 1040 1041 1042
    for(unsigned int i = 0; i < F.size(); ++i) {
        zvals[i] = F[i].first;
        onAxisField[i] = F[i].second;
    }
    double zbegin = zvals[0];
    double zend = zvals[F.size() - 1];
    std::vector<std::pair<double, double> > begin_end(numFieldmaps());
1043 1044
    std::vector<gsl_spline *> onAxisInterpolants(numFieldmaps());
    std::vector<gsl_interp_accel *> onAxisAccel(numFieldmaps());
gsell's avatar
gsell committed
1045
    size_t start = 0;
1046
    for(size_t i = 0; i < numFieldmaps(); ++ i) {
gsell's avatar
gsell committed
1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057
        begin_end[i].first = zvals[start];
        onAxisInterpolants[i] = gsl_spline_alloc(gsl_interp_cspline, length_fieldmaps[i]);
        onAxisAccel[i] = gsl_interp_accel_alloc();
        gsl_spline_init(onAxisInterpolants[i], &(zvals[start]), &(onAxisField[start]), length_fieldmaps[i]);
        start += length_fieldmaps[i];
        begin_end[i].second = zvals[start - 1];

        zbegin = std::min(zbegin, begin_end[i].first);
        zend = std::max(zend, begin_end[i].second);
    }

1058 1059
    //~ delete[] zvals;
    //~ delete[] onAxisField;
gsell's avatar
gsell committed
1060 1061 1062 1063 1064

    double z = zbegin;
    double dz = 0.5 * p / sqrt(1.0 + p * p) * cdt;
    while(z + dz < zend && z + dz > zbegin) {
        z += dz;
1065 1066
        for(size_t i = 0; i < numFieldmaps(); ++ i) {
            if(z >= begin_end[i].first && z <= begin_end[i].second) {
gsell's avatar
gsell committed
1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077
                double phase = multiFrequencies_m[i] * t + multiPhases_m[i];
                double scale = multiScales_m[i];
                double ez = scale * gsl_spline_eval(onAxisInterpolants[i], z, onAxisAccel[i]);
                p += cos(phase) * ez * q * cdt / mass;
            }
        }
        dz = 0.5 * p / sqrt(1.0 + p * p) * cdt;
        z += dz;
        t += dt;
    }

1078
    for(size_t i = 0; i < numFieldmaps(); ++ i) {
gsell's avatar
gsell committed
1079 1080 1081
        gsl_spline_free(onAxisInterpolants[i]);
        gsl_interp_accel_free(onAxisAccel[i]);
    }
1082 1083
    //~ delete[] onAxisInterpolants;
    //~ delete[] onAxisAccel;
gsell's avatar
gsell committed
1084 1085 1086 1087 1088

    const double beta = sqrt(1. - 1 / (p * p + 1.));
    const double tErr  = (z - zend) / (Physics::c * beta);

    return pair<double, double>(p, t - tErr);
1089
}