MultipoleTStraight.h

Go to the documentation of this file.

//
// Cubic Spline Interpolation to replace GSL spline
//
// Copyright (c) 2023, Paul Scherrer Institute, 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 License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
 
#ifndef ABSBEAMLINE_MULTIPOLET_STRAIGHT_H
#define ABSBEAMLINE_MULTIPOLET_STRAIGHT_H
 
#include "AbsBeamline/MultipoleT.h"
#include "AbsBeamline/MultipoleTConfig.h"
#include "BeamlineGeometry/StraightGeometry.h"
 
class MultipoleTStraight final : public MultipoleTBase {
public:
    explicit MultipoleTStraight(MultipoleT* element);
    void initialise() override;
    BGeometryBase* getGeometry() override { return &straightGeometry_m; }
    // Container-agnostic: R/E/B views come from the caller's particle container.
    void getField(
            Kokkos::View<Vector_t<double, 3>*> R, Kokkos::View<Vector_t<double, 3>*> E,
            Kokkos::View<Vector_t<double, 3>*> B, double scaling, size_t count) override;
    bool getField(
            const Vector_t<double, 3>& R, Vector_t<double, 3>& E, Vector_t<double, 3>& B,
            double scaling) override;
 
private:
    StraightGeometry straightGeometry_m;
 
    // Helpers
    KOKKOS_INLINE_FUNCTION static Vector_t<double, 3> toMagnetCoords(
            const Vector_t<double, 3>& R, const MultipoleTConfig& config);
    template <class ViewType>
    KOKKOS_INLINE_FUNCTION static bool computeBField(
            const Vector_t<double, 3>& R, Vector_t<double, 3>& B, double scaling,
            const MultipoleTConfig& config, const ViewType& tanhCoefficients);
};
 
KOKKOS_INLINE_FUNCTION Vector_t<double, 3> MultipoleTStraight::toMagnetCoords(
        const Vector_t<double, 3>& R, const MultipoleTConfig& config) {
    // Skew and entry angle
    Vector_t<double, 3> result = rotateFrame(R, config);
    // Magnet origin at the center rather than entry
    result[2] -= config.length_m / 2.0;
    return result;
}
 
template <class ViewType>
KOKKOS_INLINE_FUNCTION bool MultipoleTStraight::computeBField(
        const Vector_t<double, 3>& R, Vector_t<double, 3>& B, const double scaling,
        const MultipoleTConfig& config, const ViewType& tanhCoefficients) {
    const Vector_t<double, 3> RPrime = toMagnetCoords(R, config);
    const bool insideAperture        = Kokkos::abs(RPrime[1]) <= config.verticalAperture_m / 2.0
                                && Kokkos::abs(RPrime[0]) <= config.horizontalAperture_m / 2.0;
    const bool insideBoundingBox = config.boundingBoxLength_m == 0.0
                                   || Kokkos::abs(RPrime[2]) <= config.boundingBoxLength_m / 2.0;
    if (insideAperture && insideBoundingBox) {
        Vector_t<double, 3> myB{};
        Kokkos::Array<double, MaxDerivatives> dt;
        Kokkos::Array<double, MaxDerivatives> ds;
        calcTransverseDerivatives(
                config.transverseProfile_m, config.maxFOrder_m * 2 + 2, RPrime[0], dt);
        calcFringeDerivatives(
                config.fringeS0_m, config.fringeLambdaLeft_m, config.fringeLambdaRight_m, RPrime[2],
                tanhCoefficients, ds);
        for (unsigned int n = 0; n <= config.maxFOrder_m; n++) {
            double innerSumX{};
            double innerSumZ{};
            double innerSumS{};
            for (unsigned int i = 0; i <= n; i++) {
                const double k = factorial(n) / (factorial(i) * factorial(n - i));
                innerSumX += k * dt[2 * i + 1] * ds[2 * n - 2 * i];
                innerSumZ += k * dt[2 * i] * ds[2 * n - 2 * i];
                innerSumS += k * dt[2 * i] * ds[2 * n - 2 * i + 1];
            }
            const double negOnePowN = powerInteger(-1.0, n);
            const double xszk =
                    powerInteger(RPrime[1], 2 * n + 1) / factorial(2 * n + 1) * negOnePowN;
            const double zzk = powerInteger(RPrime[1], 2 * n) / factorial(2 * n) * negOnePowN;
            myB[0] += innerSumX * xszk;
            myB[1] += innerSumZ * zzk;
            myB[2] += innerSumS * xszk;
        }
        B += myB * scaling;
    }
    return !insideAperture;
}
 
#endif