Solenoid.cpp

Go to the documentation of this file.

// ------------------------------------------------------------------------
// $RCSfile: Solenoid.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: Solenoid
//   Defines the abstract interface for a solenoid magnet.
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 20.1.2026 $
// $Author: PSI $
//
// ------------------------------------------------------------------------
 
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/BeamlineVisitor.h"
#include "Fields/Fieldmap.h"
#include "PartBunch/PartBunch.h"
#include "Physics/Physics.h"
 
#include <cmath>
#include <fstream>
#include <iostream>
 
extern Inform* gmsg;
 
/* ============================== Constructors ============================== */
Solenoid::Solenoid() : Solenoid("") {}
 
Solenoid::Solenoid(const Solenoid& right)
    : Component(right),
      filename_m(right.filename_m),
      fieldmap_m(right.fieldmap_m),
      scale_m(right.scale_m),
      scaleError_m(right.scaleError_m),
      startField_m(right.startField_m),
      endField_m(right.endField_m),
      fast_m(right.fast_m) {}
 
Solenoid::Solenoid(const std::string& name)
    : Component(name),
      filename_m(""),
      fieldmap_m(nullptr),
      scale_m(1.0),
      scaleError_m(0.0),
      startField_m(0.0),
      endField_m(0.0),
      fast_m(true) {}
 
Solenoid::~Solenoid() {
    //    _Fieldmap::deleteFieldmap(filename_m);
}
/* ========================================================================== */
/* ============================== Apply Functions =========================== */
bool Solenoid::apply(const std::shared_ptr<ParticleContainer_t>& pc) {
    Inform m("Solenoid::apply");
    m << level5 << "Solenoid::apply() called." << endl;
 
    fieldmap_m->applyField(pc, scale_m + scaleError_m);
 
    return false;
}
 
bool Solenoid::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, P, t, E, B);
}
 
bool Solenoid::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) {
    if (R(2) >= startField_m && R(2) < endField_m) {
        Vector_t<double, 3> tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
 
        const bool outOfBounds = fieldmap_m->getFieldstrength(R, tmpE, tmpB);
        if (outOfBounds) {
            return getFlagDeleteOnTransverseExit();
        }
 
        B += (scale_m + scaleError_m) * tmpB;
    }
 
    return false;
}
 
bool Solenoid::applyToReferenceParticle(
        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) {
    if (R(2) >= startField_m && R(2) < endField_m) {
        Vector_t<double, 3> tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
 
        const bool outOfBounds = fieldmap_m->getFieldstrength(R, tmpE, tmpB);
        if (outOfBounds) return true;
 
        B += scale_m * tmpB;
    }
 
    return false;
}
/* ========================================================================== */
/* ============================== Functions ================================= */
void Solenoid::accept(BeamlineVisitor& visitor) const { visitor.visitSolenoid(*this); }
 
void Solenoid::setKS(double ks) { scale_m = ks; }
 
void Solenoid::setDKS(double ks) { scaleError_m = ks; }
 
void Solenoid::initialise(PartBunch_t* bunch, double& startField, double& endField) {
    Inform msg("Solenoid ", *gmsg);
 
    RefPartBunch_m = bunch;
 
    fieldmap_m = Fieldmap::getFieldmap(filename_m, fast_m);
 
    if (fieldmap_m != nullptr) {
        msg << level2 << getName() << " using file ";
        fieldmap_m->getInfo(&msg);
 
        double zBegin = 0.0, zEnd = 0.0;
        fieldmap_m->getFieldDimensions(zBegin, zEnd);
 
        startField_m = zBegin;
        endField_m   = zEnd;
 
        const double bodyBegin = startField;
        startField             = bodyBegin + startField_m;
        endField               = bodyBegin + endField_m;
    } else {
        startField_m = 0.0;
        endField_m   = 0.0;
        endField     = startField;
    }
}
 
void Solenoid::finalise() {}
 
bool Solenoid::bends() const { return false; }
 
void Solenoid::goOnline(const double&) {
    Fieldmap::readMap(filename_m);
    online_m = true;
}
 
void Solenoid::goOffline() {
    Fieldmap::freeMap(filename_m);
    online_m = false;
}
 
void Solenoid::setFieldMapFN(std::string fn) { filename_m = fn; }
 
void Solenoid::setFast(bool fast) { fast_m = fast; }
 
bool Solenoid::getFast() const { return fast_m; }
 
void Solenoid::getFieldExtend(double& zBegin, double& zEnd) const {
    zBegin = startField_m;
    zEnd   = endField_m;
}
 
ElementType Solenoid::getType() const { return ElementType::SOLENOID; }
 
bool Solenoid::isInside(const Vector_t<double, 3>& r) const {
    return fieldmap_m != nullptr && isInsideTransverse(r) && fieldmap_m->isInside(r);
}
 
void Solenoid::getElementDimensions(double& begin, double& end) const {
    begin = 0.0;
    end   = getElementLength();
}
 
bool Solenoid::getSupportEnvelope(double& horizontalRadius, double& verticalRadius) const {
    const auto aperture = getAperture();
    if (aperture.second.size() >= 2 && std::abs(aperture.second[0]) < 1e5
        && std::abs(aperture.second[1]) < 1e5) {
        horizontalRadius = std::abs(aperture.second[0]);
        verticalRadius   = std::abs(aperture.second[1]);
        return horizontalRadius > 0.0 && verticalRadius > 0.0;
    }
 
    if (fieldmap_m == nullptr) {
        return false;
    }
 
    try {
        double xIni = 0.0, xFinal = 0.0, yIni = 0.0, yFinal = 0.0, zIni = 0.0, zFinal = 0.0;
        fieldmap_m->getFieldDimensions(xIni, xFinal, yIni, yFinal, zIni, zFinal);
        horizontalRadius = 0.5 * std::abs(xFinal - xIni);
        verticalRadius   = 0.5 * std::abs(yFinal - yIni);
        return horizontalRadius > 0.0 && verticalRadius > 0.0;
    } catch (...) {
        return false;
    }
}