MultipoleT.cpp

Go to the documentation of this file.

/*
 *  Copyright (c) 2017, Titus Dascalu
 *  Copyright (c) 2018, Martin Duy Tat
 *  All rights reserved.
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are met:
 *  1. Redistributions of source code must retain the above copyright notice,
 *     this list of conditions and the following disclaimer.
 *  2. Redistributions in binary form must reproduce the above copyright notice,
 *     this list of conditions and the following disclaimer in the documentation
 *     and/or other materials provided with the distribution.
 *  3. Neither the name of STFC nor the names of its contributors may be used to
 *     endorse or promote products derived from this software without specific
 *     prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 */
 
#include "AbsBeamline/Component.h"
#include "BeamlineGeometry/StraightGeometry.h"
#include "Fields/BMultipoleField.h"
 
#include <gsl/gsl_math.h>
#include "MultipoleT.h"
#include "MultipoleTFunctions/CoordinateTransform.h"
#include "MultipoleTFunctions/RecursionRelation.h"
#include "MultipoleTFunctions/RecursionRelationTwo.h"
#include "MultipoleTFunctions/tanhDeriv.h"
 
using namespace endfieldmodel;
 
MultipoleT::MultipoleT(const std::string& name)
    : Component(name),
      fringeField_l(endfieldmodel::Tanh()),
      fringeField_r(endfieldmodel::Tanh()),
      maxOrder_m(5),
      maxOrderX_m(10),
      transMaxOrder_m(0),
      planarArcGeometry_m(1., 1.),
      length_m(1.0),
      angle_m(0.0),
      entranceAngle_m(0.0),
      rotation_m(0.0),
      variableRadius_m(false),
      boundingBoxLength_m(0.0),
      verticalApert_m(0.5),
      horizApert_m(0.5) {
}
 
MultipoleT::MultipoleT(const MultipoleT& right)
    : Component(right),
      fringeField_l(right.fringeField_l),
      fringeField_r(right.fringeField_r),
      maxOrder_m(right.maxOrder_m),
      maxOrderX_m(right.maxOrderX_m),
      recursion_VarRadius_m(right.recursion_VarRadius_m),
      recursion_ConstRadius_m(right.recursion_ConstRadius_m),
      transMaxOrder_m(right.transMaxOrder_m),
      transProfile_m(right.transProfile_m),
      planarArcGeometry_m(right.planarArcGeometry_m),
      length_m(right.length_m),
      angle_m(right.angle_m),
      entranceAngle_m(right.entranceAngle_m),
      rotation_m(right.rotation_m),
      variableRadius_m(right.variableRadius_m),
      boundingBoxLength_m(right.boundingBoxLength_m),
      verticalApert_m(right.verticalApert_m),
      horizApert_m(right.horizApert_m),
      dummy() {
    RefPartBunch_m = right.RefPartBunch_m;
}
 
MultipoleT::~MultipoleT() {
}
 
ElementBase* MultipoleT::clone() const {
    MultipoleT* newMultipole = new MultipoleT(*this);
    newMultipole->initialise();
    return newMultipole;
}
 
void MultipoleT::finalise() {
    RefPartBunch_m = nullptr;
}
 
bool MultipoleT::apply(
    const Vector_t<double, 3>& R, const Vector_t<double, 3>& /*P*/, const double& /*t*/,
    Vector_t<double, 3>& /*E*/, Vector_t<double, 3>& B) {
    Vector_t<double, 3> R_prime = rotateFrame(R);
    // R_prime[2] *= -1; //OPAL uses different sign convention...
    Vector_t<double, 3> X = R_prime;
    R_prime               = transformCoords(X);
    if (insideAperture(R_prime)) {
        double theta;
        double prefactor = (length_m / angle_m)
                           * (tanh((fringeField_l.getX0()) / fringeField_l.getLambda())
                              + tanh((fringeField_r.getX0()) / fringeField_r.getLambda()));
        if (angle_m == 0.0) {
            theta = 0.0;
        } else if (!variableRadius_m) {
            theta = R_prime[2] * angle_m / length_m;
        } else {  // variableRadius_m == true
            theta =
                fringeField_l.getLambda()
                    * log(cosh((R_prime[2] + fringeField_l.getX0()) / fringeField_l.getLambda()))
                - fringeField_r.getLambda()
                      * log(cosh((R_prime[2] - fringeField_r.getX0()) / fringeField_r.getLambda()));
            theta /= prefactor;
        }
        double Bx = getBx(R_prime);
        double Bs = getBs(R_prime);
        B[0]      = Bx * cos(theta) - Bs * sin(theta);
        B[2]      = Bx * sin(theta) + Bs * cos(theta);
        B[1]      = getBz(R_prime);
        // B[2] *= -1; //OPAL uses different sign convention
        return false;
    } else {
        for (int i = 0; i < 3; i++) {
            B[i] = 0.0;
        }
        return getFlagDeleteOnTransverseExit();
    }
}
 
bool MultipoleT::apply(
    const size_t& i, const double& t, Vector_t<double, 3>& E, Vector_t<double, 3>& B) {
    std::shared_ptr<ParticleContainer_t> pc = RefPartBunch_m->getParticleContainer();
    auto Rview                              = pc->R.getView();
    auto Pview                              = pc->P.getView();
 
    const Vector_t<double, 3> R = Rview(i);
    const Vector_t<double, 3> P = Pview(i);
 
    return apply(R(i), P(i), t, E, B);
}
 
Vector_t<double, 3> MultipoleT::rotateFrame(const Vector_t<double, 3>& R) {
    Vector_t<double, 3> R_prime(3), R_pprime(3);
    // 1st rotation
    R_prime[0] = R[0] * cos(rotation_m) + R[1] * sin(rotation_m);
    R_prime[1] = -R[0] * sin(rotation_m) + R[1] * cos(rotation_m);
    R_prime[2] = R[2];
    // 2nd rotation
    R_pprime[0] = R_prime[2] * sin(entranceAngle_m) + R_prime[0] * cos(entranceAngle_m);
    R_pprime[1] = R_prime[1];
    R_pprime[2] = R_prime[2] * cos(entranceAngle_m) - R_prime[0] * sin(entranceAngle_m);
    return R_pprime;
}
 
Vector_t<double, 3> MultipoleT::rotateFrameInverse(Vector_t<double, 3>& B) {
    Vector_t<double, 3> B_prime(3);
    B_prime[0] = B[0] * cos(rotation_m) + B[1] * sin(rotation_m);
    B_prime[1] = -B[0] * sin(rotation_m) + B[1] * cos(rotation_m);
    B_prime[2] = B[2];
    return B_prime;
}
 
bool MultipoleT::insideAperture(const Vector_t<double, 3>& R) {
    if (std::abs(R[1]) <= verticalApert_m / 2. && std::abs(R[0]) <= horizApert_m / 2.) {
        return true;
    } else {
        return false;
    }
}
 
Vector_t<double, 3> MultipoleT::transformCoords(const Vector_t<double, 3>& R) {
    Vector_t<double, 3> X(3);
    // If radius is constant
    if (!variableRadius_m) {
        double radius = length_m / angle_m;
        // Transform from Cartesian to Frenet-Serret along the magnet
        double alpha = atan(R[2] / (R[0] + radius));
        if (alpha != 0.0 && angle_m != 0.0) {
            X[0] = R[2] / sin(alpha) - radius;
            X[1] = R[1];
            X[2] = radius * alpha;  // + boundingBoxLength_m;
        } else {
            X[0] = R[0];
            X[1] = R[1];
            X[2] = R[2];  // + boundingBoxLength_m;
        }
    } else {
        // If radius is variable
        coordinatetransform::CoordinateTransform t(
            R[0], R[1], R[2], fringeField_l.getX0(), fringeField_l.getLambda(),
            fringeField_r.getLambda(), (length_m / angle_m));
        std::vector<double> r = t.getTransformation();
        X[0]                  = r[0];
        X[1]                  = r[1];
        X[2]                  = r[2];
    }
    return X;
}
 
void MultipoleT::setMaxOrder(std::size_t maxOrder) {
    if (variableRadius_m && angle_m != 0.0) {
        std::size_t N = recursion_VarRadius_m.size();
        while (maxOrder >= N) {
            polynomial::RecursionRelationTwo r(N, 2 * (N + maxOrderX_m + 1));
            r.resizeX(transMaxOrder_m);
            r.truncate(maxOrderX_m);
            recursion_VarRadius_m.push_back(r);
            N = recursion_VarRadius_m.size();
        }
    } else if (!variableRadius_m && angle_m != 0.0) {
        std::size_t N = recursion_ConstRadius_m.size();
        while (maxOrder >= N) {
            polynomial::RecursionRelation r(N, 2 * (N + maxOrderX_m + 1));
            r.resizeX(transMaxOrder_m);
            r.truncate(maxOrderX_m);
            recursion_ConstRadius_m.push_back(r);
            N = recursion_ConstRadius_m.size();
        }
    }
    maxOrder_m = maxOrder;
}
 
void MultipoleT::setTransProfile(std::size_t n, double dTn) {
    if (n > transMaxOrder_m) {
        transMaxOrder_m = n;
        transProfile_m.resize(n + 1, 0.0);
    }
    transProfile_m[n] = dTn;
}
 
bool MultipoleT::setFringeField(double s0, double lambda_l, double lambda_r) {
    fringeField_l.Tanh::setLambda(lambda_l);
    fringeField_l.Tanh::setX0(s0);
    fringeField_l.Tanh::setTanhDiffIndices(2 * maxOrder_m + 1);
 
    fringeField_r.Tanh::setLambda(lambda_r);
    fringeField_r.Tanh::setX0(s0);
    fringeField_r.Tanh::setTanhDiffIndices(2 * maxOrder_m + 1);
 
    return true;
}
 
double MultipoleT::getBz(const Vector_t<double, 3>& R) {
    double Bz = 0.0;
    if (angle_m == 0.0) {
        // Straight geometry -> Use corresponding field expansion directly
        for (std::size_t n = 0; n <= maxOrder_m; n++) {
            double f_n = 0.0;
            for (std::size_t i = 0; i <= n; i++) {
                f_n += gsl_sf_choose(n, i) * getTransDeriv(2 * i, R[0])
                       * getFringeDeriv(2 * n - 2 * i, R[2]);
            }
            f_n *= gsl_sf_pow_int(-1.0, n);
            Bz += gsl_sf_pow_int(R[1], 2 * n) / gsl_sf_fact(2 * n) * f_n;
        }
    } else {
        if (variableRadius_m == true and getFringeDeriv(0, R[2]) < 1.0e-12) {
            // Return 0 if end of fringe field is reached
            // This is to avoid functions being called at infinite radius
            return 0.0;
        }
        // Curved geometry -> Use full machinery of differential operators
        for (std::size_t n = 0; n <= maxOrder_m; n++) {
            double f_n = getFn(n, R[0], R[2]);
            Bz += gsl_sf_pow_int(R[1], 2 * n) / gsl_sf_fact(2 * n) * f_n;
        }
    }
    return Bz;
}
 
double MultipoleT::getBx(const Vector_t<double, 3>& R) {
    double Bx = 0.0;
    if (angle_m == 0.0) {
        // Straight geometry -> Use corresponding field expansion directly
        for (std::size_t n = 0; n <= maxOrder_m; n++) {
            double f_n = 0.0;
            for (std::size_t i = 0; i <= n; i++) {
                f_n += gsl_sf_choose(n, i) * getTransDeriv(2 * i + 1, R[0])
                       * getFringeDeriv(2 * n - 2 * i, R[2]);
            }
            f_n *= gsl_sf_pow_int(-1.0, n);
            Bx += gsl_sf_pow_int(R[1], 2 * n + 1) / gsl_sf_fact(2 * n + 1) * f_n;
        }
    } else {
        if (variableRadius_m == true and getFringeDeriv(0, R[2]) < 1.0e-12) {
            // Return 0 if end of fringe field is reached
            // This is to avoid functions being called at infinite radius
            return 0.0;
        }
        // Curved geometry -> Use full machinery of differential operators
        for (std::size_t n = 0; n <= maxOrder_m; n++) {
            double partialX_fn = getFnDerivX(n, R[0], R[2]);
            Bx += partialX_fn * gsl_sf_pow_int(R[1], 2 * n + 1) / gsl_sf_fact(2 * n + 1);
        }
    }
    return Bx;
}
 
double MultipoleT::getBs(const Vector_t<double, 3>& R) {
    double Bs = 0.0;
    if (angle_m == 0.0) {
        // Straight geometry -> Use corresponding field expansion directly
        for (std::size_t n = 0; n <= maxOrder_m; n++) {
            double f_n = 0.0;
            for (std::size_t i = 0; i <= n; i++) {
                f_n += gsl_sf_choose(n, i) * getTransDeriv(2 * i, R[0])
                       * getFringeDeriv(2 * n - 2 * i + 1, R[2]);
            }
            f_n *= gsl_sf_pow_int(-1.0, n);
            Bs += gsl_sf_pow_int(R[1], 2 * n + 1) / gsl_sf_fact(2 * n + 1) * f_n;
        }
    } else {
        if (variableRadius_m == true and getFringeDeriv(0, R[2]) < 1.0e-12) {
            // Return 0 if end of fringe field is reached
            // This is to avoid functions being called at infinite radius
            return 0.0;
        }
        // Curved geometry -> Use full machinery of differential operators
        for (std::size_t n = 0; n <= maxOrder_m; n++) {
            double partialS_fn = getFnDerivS(n, R[0], R[2]);
            Bs += partialS_fn * gsl_sf_pow_int(R[1], 2 * n + 1) / gsl_sf_fact(2 * n + 1);
        }
        Bs /= getScaleFactor(R[0], R[2]);
    }
    return Bs;
}
 
double MultipoleT::getFringeDeriv(int n, double s) {
    if (n <= 10) {
        return (fringeField_l.getTanh(s, n) - fringeField_r.getNegTanh(s, n)) / 2;
    } else {
        return tanhderiv::integrate(
            s, fringeField_l.Tanh::getX0(), fringeField_l.Tanh::getLambda(),
            fringeField_r.Tanh::getLambda(), n);
    }
}
 
double MultipoleT::getTransDeriv(std::size_t n, double x) {
    double func              = 0;
    std::vector<double> temp = transProfile_m;
    if (n <= transMaxOrder_m) {
        if (n != 0) {
            for (std::size_t i = 1; i <= n; i++) {
                for (std::size_t j = 0; j <= transMaxOrder_m; j++) {
                    if (j <= transMaxOrder_m - i) {
                        // Move terms to the left and multiply by power
                        temp[j] = temp[j + 1] * (j + 1);
                    } else {
                        // Put 0 at the end for missing higher powers
                        temp[j] = 0.0;
                    }
                }
            }
        }
        // Now use the vector to calculate value of the function
        for (std::size_t k = 0; k <= transMaxOrder_m; k++) {
            func += temp[k] * gsl_sf_pow_int(x, k);
        }
    }
    return func;
}
 
void MultipoleT::setDipoleConstant(double B0) {
    if (transMaxOrder_m < 1) {
        transProfile_m.resize(1, 0.);
    }
    transProfile_m[0] = B0;
}
 
void MultipoleT::accept(BeamlineVisitor& visitor) const {
    visitor.visitMultipoleT(*this);
}
 
void MultipoleT::getDimensions(double& /*zBegin*/, double& /*zEnd*/) const {
}
 
void MultipoleT::setAperture(double vertAp, double horizAp) {
    verticalApert_m = vertAp;
    horizApert_m    = horizAp;
}
 
std::vector<double> MultipoleT::getAperture() const {
    std::vector<double> temp(2, 0.0);
    temp[0] = verticalApert_m;
    temp[1] = horizApert_m;
    return temp;
}
 
std::vector<double> MultipoleT::getFringeLength() const {
    std::vector<double> temp(2, 0.0);
    temp[0] = fringeField_l.getLambda();
    temp[1] = fringeField_r.getLambda();
    return temp;
}
 
double MultipoleT::getRadius(double s) {
    double centralRadius = length_m / angle_m;
    if (!variableRadius_m) {
        return centralRadius;
    }
    if (getFringeDeriv(0, s) != 0.0) {
        return centralRadius * getFringeDeriv(0, 0) / getFringeDeriv(0, s);
    } else {
        return -1;  // Return -1 if radius is infinite
    }
}
 
double MultipoleT::getScaleFactor(double x, double s) {
    if (!variableRadius_m) {
        return 1 + x * angle_m / length_m;
    } else {
        return 1 + x / getRadius(s);
    }
}
 
double MultipoleT::getFnDerivX(std::size_t n, double x, double s) {
    if (n == 0) {
        return getTransDeriv(1, x) * getFringeDeriv(0, s);
    }
    double deriv    = 0.0;
    double stepSize = 1e-3;
    deriv += 1. * getFn(n, x - 2. * stepSize, s);
    deriv += -8. * getFn(n, x - stepSize, s);
    deriv += 8. * getFn(n, x + stepSize, s);
    deriv += -1. * getFn(n, x + 2. * stepSize, s);
    deriv /= 12 * stepSize;
    return deriv;
}
 
double MultipoleT::getFnDerivS(std::size_t n, double x, double s) {
    if (n == 0) {
        return getTransDeriv(0, x) * getFringeDeriv(1, s);
    }
    double deriv    = 0.0;
    double stepSize = 1e-3;
    deriv += 1. * getFn(n, x, s - 2. * stepSize);
    deriv += -8. * getFn(n, x, s - stepSize);
    deriv += 8. * getFn(n, x, s + stepSize);
    deriv += -1. * getFn(n, x, s + 2. * stepSize);
    deriv /= 12 * stepSize;
    return deriv;
}
 
double MultipoleT::getFn(std::size_t n, double x, double s) {
    if (n == 0) {
        return getTransDeriv(0, x) * getFringeDeriv(0, s);
    }
    if (!variableRadius_m) {
        double rho  = length_m / angle_m;
        double func = 0.0;
 
        for (std::size_t j = 0; j <= recursion_ConstRadius_m.at(n).getMaxSDerivatives(); j++) {
            double FringeDerivj = getFringeDeriv(2 * j, s);
            for (std::size_t i = 0; i <= recursion_ConstRadius_m.at(n).getMaxXDerivatives(); i++) {
                if (recursion_ConstRadius_m.at(n).isPolynomialZero(i, j)) {
                    continue;
                }
                func += (recursion_ConstRadius_m.at(n).evaluatePolynomial(x / rho, i, j)
                         * getTransDeriv(i, x) * FringeDerivj)
                        / gsl_sf_pow_int(rho, 2 * n - i - 2 * j);
            }
        }
        func *= gsl_sf_pow_int(-1.0, n);
        return func;
    } else {
        double rho  = length_m / angle_m;
        double S_0  = getFringeDeriv(0, 0);
        double y    = getFringeDeriv(0, s) / (S_0 * rho);
        double func = 0.0;
        std::vector<double> fringeDerivatives;
        for (std::size_t j = 0; j <= recursion_VarRadius_m.at(n).getMaxSDerivatives(); j++) {
            fringeDerivatives.push_back(getFringeDeriv(j, s) / (S_0 * rho));
        }
        for (std::size_t i = 0; i <= recursion_VarRadius_m.at(n).getMaxXDerivatives(); i++) {
            double temp = 0.0;
            for (std::size_t j = 0; j <= recursion_VarRadius_m.at(n).getMaxSDerivatives(); j++) {
                temp +=
                    recursion_VarRadius_m.at(n).evaluatePolynomial(x, y, i, j, fringeDerivatives)
                    * fringeDerivatives.at(j);
            }
            func += temp * getTransDeriv(i, x);
        }
        func *= gsl_sf_pow_int(-1.0, n) * S_0 * rho;
        return func;
    }
}
 
void MultipoleT::initialise() {
    planarArcGeometry_m.setElementLength(2 * boundingBoxLength_m);
    planarArcGeometry_m.setCurvature(angle_m / length_m);
}
 
void MultipoleT::initialise(PartBunch_t* bunch, double& /*startField*/, double& /*endField*/) {
    RefPartBunch_m = bunch;
    initialise();
}
 
bool MultipoleT::bends() const {
    return (transProfile_m[0] != 0);
}
 
PlanarArcGeometry& MultipoleT::getGeometry() {
    return planarArcGeometry_m;
}
 
const PlanarArcGeometry& MultipoleT::getGeometry() const {
    return planarArcGeometry_m;
}
 
EMField& MultipoleT::getField() {
    return dummy;
}
 
const EMField& MultipoleT::getField() const {
    return dummy;
}