CavityAutophaser.cpp

Go to the documentation of this file.

//
// Class CavityAutophaser
//
// This class determines the phase of an RF cavity for which the reference particle
// is accelerated to the highest energy.
//
// Copyright (c) 2016,       Christof Metzger-Kraus, Helmholtz-Zentrum Berlin, Germany
//               2017 - 2020 Christof Metzger-Kraus
//
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "Ippl.h"
#include "OPALTypes.h"
 
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/TravelingWave.h"
#include "AbstractObjects/OpalData.h"
#include "Algorithms/CavityAutophaser.h"
//
#include "Physics/Units.h"
#include "Utilities/OpalException.h"
#include "Utilities/Options.h"
#include "Utilities/Util.h"
 
#include <fstream>
#include <iostream>
#include <string>
 
extern Inform* gmsg;
 
CavityAutophaser::CavityAutophaser(const PartData& ref, std::shared_ptr<Component> cavity)
    : itsReference_m(ref), itsCavity_m(cavity) {
    double zbegin = 0.0, zend = 0.0;
    cavity->getDimensions(zbegin, zend);
    initialR_m = Vector_t<double, 3>(0, 0, zbegin);
}
 
CavityAutophaser::~CavityAutophaser() {
}
 
double CavityAutophaser::getPhaseAtMaxEnergy(
    const Vector_t<double, 3>& R, const Vector_t<double, 3>& P, double t, double dt) {
    if (!(itsCavity_m->getType() == ElementType::TRAVELINGWAVE
          || itsCavity_m->getType() == ElementType::RFCAVITY)) {
        throw OpalException(
            "CavityAutophaser::getPhaseAtMaxEnergy()", "given element is not a cavity");
    }
 
    initialP_m = P;  // \todo need to check ... Vector_t<double, 3>(0, 0, std::sqrt(dot(P, P)));
 
    RFCavity* element    = static_cast<RFCavity*>(itsCavity_m.get());
    bool apVeto          = element->getAutophaseVeto();
    bool isDCGun         = false;
    double originalPhase = element->getPhasem();
    double tErr = (initialR_m(2) - R(2)) * std::sqrt(dot(P, P) + 1.0) / (P(2) * Physics::c);
    double optimizedPhase = 0.0;
    double finalEnergy    = 0.0;
    double newPhase       = 0.0;
    double amplitude      = element->getAmplitudem();
    double basePhase      = std::fmod(element->getFrequencym() * (t + tErr), Physics::two_pi);
    double frequency      = element->getFrequencym();
 
    if ((!apVeto) && frequency <= (1.0 + 1e-6) * Physics::two_pi) {  // DC gun
        optimizedPhase = (amplitude * itsReference_m.getQ() > 0.0 ? 0.0 : Physics::pi);
        element->setPhasem(optimizedPhase + originalPhase);
        element->setAutophaseVeto();
 
        originalPhase += optimizedPhase;
        OpalData::getInstance()->setMaxPhase(itsCavity_m->getName(), originalPhase);
 
        apVeto  = true;
        isDCGun = true;
    }
 
    std::stringstream ss;
    for (char c : itsCavity_m->getName()) {
        ss << std::setw(2) << std::left << c;
    }
    *ippl::Info << level1 << "\n* ************* " << std::left << std::setw(68) << std::setfill('*')
                << ss.str() << std::setfill(' ') << endl;
    if (!apVeto) {
        double initialEnergy = Util::getKineticEnergy(P, itsReference_m.getM()) * Units::eV2MeV;
        double AstraPhase    = 0.0;
        double designEnergy  = element->getDesignEnergy();
 
        if (amplitude < 0.0) {
            amplitude = -amplitude;
            element->setAmplitudem(amplitude);
        }
 
        double initialPhase = guessCavityPhase(t + tErr);
 
        if (amplitude == 0.0 && designEnergy <= 0.0) {
            throw OpalException(
                "CavityAutophaser::getPhaseAtMaxEnergy()",
                "neither amplitude or design energy given to cavity " + element->getName());
        }
 
        if (designEnergy > 0.0) {
            const double length = itsCavity_m->getElementLength();
            if (length <= 0.0) {
                throw OpalException(
                    "CavityAutophaser::getPhaseAtMaxEnergy()",
                    "length of cavity " + element->getName() + " is zero");
            }
 
            amplitude =
                2 * (designEnergy - initialEnergy) / (std::abs(itsReference_m.getQ()) * length);
 
            element->setAmplitudem(amplitude);
 
            int count = 0;
            while (count < 1000) {
                initialPhase = guessCavityPhase(t + tErr);
                auto status  = optimizeCavityPhase(initialPhase, t + tErr, dt);
 
                optimizedPhase = status.first;
                finalEnergy    = status.second;
 
                if (std::abs(designEnergy - finalEnergy) < 1e-7)
                    break;
 
                amplitude *= std::abs(designEnergy / finalEnergy);
                element->setAmplitudem(amplitude);
                initialPhase = optimizedPhase;
 
                ++count;
            }
        }
 
        auto status = optimizeCavityPhase(initialPhase, t + tErr, dt);
 
        optimizedPhase = status.first;
        finalEnergy    = status.second;
 
        AstraPhase = std::fmod(optimizedPhase + Physics::pi / 2 + Physics::two_pi, Physics::two_pi);
        newPhase   = std::fmod(originalPhase + optimizedPhase + Physics::two_pi, Physics::two_pi);
        element->setPhasem(newPhase);
        element->setAutophaseVeto();
 
        auto opal = OpalData::getInstance();
 
        opal->setMaxPhase(itsCavity_m->getName(), newPhase);
 
        newPhase = std::fmod(newPhase + basePhase, Physics::two_pi);
 
        if (!opal->isOptimizerRun()) {
            std::string fname = Util::combineFilePath(
                {opal->getAuxiliaryOutputDirectory(), itsCavity_m->getName() + "_AP.dat"});
            std::ofstream out(fname);
            track(t + tErr, dt, newPhase, &out);
            out.close();
        } else {
            track(t + tErr, dt, newPhase, nullptr);
        }
 
        *ippl::Info << level1 << std::fixed << std::setprecision(4) << itsCavity_m->getName()
                    << "_phi1 = " << newPhase * Units::rad2deg << " [deg], "
                    << "corresp. in Astra = " << AstraPhase * Units::rad2deg << " [deg],\n"
                    << "E = " << finalEnergy << " [MeV], "
                    << "phi_nom = " << originalPhase * Units::rad2deg << " [deg]\n"
                    << "Ez_0 = " << amplitude << " [MV/m]"
                    << "\n"
                    << "time = " << (t + tErr) * Units::s2ns << " [ns], dt = " << dt * Units::s2ps
                    << " [ps]" << endl;
 
    } else {
        auto status = optimizeCavityPhase(originalPhase, t + tErr, dt);
 
        finalEnergy = status.second;
 
        originalPhase = std::fmod(originalPhase, Physics::two_pi);
        double AstraPhase =
            std::fmod(optimizedPhase + Physics::pi / 2 + Physics::two_pi, Physics::two_pi);
 
        if (!isDCGun) {
            *ippl::Info << level1 << ">>>>>> APVETO >>>>>> " << endl;
        }
        *ippl::Info << level1 << std::fixed << std::setprecision(4) << itsCavity_m->getName()
                    << "_phi2 = " << originalPhase * Units::rad2deg << " [deg], "
                    << "corresp. in Astra = " << AstraPhase * Units::rad2deg << " [deg],\n"
                    << "E = " << finalEnergy << " [MeV], "
                    << "phi_nom = " << originalPhase * Units::rad2deg << " [deg]\n"
                    << "Ez_0 = " << amplitude << " [MV/m]"
                    << "\n"
                    << "time = " << (t + tErr) * Units::s2ns << " [ns], dt = " << dt * Units::s2ps
                    << " [ps]" << endl;
        if (!isDCGun) {
            *ippl::Info << level1 << " <<<<<< APVETO <<<<<< " << endl;
        }
 
        optimizedPhase = originalPhase;
    }
    *ippl::Info << level1 << "* " << std::right << std::setw(83) << std::setfill('*') << "*\n"
                << std::setfill(' ') << endl;
 
    return optimizedPhase;
}
 
double CavityAutophaser::guessCavityPhase(double t) {
    const Vector_t<double, 3>& refP = initialP_m;
    double Phimax                   = 0.0;
    bool apVeto;
    RFCavity* element = static_cast<RFCavity*>(itsCavity_m.get());
    double orig_phi   = element->getPhasem();
    apVeto            = element->getAutophaseVeto();
    if (apVeto) {
        return orig_phi;
    }
 
    Phimax = element->getAutoPhaseEstimate(
        Util::getKineticEnergy(refP, itsReference_m.getM()) * Units::eV2MeV, t,
        itsReference_m.getQ(), itsReference_m.getM() * Units::eV2MeV);
 
    return std::fmod(Phimax + Physics::two_pi, Physics::two_pi);
}
 
std::pair<double, double> CavityAutophaser::optimizeCavityPhase(
    double initialPhase, double t, double dt) {
    RFCavity* element    = static_cast<RFCavity*>(itsCavity_m.get());
    double originalPhase = element->getPhasem();
 
    if (element->getAutophaseVeto()) {
        double basePhase = std::fmod(element->getFrequencym() * t, Physics::two_pi);
        double phase     = std::fmod(originalPhase - basePhase + Physics::two_pi, Physics::two_pi);
        double E         = track(t, dt, phase);
        std::pair<double, double> status(originalPhase, E);  //-basePhase, E);
        return status;
    }
 
    double Phimax            = initialPhase;
    double phi               = initialPhase;
    double dphi              = Physics::pi / 360.0;
    const int numRefinements = Options::autoPhase;
    int j                    = -1;
 
    double E    = track(t, dt, phi);
    double Emax = E;
 
    do {
        j++;
        Emax         = E;
        initialPhase = phi;
        phi -= dphi;
        E = track(t, dt, phi);
    } while (E > Emax);
 
    if (j == 0) {
        phi = initialPhase;
        E   = Emax;
        // j = -1;
        do {
            // j ++;
            Emax         = E;
            initialPhase = phi;
            phi += dphi;
            E = track(t, dt, phi);
        } while (E > Emax);
    }
 
    for (int rl = 0; rl < numRefinements; ++rl) {
        dphi /= 2.;
        phi = initialPhase - dphi;
        E   = track(t, dt, phi);
        if (E > Emax) {
            initialPhase = phi;
            Emax         = E;
        } else {
            phi = initialPhase + dphi;
            E   = track(t, dt, phi);
            if (E > Emax) {
                initialPhase = phi;
                Emax         = E;
            }
        }
    }
    Phimax = std::fmod(initialPhase + Physics::two_pi, Physics::two_pi);
    E      = track(t, dt, Phimax + originalPhase);
    std::pair<double, double> status(Phimax, E);
 
    return status;
}
 
double CavityAutophaser::track(
    double t, const double dt, const double phase, std::ofstream* out) const {
    const Vector_t<double, 3>& refP = initialP_m;
 
    RFCavity* rfc       = static_cast<RFCavity*>(itsCavity_m.get());
    double initialPhase = rfc->getPhasem();
    rfc->setPhasem(phase);
 
    std::pair<double, double> pe = rfc->trackOnAxisParticle(
        refP(2), t, dt, itsReference_m.getQ(), itsReference_m.getM() * Units::eV2MeV, out);
    rfc->setPhasem(initialPhase);
 
    double finalKineticEnergy = Util::getKineticEnergy(
        Vector_t<double, 3>(0.0, 0.0, pe.first), itsReference_m.getM() * Units::eV2MeV);
    return finalKineticEnergy;
}