Multipole.cpp

Go to the documentation of this file.

//
// Class Multipole
//   The MULTIPOLE element defines a thick multipole.
//
// Copyright (c) 2012-2021, Paul Scherrer Institut, Villigen PSI, Switzerland
// 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 "AbsBeamline/Multipole.h"
 
#include "AbsBeamline/BeamlineVisitor.h"
#include "Fields/Fieldmap.h"
#include "PartBunch/PartBunch.h"
#include "Physics/Physics.h"
#include "Utilities/GeneralClassicException.h"
 
namespace {
    enum { DIPOLE, QUADRUPOLE, SEXTUPOLE, OCTUPOLE, DECAPOLE };
}
 
namespace {
    unsigned int factorial(unsigned int n) {
        static int fact[] = {1,    1,     2,      6,       24,       120,      720,
                             5040, 40320, 362880, 3628800, 39916800, 479001600};
        if (n > 12) {
            throw GeneralClassicException(
                "Multipole::factorial", "Factorial can't be larger than 12");
        }
        return fact[n];
    }
}  // namespace
 
Multipole::Multipole() : Multipole("") {
}
 
Multipole::Multipole(const Multipole& right)
    : Component(right),
      NormalComponents(right.NormalComponents),
      NormalComponentErrors(right.NormalComponentErrors),
      SkewComponents(right.SkewComponents),
      SkewComponentErrors(right.SkewComponentErrors),
      max_SkewComponent_m(right.max_SkewComponent_m),
      max_NormalComponent_m(right.max_NormalComponent_m),
      nSlices_m(right.nSlices_m) {
}
 
Multipole::Multipole(const std::string& name)
    : Component(name),
      NormalComponents(1, 0.0),
      NormalComponentErrors(1, 0.0),
      SkewComponents(1, 0.0),
      SkewComponentErrors(1, 0.0),
      max_SkewComponent_m(1),
      max_NormalComponent_m(1),
      nSlices_m(1) {
}
 
Multipole::~Multipole() {
}
 
void Multipole::accept(BeamlineVisitor& visitor) const {
    visitor.visitMultipole(*this);
}
 
double Multipole::getNormalComponent(int n) const {
    if (n < max_NormalComponent_m) {
        return NormalComponents[n];
    }
    return 0.0;
}
 
double Multipole::getSkewComponent(int n) const {
    if (n < max_SkewComponent_m) {
        return SkewComponents[n];
    }
    return 0.0;
}
 
void Multipole::setNormalComponent(int n, double v, double vError) {
    //   getField().setNormalComponent(n, v);
    PAssert_GE(n, 1);
 
    if (n > max_NormalComponent_m) {
        max_NormalComponent_m = n;
        NormalComponents.resize(max_NormalComponent_m, 0.0);
        NormalComponentErrors.resize(max_NormalComponent_m, 0.0);
    }
    switch (n - 1) {
        case DIPOLE:
            NormalComponents[n - 1]      = (v + vError) / 2;
            NormalComponentErrors[n - 1] = NormalComponents[n - 1];
            break;
        default:
            NormalComponents[n - 1]      = (v + vError) / factorial(n - 1);
            NormalComponentErrors[n - 1] = vError / factorial(n - 1);
    }
}
 
void Multipole::setSkewComponent(int n, double v, double vError) {
    //   getField().setSkewComponent(n, v);
    PAssert_GT(n, 1);
 
    if (n > max_SkewComponent_m) {
        max_SkewComponent_m = n;
        SkewComponents.resize(max_SkewComponent_m, 0.0);
        SkewComponentErrors.resize(max_SkewComponent_m, 0.0);
    }
    switch (n - 1) {
        case DIPOLE:
            SkewComponents[n - 1]      = (v + vError) / 2;
            SkewComponentErrors[n - 1] = SkewComponents[n - 1];
            break;
        default:
            SkewComponents[n - 1]      = (v + vError) / factorial(n - 1);
            SkewComponentErrors[n - 1] = vError / factorial(n - 1);
    }
}
 
// set the number of slices for map tracking
void Multipole::setNSlices(const std::size_t& nSlices) {
    nSlices_m = nSlices;
}
 
// get the number of slices for map tracking
std::size_t Multipole::getNSlices() const {
    return nSlices_m;
}
 
void Multipole::computeField(
    Vector_t<double, 3> R, Vector_t<double, 3>& /*E*/, Vector_t<double, 3>& B) {
    {
        std::vector<Vector_t<double, 3>> Rn(max_NormalComponent_m + 1);
        std::vector<double> fact(max_NormalComponent_m + 1);
        Rn[0]   = Vector_t<double, 3>(1.0);
        fact[0] = 1;
        for (int i = 0; i < max_NormalComponent_m; ++i) {
            switch (i) {
                case DIPOLE:
                    B(1) += NormalComponents[i];
                    break;
 
                case QUADRUPOLE:
                    B(0) += NormalComponents[i] * R(1);
                    B(1) += NormalComponents[i] * R(0);
                    break;
 
                case SEXTUPOLE:
                    B(0) += 2 * NormalComponents[i] * R(0) * R(1);
                    B(1) += NormalComponents[i] * (Rn[2](0) - Rn[2](1));
                    break;
 
                case OCTUPOLE:
                    B(0) += NormalComponents[i] * (3 * Rn[2](0) * Rn[1](1) - Rn[3](1));
                    B(1) += NormalComponents[i] * (Rn[3](0) - 3 * Rn[1](0) * Rn[2](1));
                    break;
 
                case DECAPOLE:
                    B(0) += 4 * NormalComponents[i] * (Rn[3](0) * Rn[1](1) - Rn[1](0) * Rn[3](1));
                    B(1) += NormalComponents[i] * (Rn[4](0) - 6 * Rn[2](0) * Rn[2](1) + Rn[4](1));
                    break;
 
                default: {
                    double powMinusOne = 1;
                    double Bx = 0.0, By = 0.0;
                    for (int j = 1; j <= (i + 1) / 2; ++j) {
                        Bx += powMinusOne * NormalComponents[i]
                              * (Rn[i - 2 * j + 1](0) * fact[i - 2 * j + 1] * Rn[2 * j - 1](1)
                                 * fact[2 * j - 1]);
                        By += powMinusOne * NormalComponents[i]
                              * (Rn[i - 2 * j + 2](0) * fact[i - 2 * j + 2] * Rn[2 * j - 2](1)
                                 * fact[2 * j - 2]);
                        powMinusOne *= -1;
                    }
 
                    if ((i + 1) / 2 == i / 2) {
                        int j = (i + 2) / 2;
                        By += powMinusOne * NormalComponents[i]
                              * (Rn[i - 2 * j + 2](0) * fact[i - 2 * j + 2] * Rn[2 * j - 2](1)
                                 * fact[2 * j - 2]);
                    }
                    B(0) += Bx;
                    B(1) += By;
                }
            }
 
            Rn[i + 1](0) = Rn[i](0) * R(0);
            Rn[i + 1](1) = Rn[i](1) * R(1);
            fact[i + 1]  = fact[i] / (i + 1);
        }
    }
 
    {
        std::vector<Vector_t<double, 3>> Rn(max_SkewComponent_m + 1);
        std::vector<double> fact(max_SkewComponent_m + 1);
        Rn[0]   = Vector_t<double, 3>(1.0);
        fact[0] = 1;
        for (int i = 0; i < max_SkewComponent_m; ++i) {
            switch (i) {
                case DIPOLE:
                    B(0) -= SkewComponents[i];
                    break;
 
                case QUADRUPOLE:
                    B(0) -= SkewComponents[i] * R(0);
                    B(1) += SkewComponents[i] * R(1);
                    break;
 
                case SEXTUPOLE:
                    B(0) -= SkewComponents[i] * (Rn[2](0) - Rn[2](1));
                    B(1) += 2 * SkewComponents[i] * R(0) * R(1);
                    break;
 
                case OCTUPOLE:
                    B(0) -= SkewComponents[i] * (Rn[3](0) - 3 * Rn[1](0) * Rn[2](1));
                    B(1) += SkewComponents[i] * (3 * Rn[2](0) * Rn[1](1) - Rn[3](1));
                    break;
 
                case DECAPOLE:
                    B(0) -= SkewComponents[i] * (Rn[4](0) - 6 * Rn[2](0) * Rn[2](1) + Rn[4](1));
                    B(1) += 4 * SkewComponents[i] * (Rn[3](0) * Rn[1](1) - Rn[1](0) * Rn[3](1));
                    break;
 
                default: {
                    double powMinusOne = 1;
                    double Bx = 0, By = 0;
                    for (int j = 1; j <= (i + 1) / 2; ++j) {
                        Bx -= powMinusOne * SkewComponents[i]
                              * (Rn[i - 2 * j + 2](0) * fact[i - 2 * j + 2] * Rn[2 * j - 2](1)
                                 * fact[2 * j - 2]);
                        By += powMinusOne * SkewComponents[i]
                              * (Rn[i - 2 * j + 1](0) * fact[i - 2 * j + 1] * Rn[2 * j - 1](1)
                                 * fact[2 * j - 1]);
                        powMinusOne *= -1;
                    }
 
                    if ((i + 1) / 2 == i / 2) {
                        int j = (i + 2) / 2;
                        Bx -= powMinusOne * SkewComponents[i]
                              * (Rn[i - 2 * j + 2](0) * fact[i - 2 * j + 2] * Rn[2 * j - 2](1)
                                 * fact[2 * j - 2]);
                    }
 
                    B(0) += Bx;
                    B(1) += By;
                }
            }
 
            Rn[i + 1](0) = Rn[i](0) * R(0);
            Rn[i + 1](1) = Rn[i](1) * R(1);
            fact[i + 1]  = fact[i] / (i + 1);
        }
    }
}
 
bool Multipole::apply(
    const size_t& i, const double&, Vector_t<double, 3>& E, Vector_t<double, 3>& B) {
    //
    // \todo this is most inefficient
    std::shared_ptr<ParticleContainer_t> pc = RefPartBunch_m->getParticleContainer();
    auto Rview                              = pc->R.getView();
    const Vector_t<double, 3> R             = Rview(i);
 
    if (R(2) < 0.0 || R(2) > getElementLength())
        return false;
    if (!isInsideTransverse(R))
        return getFlagDeleteOnTransverseExit();
 
    Vector_t<double, 3> Ef(0.0), Bf(0.0);
    computeField(R, Ef, Bf);
 
    for (unsigned int d = 0; d < 3; ++d) {
        E[d] += Ef(d);
        B[d] += Bf(d);
    }
 
    return false;
}
 
bool Multipole::apply(
    const Vector_t<double, 3>& R, const Vector_t<double, 3>&, const double&, Vector_t<double, 3>& E,
    Vector_t<double, 3>& B) {
    if (R(2) < 0.0 || R(2) > getElementLength())
        return false;
    if (!isInsideTransverse(R))
        return getFlagDeleteOnTransverseExit();
 
    computeField(R, E, B);
 
    return false;
}
 
bool Multipole::applyToReferenceParticle(
    const Vector_t<double, 3>& R, const Vector_t<double, 3>&, const double&, Vector_t<double, 3>& E,
    Vector_t<double, 3>& B) {
    if (R(2) < 0.0 || R(2) > getElementLength())
        return false;
    if (!isInsideTransverse(R))
        return true;
 
    for (int i = 0; i < max_NormalComponent_m; ++i) {
        NormalComponents[i] -= NormalComponentErrors[i];
    }
    for (int i = 0; i < max_SkewComponent_m; ++i) {
        SkewComponents[i] -= SkewComponentErrors[i];
    }
 
    computeField(R, E, B);
 
    for (int i = 0; i < max_NormalComponent_m; ++i) {
        NormalComponents[i] += NormalComponentErrors[i];
    }
    for (int i = 0; i < max_SkewComponent_m; ++i) {
        SkewComponents[i] += SkewComponentErrors[i];
    }
 
    return false;
}
 
void Multipole::initialise(PartBunch_t* bunch, double& startField, double& endField) {
    RefPartBunch_m = bunch;
    endField       = startField + getElementLength();
    online_m       = true;
}
 
void Multipole::finalise() {
    online_m = false;
}
 
bool Multipole::bends() const {
    return false;
}
 
void Multipole::getDimensions(double& zBegin, double& zEnd) const {
    zBegin = 0.0;
    zEnd   = getElementLength();
}
 
ElementType Multipole::getType() const {
    return ElementType::MULTIPOLE;
}
 
bool Multipole::isInside(const Vector_t<double, 3>& r) const {
    if (r(2) >= 0.0 && r(2) < getElementLength()) {
        return isInsideTransverse(r);
    }
 
    return false;
}
 
bool Multipole::isFocusing(unsigned int component) const {
    if (component >= NormalComponents.size())
        throw GeneralClassicException("Multipole::isFocusing", "component too big");
 
    return NormalComponents[component] * std::pow(-1, component + 1)
               * RefPartBunch_m->getChargePerParticle()
           > 0.0;
}