VerticalFFAMagnet.cpp

Go to the documentation of this file.

//
// Source file for VerticalFFAMagnet Component
//
// Copyright (c) 2019 Chris Rogers
// All rights reserved.
//
// OPAL is licensed under GNU GPL version 3.
//
 
#include "AbsBeamline/BeamlineVisitor.h"
#include "AbsBeamline/EndFieldModel/EndFieldModel.h"
 
#include "AbsBeamline/VerticalFFAMagnet.h"
 
#include <cmath>
 
VerticalFFAMagnet::VerticalFFAMagnet(const std::string &name)
    : Component(name), straightGeometry_m(1.) {
}
 
VerticalFFAMagnet::VerticalFFAMagnet(const VerticalFFAMagnet &right) 
    : Component(right),
    straightGeometry_m(right.straightGeometry_m),
    dummy(right.dummy),
    maxOrder_m(right.maxOrder_m),
    k_m(right.k_m),
    Bz_m(right.Bz_m),
    zNegExtent_m(right.zNegExtent_m),
    zPosExtent_m(right.zPosExtent_m),
    halfWidth_m(right.halfWidth_m),
    bbLength_m(right.bbLength_m),
    endField_m(right.endField_m->clone()),
    dfCoefficients_m(right.dfCoefficients_m) {
        RefPartBunch_m = right.RefPartBunch_m;
}
 
 
VerticalFFAMagnet::~VerticalFFAMagnet() {
}
 
ElementBase* VerticalFFAMagnet::clone() const {
    VerticalFFAMagnet* magnet = new VerticalFFAMagnet(*this);
    magnet->initialise();
    return magnet;
}
 
EMField &VerticalFFAMagnet::getField() {
    return dummy;
}
 
const EMField &VerticalFFAMagnet::getField() const {
    return dummy;
}
 
void VerticalFFAMagnet::initialise() {
    calculateDfCoefficients();
    straightGeometry_m.setElementLength(bbLength_m); // length = phi r
}
 
void VerticalFFAMagnet::initialise(PartBunch_t *bunch, double &/*startField*/, double &/*endField*/) {
    RefPartBunch_m = bunch;
    initialise();
}
 
void VerticalFFAMagnet::finalise() {
    RefPartBunch_m = nullptr;
}
 
BGeometryBase& VerticalFFAMagnet::getGeometry() {
    return straightGeometry_m;
}
 
const BGeometryBase& VerticalFFAMagnet::getGeometry() const {
    return straightGeometry_m;
}
 
void VerticalFFAMagnet::accept(BeamlineVisitor& visitor) const {
    visitor.visitVerticalFFAMagnet(*this);
}
 
 
bool VerticalFFAMagnet::getFieldValue(const Vector_t<double, 3> &R, Vector_t<double, 3> &B) const {
    if (std::abs(R[0]) > halfWidth_m ||
        R[2] < 0. || R[2] > bbLength_m ||
        R[1] < -zNegExtent_m || R[1] > zPosExtent_m) {
        return true;
    }
    std::vector<double> fringeDerivatives(maxOrder_m+2, 0.);
    double zRel = R[2]-bbLength_m/2.; // z relative to centre of magnet
    for (size_t i = 0; i < fringeDerivatives.size(); ++i) {
        fringeDerivatives[i] = endField_m->function(zRel, i); // d^i_phi f
    }
 
    std::vector<double> x_n(maxOrder_m+1); // x^n
    x_n[0] = 1.; // x^0
    for (size_t i = 1; i < x_n.size(); ++i) {
        x_n[i] = x_n[i-1]*R[0];
    }
 
    // note that the last element is always 0, because dfCoefficients_m is 
    // of size maxOrder_m+1. This leads to better Maxwellianness in testing.
    std::vector<double> f_n(maxOrder_m+2, 0.);
    std::vector<double> dz_f_n(maxOrder_m+1, 0.);
    for (size_t n = 0; n < dfCoefficients_m.size(); ++n) {
        const std::vector<double>& coefficients = dfCoefficients_m[n];
        for (size_t i = 0; i < coefficients.size(); ++i) {
            f_n[n] += coefficients[i]*fringeDerivatives[i];
            dz_f_n[n] += coefficients[i]*fringeDerivatives[i+1];
        }
    }
    double bref = Bz_m*exp(k_m*R[1]);
    B[0] = 0.;
    B[1] = 0.;
    B[2] = 0.;
    for (size_t n = 0; n < x_n.size(); ++n) {
        B[0] += bref*f_n[n+1]*(n+1)/k_m*x_n[n];
        B[1] += bref*f_n[n]*x_n[n];
        B[2] += bref*dz_f_n[n]/k_m*x_n[n];
    }
    return false;
}
 
void VerticalFFAMagnet::calculateDfCoefficients() {
    dfCoefficients_m = std::vector< std::vector<double> >(maxOrder_m+1);
    dfCoefficients_m[0] = std::vector<double>(1, 1.);
    if (maxOrder_m > 0) {
        dfCoefficients_m[1] = std::vector<double>();
    }
    // n indexes like the polynomial order of the midplane expansion
    // e.g. Bz = exp(mz) f_n y^n
    // where y is distance from the midplane and z is height
    for (size_t n = 2; n < dfCoefficients_m.size(); n+=2) {
        const std::vector<double>& oldCoefficients = dfCoefficients_m[n-2];
        std::vector<double> coefficients(oldCoefficients.size()+2, 0);
        // j indexes the derivative of f_0
        for (size_t j = 0; j < oldCoefficients.size(); ++j) {
            coefficients[j] += -1./(n)/(n-1)*k_m*k_m*oldCoefficients[j];
            coefficients[j+2] += -1./(n)/(n-1)*oldCoefficients[j];
        }
        dfCoefficients_m[n] = coefficients;
    }
}
 
void VerticalFFAMagnet::setEndField(endfieldmodel::EndFieldModel* endField) {
    endField_m.reset(endField);
    endField_m->setMaximumDerivative(maxOrder_m);
}
 
void VerticalFFAMagnet::setMaxOrder(size_t maxOrder) {
    if (endField_m.get()) {
        endField_m->setMaximumDerivative(maxOrder);
    }
    maxOrder_m = maxOrder;
}