AutophaseTracker.cpp 20.2 KB
Newer Older
kraus's avatar
kraus committed
1
#include "Algorithms/AutophaseTracker.h"
kraus's avatar
kraus committed
2
#include "Utilities/ClassicField.h"
kraus's avatar
kraus committed
3 4
#include "AbstractObjects/OpalData.h"
#include "Utilities/OpalOptions.h"
kraus's avatar
kraus committed
5
#include "Utilities/Timer.h"
kraus's avatar
kraus committed
6

7 8 9
#include <iostream>
#include <fstream>

kraus's avatar
kraus committed
10 11
#define DTFRACTION 2

kraus's avatar
kraus committed
12 13
extern Inform *gmsg;

kraus's avatar
kraus committed
14 15 16 17 18
AutophaseTracker::AutophaseTracker(const Beamline &beamline,
                                   const PartData& data,
                                   const double &T0,
                                   double initialR,
                                   double initialP):
kraus's avatar
kraus committed
19 20
    DefaultVisitor(beamline, false, false),
    itsBunch_m(&data),
kraus's avatar
kraus committed
21 22
    itsPusher_m(data),
    particleLayout_m(NULL)
kraus's avatar
kraus committed
23
{
kraus's avatar
kraus committed
24
    particleLayout_m = new Layout_t();
kraus's avatar
kraus committed
25 26 27
    timeIntegrationTimer_m = IpplTimings::getTimer("AP: time integration");
    fieldEvaluationTimer_m = IpplTimings::getTimer("AP: field evaluation");

kraus's avatar
kraus committed
28
    itsBunch_m.initialize(particleLayout_m);
kraus's avatar
kraus committed
29
    itsBunch_m.setT(T0);
kraus's avatar
kraus committed
30

kraus's avatar
kraus committed
31 32 33 34 35 36
    if(Ippl::myNode() == 0) {
        itsBunch_m.create(1);
        itsBunch_m.R[0] = Vector_t(0.0, 0.0, initialR);
        itsBunch_m.P[0] = Vector_t(0.0, 0.0, initialP);
        itsBunch_m.Bin[0] = 0;
        itsBunch_m.Q[0] = itsBunch_m.getChargePerParticle();
kraus's avatar
kraus committed
37
        itsBunch_m.PType[0] = ParticleType::REGULAR;
kraus's avatar
kraus committed
38 39 40
        itsBunch_m.LastSection[0] = 0;
    }

kraus's avatar
kraus committed
41
    visitBeamline(beamline);
kraus's avatar
kraus committed
42
    itsOpalBeamline_m.prepareSections();
kraus's avatar
kraus committed
43 44 45 46 47 48 49 50 51
}

AutophaseTracker::~AutophaseTracker()
{ }

void AutophaseTracker::execute(const std::queue<double> &dtAllTracks,
                               const std::queue<double> &maxZ,
                               const std::queue<unsigned long long> &maxTrackSteps) {

52 53 54 55
    if (Ippl::myNode() != 0) {
        receiveCavityPhases();
        return;
    }
kraus's avatar
kraus committed
56 57
    OPALTimer::Timer clock;
    *gmsg << level1 << "\n\nstart autophasing at " << clock.time() << "\n" << endl;
kraus's avatar
kraus committed
58

kraus's avatar
kraus committed
59 60 61 62 63 64
    Vector_t rmin, rmax;
    Vector_t &refR = itsBunch_m.R[0];
    Vector_t &refP = itsBunch_m.P[0];
    std::queue<double> dtAutoPhasing(dtAllTracks);
    std::queue<double> maxSPosAutoPhasing(maxZ);
    std::queue<unsigned long long> maxStepsAutoPhasing(maxTrackSteps);
kraus's avatar
kraus committed
65

kraus's avatar
kraus committed
66
    maxSPosAutoPhasing.back() = itsOpalBeamline_m.calcBeamlineLenght();
kraus's avatar
kraus committed
67

kraus's avatar
kraus committed
68
    size_t step = 0;
kraus's avatar
kraus committed
69 70 71
    const int dtfraction = DTFRACTION;
    double newDt = dtAutoPhasing.front() / dtfraction;
    double scaleFactor = newDt * Physics::c;
kraus's avatar
kraus committed
72

kraus's avatar
kraus committed
73 74
    itsBunch_m.LastSection[0] -= 1;
    itsOpalBeamline_m.getSectionIndexAt(refR, itsBunch_m.LastSection[0]);
kraus's avatar
kraus committed
75
    itsOpalBeamline_m.print(*gmsg);
kraus's avatar
kraus committed
76

kraus's avatar
kraus committed
77
    itsOpalBeamline_m.switchAllElements();
kraus's avatar
kraus committed
78

kraus's avatar
kraus committed
79
    std::shared_ptr<Component> cavity = NULL;
kraus's avatar
kraus committed
80
    while (true) {
kraus's avatar
kraus committed
81
        std::shared_ptr<Component> next = getNextCavity(cavity);
kraus's avatar
kraus committed
82 83
        if (next == NULL) break;

kraus's avatar
kraus committed
84 85
        double endNext = getEndCavity(next);
        if (refR(2) < endNext) break;
kraus's avatar
kraus committed
86 87 88

        cavity = next;
    }
kraus's avatar
kraus committed
89

kraus's avatar
kraus committed
90 91
    while (maxStepsAutoPhasing.size() > 0) {
        maxStepsAutoPhasing.front() = maxStepsAutoPhasing.front() * dtfraction + step;
kraus's avatar
kraus committed
92 93
        newDt = dtAutoPhasing.front() / dtfraction;
        scaleFactor = newDt * Physics::c;
kraus's avatar
kraus committed
94 95 96
        Vector_t vscaleFactor = Vector_t(scaleFactor);
        itsBunch_m.setdT(newDt);
        itsBunch_m.dt = newDt;
kraus's avatar
kraus committed
97

kraus's avatar
kraus committed
98 99
        INFOMSG(level2 << "\n"
                << "*************************************************\n"
kraus's avatar
kraus committed
100
                << "new set of dt, zstop and max number of time steps\n"
kraus's avatar
kraus committed
101
                << "*************************************************\n\n"
kraus's avatar
kraus committed
102
                << "start at t = " << itsBunch_m.getT() << " [s], "
kraus's avatar
kraus committed
103
                << "zstop at z = " << maxSPosAutoPhasing.front() << " [m],\n"
kraus's avatar
kraus committed
104 105 106
                << "dt = " << itsBunch_m.dt[0] << " [s], "
                << "step = " << step << ", "
                << "R =  " << refR << " [m]" << endl);
kraus's avatar
kraus committed
107

kraus's avatar
kraus committed
108
        for(; step < maxStepsAutoPhasing.front(); ++ step) {
kraus's avatar
kraus committed
109
            std::shared_ptr<Component> next = getNextCavity(cavity);
kraus's avatar
kraus committed
110 111 112 113 114 115 116 117 118 119 120 121 122 123
            if (next == NULL) break;

            double beginNext, endNext;
            next->getDimensions(beginNext, endNext);
            double sStop = beginNext;
            bool reachesNextCavity = true;
            if (sStop > maxSPosAutoPhasing.front()) {
                sStop = maxSPosAutoPhasing.front();
                reachesNextCavity = false;
            }

            advanceTime(step, maxStepsAutoPhasing.front(), sStop);
            if (step >= maxStepsAutoPhasing.front() || !reachesNextCavity) break;

kraus's avatar
kraus committed
124 125
            INFOMSG(level2 << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
            INFOMSG(level2 << "\nFound " << next->getName()
kraus's avatar
kraus committed
126 127
                    << " at " << sStop << " [m], "
                    << "step  " << step << ", "
kraus's avatar
kraus committed
128 129
                    << "t = " << itsBunch_m.getT() << " [s], "
                    << "E = " << getEnergyMeV(refP) << " [MeV]\n\n"
130
                    << "start phase scan ... " << endl;);
kraus's avatar
kraus committed
131 132 133 134 135 136

            double guess = guessCavityPhase(next);
            optimizeCavityPhase(next,
                                guess,
                                step,
                                newDt);
kraus's avatar
kraus committed
137

138
            INFOMSG(level2 << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl);
kraus's avatar
kraus committed
139 140 141 142 143 144 145 146 147

            track(endNext,
                  step,
                  dtAutoPhasing,
                  maxSPosAutoPhasing,
                  maxStepsAutoPhasing);

            cavity = next;
        }
kraus's avatar
kraus committed
148

kraus's avatar
kraus committed
149
        if (dtAutoPhasing.size() == 0) break;
kraus's avatar
kraus committed
150

kraus's avatar
kraus committed
151 152 153 154
        dtAutoPhasing.pop();
        maxSPosAutoPhasing.pop();
        maxStepsAutoPhasing.pop();
    }
155

156 157 158 159
    INFOMSG(level2 << "step = " << step << ", spos = " << refR(2) << " [m], "
            << "t= " << itsBunch_m.getT() << " [s], " << "E = " << getEnergyMeV(refP) << " [MeV] "
            << endl);

160
    printCavityPhases();
161 162
    sendCavityPhases();

163
    *gmsg << level1 << "\n\nfinished autophasing " << clock.time() << "\n\n"  << endl;
kraus's avatar
kraus committed
164 165 166 167 168 169 170 171 172
}

void AutophaseTracker::track(double uptoSPos,
                             size_t &step,
                             std::queue<double> &dtAutoPhasing,
                             std::queue<double> &maxSPosAutoPhasing,
                             std::queue<unsigned long long> &maxStepsAutoPhasing) {
    Vector_t &refR = itsBunch_m.R[0];
    Vector_t &refP = itsBunch_m.P[0];
173

kraus's avatar
kraus committed
174 175
    double t = itsBunch_m.getT();
    const int dtfraction = DTFRACTION;
kraus's avatar
kraus committed
176

kraus's avatar
kraus committed
177 178 179 180 181 182 183
    while (maxStepsAutoPhasing.size() > 0) {
        maxStepsAutoPhasing.front() = maxStepsAutoPhasing.front() * dtfraction + step;
        double newDt = dtAutoPhasing.front() / dtfraction;
        double scaleFactor = newDt * Physics::c;
        Vector_t vscaleFactor = Vector_t(scaleFactor);
        itsBunch_m.setdT(newDt);
        itsBunch_m.dt[0] = newDt;
kraus's avatar
kraus committed
184

kraus's avatar
kraus committed
185
        for(; step < maxStepsAutoPhasing.front(); ++ step) {
kraus's avatar
kraus committed
186
            IpplTimings::startTimer(timeIntegrationTimer_m);
kraus's avatar
kraus committed
187

kraus's avatar
kraus committed
188
            refR /= vscaleFactor;
kraus's avatar
kraus committed
189
            itsPusher_m.push(refR, refP, itsBunch_m.dt[0]);
kraus's avatar
kraus committed
190
            refR *= vscaleFactor;
kraus's avatar
kraus committed
191

kraus's avatar
kraus committed
192
            IpplTimings::stopTimer(timeIntegrationTimer_m);
kraus's avatar
kraus committed
193

kraus's avatar
kraus committed
194
            evaluateField();
kraus's avatar
kraus committed
195

kraus's avatar
kraus committed
196
            IpplTimings::startTimer(timeIntegrationTimer_m);
kraus's avatar
kraus committed
197

kraus's avatar
kraus committed
198
            itsPusher_m.kick(refR, refP, itsBunch_m.Ef[0], itsBunch_m.Bf[0], itsBunch_m.dt[0]);
kraus's avatar
kraus committed
199

kraus's avatar
kraus committed
200
            refR /= vscaleFactor;
kraus's avatar
kraus committed
201
            itsPusher_m.push(refR, refP, itsBunch_m.dt[0]);
kraus's avatar
kraus committed
202
            refR *= vscaleFactor;
kraus's avatar
kraus committed
203

kraus's avatar
kraus committed
204 205 206
            t += itsBunch_m.getdT();
            itsBunch_m.setT(t);
            IpplTimings::stopTimer(timeIntegrationTimer_m);
kraus's avatar
kraus committed
207

kraus's avatar
kraus committed
208
            if (refR(2) > uptoSPos) return;
kraus's avatar
kraus committed
209

kraus's avatar
kraus committed
210 211
            if(refR(2) > maxSPosAutoPhasing.front())
                maxStepsAutoPhasing.front() = step;
kraus's avatar
kraus committed
212

kraus's avatar
kraus committed
213
            if(step % 5000 == 0) {
214 215 216
                INFOMSG(level2 << "step = " << step << ", spos = " << refR(2) << " [m], "
                        << "t= " << itsBunch_m.getT() << " [s], " << "E = " << getEnergyMeV(refP) << " [MeV] "
                        << endl);
kraus's avatar
kraus committed
217 218
            }
        }
kraus's avatar
kraus committed
219

kraus's avatar
kraus committed
220 221 222 223
        dtAutoPhasing.pop();
        maxSPosAutoPhasing.pop();
        maxStepsAutoPhasing.pop();
    }
kraus's avatar
kraus committed
224 225
}

kraus's avatar
kraus committed
226 227
std::shared_ptr<Component> AutophaseTracker::getNextCavity(const std::shared_ptr<Component> &current) {
    auto cavities = itsOpalBeamline_m.getSuccessors(current);
kraus's avatar
kraus committed
228 229

    if (cavities.size() > 0) {
kraus's avatar
kraus committed
230
        return cavities.front();
kraus's avatar
kraus committed
231 232
    }

kraus's avatar
kraus committed
233 234
    std::shared_ptr<Component> noCavity;
    return noCavity;
kraus's avatar
kraus committed
235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250
}

void AutophaseTracker::advanceTime(size_t &step, const size_t maxSteps, const double beginNextCavity) {
    Vector_t &refR = itsBunch_m.R[0];
    Vector_t &refP = itsBunch_m.P[0];

    const double distance = beginNextCavity - refR(2);
    if (distance < 0.0) return;

    const double beta = refP(2) / sqrt(dot(refP, refP) + 1.0);
    size_t  numSteps = std::floor(distance / (beta * Physics::c * itsBunch_m.dt[0]));
    numSteps = std::min(numSteps, maxSteps - step);
    double timeDifference = numSteps * itsBunch_m.dt[0];

    itsBunch_m.setT(timeDifference + itsBunch_m.getT());
    step += numSteps;
kraus's avatar
kraus committed
251
    refR = Vector_t(0.0, 0.0, refR(2) + beta * Physics::c * timeDifference);
kraus's avatar
kraus committed
252 253 254
    refP = Vector_t(0.0, 0.0, sqrt(dot(refP, refP)));
}

kraus's avatar
kraus committed
255
double AutophaseTracker::guessCavityPhase(const std::shared_ptr<Component> &cavity) {
kraus's avatar
kraus committed
256 257 258 259 260 261
    Vector_t &refR = itsBunch_m.R[0];
    Vector_t &refP = itsBunch_m.P[0];

    double cavityStart = getBeginCavity(cavity);
    double orig_phi = 0.0, Phimax = 0.0;

kraus's avatar
kraus committed
262
    const double beta = std::max(refP(2) / sqrt(dot(refP, refP) + 1.0), 1e-6);
kraus's avatar
kraus committed
263 264 265
    const double tErr  = (cavityStart - refR(2)) / (Physics::c * beta);

    bool apVeto;
266
    if(cavity->getType() == ElementBase::TRAVELINGWAVE) {
kraus's avatar
kraus committed
267
        TravelingWave *element = static_cast<TravelingWave *>(cavity.get());
kraus's avatar
kraus committed
268 269 270
        orig_phi = element->getPhasem();
        apVeto = element->getAutophaseVeto();
        if(apVeto) {
271
            INFOMSG(level2 << " ----> APVETO -----> "
kraus's avatar
kraus committed
272 273 274
                    << element->getName() <<  endl);
            return orig_phi;
        }
275 276
        INFOMSG(level2 << cavity->getName() << ", "
                << "phase = " << orig_phi << " rad" << endl);
kraus's avatar
kraus committed
277 278 279 280 281 282

        Phimax = element->getAutoPhaseEstimate(getEnergyMeV(refP),
                                               itsBunch_m.getT() + tErr,
                                               itsBunch_m.getQ(),
                                               itsBunch_m.getM() * 1e-6);
    } else {
kraus's avatar
kraus committed
283
        RFCavity *element = static_cast<RFCavity *>(cavity.get());
kraus's avatar
kraus committed
284 285 286
        orig_phi = element->getPhasem();
        apVeto = element->getAutophaseVeto();
        if(apVeto) {
287
            INFOMSG(level2 << " ----> APVETO -----> "
kraus's avatar
kraus committed
288 289 290
                    << element->getName() << endl);
            return orig_phi;
        }
291 292
        INFOMSG(level2 << cavity->getName() << ", "
                << "phase = " << orig_phi << " rad" << endl);
kraus's avatar
kraus committed
293 294 295 296 297 298 299 300 301 302

        Phimax = element->getAutoPhaseEstimate(getEnergyMeV(refP),
                                               itsBunch_m.getT() + tErr,
                                               itsBunch_m.getQ(),
                                               itsBunch_m.getM() * 1e-6);
    }

    return Phimax;
}

kraus's avatar
kraus committed
303
double AutophaseTracker::optimizeCavityPhase(const std::shared_ptr<Component> &cavity,
kraus's avatar
kraus committed
304 305 306 307
                                             double initialPhase,
                                             size_t currentStep,
                                             double dt) {

308
    if(cavity->getType() == ElementBase::TRAVELINGWAVE) {
kraus's avatar
kraus committed
309
        if (static_cast<TravelingWave *>(cavity.get())->getAutophaseVeto()) return initialPhase;
kraus's avatar
kraus committed
310
    } else {
kraus's avatar
kraus committed
311
        if (static_cast<RFCavity *>(cavity.get())->getAutophaseVeto()) return initialPhase;
kraus's avatar
kraus committed
312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363
    }

    itsBunch_m.setdT(dt);
    itsBunch_m.dt[0] = dt;
    double Phimax = initialPhase;
    double phi = initialPhase;
    double originalPhase = 0.0, PhiAstra = 0.0;
    double dphi = Physics::pi / 360.0;
    double cavityStart = getBeginCavity(cavity);
    const int numRefinements = Options::autoPhase;

    int j = -1;

    double E = APtrack(cavity, cavityStart, phi);
    double Emax = E;
    do {
        j ++;
        Emax = E;
        initialPhase = phi;
        phi -= dphi;
        E = APtrack(cavity, cavityStart, phi);
    } while(E > Emax);

    if(j == 0) {
        phi = initialPhase;
        E = Emax;
        j = -1;
        do {
            j ++;
            Emax = E;
            initialPhase = phi;
            phi += dphi;
            E = APtrack(cavity, cavityStart, phi);
        } while(E > Emax);
    }

    for(int rl = 0; rl < numRefinements; ++ rl) {
        dphi /= 2.;
        phi = initialPhase - dphi;
        E = APtrack(cavity, cavityStart, phi);
        if(E > Emax) {
            initialPhase = phi;
            Emax = E;
        } else {
            phi = initialPhase + dphi;
            E = APtrack(cavity, cavityStart, phi);
            if(E > Emax) {
                initialPhase = phi;
                Emax = E;
            }
        }
    }
364
    Phimax = std::fmod(initialPhase, 2 * Physics::pi);
kraus's avatar
kraus committed
365

366
    if(cavity->getType() == ElementBase::TRAVELINGWAVE) {
kraus's avatar
kraus committed
367
        TravelingWave *element = static_cast<TravelingWave *>(cavity.get());
kraus's avatar
kraus committed
368
        originalPhase = element->getPhasem();
369 370
        double newPhase = std::fmod(originalPhase + Phimax, 2 * Physics::pi);
        element->updatePhasem(newPhase);
kraus's avatar
kraus committed
371
    } else {
kraus's avatar
kraus committed
372
        RFCavity *element = static_cast<RFCavity *>(cavity.get());
kraus's avatar
kraus committed
373
        originalPhase = element->getPhasem();
374 375
        double newPhase = std::fmod(originalPhase + Phimax, 2 * Physics::pi);
        element->updatePhasem(newPhase);
kraus's avatar
kraus committed
376 377
    }

378
    PhiAstra = std::fmod((Phimax * Physics::rad2deg) + 90.0, 360.0);
kraus's avatar
kraus committed
379 380 381

    INFOMSG(cavity->getName() << "_phi = "  << Phimax << " rad / "
            << Phimax *Physics::rad2deg <<  " deg, AstraPhi = " << PhiAstra << " deg,\n"
kraus's avatar
kraus committed
382
            << "E = " << Emax << " (MeV), " << "phi_nom = " << originalPhase *Physics::rad2deg << "\n");
kraus's avatar
kraus committed
383 384 385

    OpalData::getInstance()->setMaxPhase(cavity->getName(), Phimax);

386 387 388 389
    double cavBegin, cavEnd;
    cavity->getDimensions(cavBegin, cavEnd);
    OpalData::getInstance()->addEnergyData(cavEnd, Emax * 1e6);

kraus's avatar
kraus committed
390
    return Phimax + originalPhase;
kraus's avatar
kraus committed
391 392
}

kraus's avatar
kraus committed
393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413
void AutophaseTracker::evaluateField() {
    IpplTimings::startTimer(fieldEvaluationTimer_m);

    Vector_t &refR = itsBunch_m.R[0];

    itsBunch_m.Ef = Vector_t(0.0);
    itsBunch_m.Bf = Vector_t(0.0);

    long ls = itsBunch_m.LastSection[0];
    itsOpalBeamline_m.getSectionIndexAt(refR, ls);
    itsBunch_m.LastSection[0] = ls;

    Vector_t externalE, externalB;
    itsOpalBeamline_m.getFieldAt(0, refR, ls, itsBunch_m.getT() + itsBunch_m.dt[0] / 2., externalE, externalB);

    itsBunch_m.Ef[0] += externalE;
    itsBunch_m.Bf[0] += externalB;

    IpplTimings::stopTimer(fieldEvaluationTimer_m);
}

kraus's avatar
kraus committed
414
double AutophaseTracker::APtrack(const std::shared_ptr<Component> &cavity, double cavityStart, const double &phi) const {
kraus's avatar
kraus committed
415 416
    const Vector_t &refR = itsBunch_m.R[0];
    const Vector_t &refP = itsBunch_m.P[0];
kraus's avatar
kraus committed
417

kraus's avatar
kraus committed
418 419 420
    double beta = refP(2) / std::sqrt(dot(refP, refP) + 1.0);
    beta = std::max(beta, 0.0001);
    double tErr  = (cavityStart - refR(2)) / (Physics::c * beta);
kraus's avatar
kraus committed
421

kraus's avatar
kraus committed
422
    double initialPhase = 0.0;
kraus's avatar
kraus committed
423
    double finalMomentum = 0.0;
424
    if(cavity->getType() == ElementBase::TRAVELINGWAVE) {
kraus's avatar
kraus committed
425 426
        TravelingWave *tws = static_cast<TravelingWave *>(cavity.get());
        initialPhase = tws->getPhasem();
kraus's avatar
kraus committed
427
        tws->updatePhasem(phi);
kraus's avatar
kraus committed
428
        std::pair<double, double> pe = tws->trackOnAxisParticle(refP(2),
kraus's avatar
kraus committed
429 430 431 432 433
                                                                itsBunch_m.getT() + tErr,
                                                                itsBunch_m.dt[0],
                                                                itsBunch_m.getQ(),
                                                                itsBunch_m.getM() * 1e-6);
        finalMomentum = pe.first;
kraus's avatar
kraus committed
434
        tws->updatePhasem(initialPhase);
kraus's avatar
kraus committed
435
    } else {
kraus's avatar
kraus committed
436 437
        RFCavity *rfc = static_cast<RFCavity *>(cavity.get());
        initialPhase = rfc->getPhasem();
kraus's avatar
kraus committed
438 439
        rfc->updatePhasem(phi);

kraus's avatar
kraus committed
440
        std::pair<double, double> pe = rfc->trackOnAxisParticle(refP(2),
kraus's avatar
kraus committed
441 442 443 444 445
                                                                itsBunch_m.getT() + tErr,
                                                                itsBunch_m.dt[0],
                                                                itsBunch_m.getQ(),
                                                                itsBunch_m.getM() * 1e-6);
        finalMomentum = pe.first;
kraus's avatar
kraus committed
446
        rfc->updatePhasem(initialPhase);
kraus's avatar
kraus committed
447 448 449 450 451 452 453
    }
    double finalGamma = sqrt(1.0 + finalMomentum * finalMomentum);
    double finalKineticEnergy = (finalGamma - 1.0) * itsBunch_m.getM() * 1e-6;
    return finalKineticEnergy;
}


454 455 456 457 458
void AutophaseTracker::sendCavityPhases() {

    int tag = 101;
    Message *mess = new Message();
    putMessage(*mess, OpalData::getInstance()->getNumberOfMaxPhases());
kraus's avatar
kraus committed
459

460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476
    {
        typedef std::vector<MaxPhasesT>::iterator iterator_t;
        iterator_t it = OpalData::getInstance()->getFirstMaxPhases();
        iterator_t end = OpalData::getInstance()->getLastMaxPhases();
        for(; it != end; ++ it) {
            putMessage(*mess, (*it).first);
            putMessage(*mess, (*it).second);
        }
    }
    {
        typedef energyEvolution_t::iterator iterator_t;
        iterator_t it = OpalData::getInstance()->getFirstEnergyData();
        iterator_t end = OpalData::getInstance()->getLastEnergyData();
        for(; it != end; ++ it) {
            putMessage(*mess, (*it).first);
            putMessage(*mess, (*it).second);
        }
477 478
    }
    Ippl::Comm->broadcast_all(mess, tag);
kraus's avatar
kraus committed
479

480
    delete mess;
kraus's avatar
kraus committed
481 482
}

483 484 485 486 487 488 489 490
void AutophaseTracker::receiveCavityPhases() {

    int parent = 0;
    int tag = 101;
    int nData = 0;
    Message *mess = Ippl::Comm->receive_block(parent, tag);
    getMessage(*mess, nData);

491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510
    {
        typedef std::vector<MaxPhasesT>::iterator iterator_t;
        for(int i = 0; i < nData; i++) {
            std::string elName;
            double maxPhi;
            getMessage(*mess, elName);
            getMessage(*mess, maxPhi);
            OpalData::getInstance()->setMaxPhase(elName, maxPhi);
        }
    }
    {
        typedef energyEvolution_t::iterator iterator_t;
        for(int i = 0; i < nData; i++) {
            double spos;
            double maxEnergy;
            getMessage(*mess, spos);
            getMessage(*mess, maxEnergy);
            OpalData::getInstance()->addEnergyData(spos, maxEnergy);
        }
    }
511
    delete mess;
512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 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
}

void AutophaseTracker::printCavityPhases() {
    std::shared_ptr<Component> cavity = NULL;
    const double RADDEG = 180.0 / Physics::pi;
    const double globalTimeShift = OpalData::getInstance()->getGlobalPhaseShift();

    *gmsg << level1 << "\n-------------------------------------------------------------------------------------\n";
    while (true) {
        std::shared_ptr<Component> next = getNextCavity(cavity);
        if (next == NULL) break;

        std::string name = next->getName();
        double frequency;
        double phase;

        if (next->getType() == ElementBase::TRAVELINGWAVE) {
            phase = static_cast<TravelingWave *>(next.get())->getPhasem();
	    frequency = static_cast<TravelingWave *>(next.get())->getFrequencym();
        } else {
            phase = static_cast<RFCavity *>(next.get())->getPhasem();
	    frequency = static_cast<RFCavity *>(next.get())->getFrequencym();
        }

        phase -= globalTimeShift * frequency;

        *gmsg << (cavity == NULL? "": "\n")
            << name
            << ": phi = phi_nom + phi_maxE + global phase shift = " << phase * RADDEG << " degree, "
            << "(global phase shift = " << -globalTimeShift *frequency *RADDEG << " degree) \n";

        cavity = next;
    }
    *gmsg << "-------------------------------------------------------------------------------------\n"
          << endl;
}

void AutophaseTracker::save(const std::string &fname) {
    std::ofstream out(fname);

    std::shared_ptr<Component> cavity = NULL;
    const double globalTimeShift = OpalData::getInstance()->getGlobalPhaseShift();

    while (true) {
        std::shared_ptr<Component> next = getNextCavity(cavity);
        if (next == NULL) break;

        std::string name = next->getName();
        double frequency;
        double phase;

        if (next->getType() == ElementBase::TRAVELINGWAVE) {
            phase = static_cast<TravelingWave *>(next.get())->getPhasem();
	    frequency = static_cast<TravelingWave *>(next.get())->getFrequencym();
        } else {
            phase = static_cast<RFCavity *>(next.get())->getPhasem();
	    frequency = static_cast<RFCavity *>(next.get())->getFrequencym();
        }

        phase -= globalTimeShift * frequency;

        out << name << "_lag = " << phase << ";\n";

        cavity = next;
    }

    out.close();
579
}