OpalTrimCoil.cpp

Go to the documentation of this file.

//
// Class OpalTrimCoil
//   A TRIMCOIL definition is used to define a trim coil which can be applied
//   to a Cyclotron.
//
// Copyright (c) 2018 - 2019, Matthias Frey and Jochem Snuverink,
//                            Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// Implemented as part of the PhD thesis
// "Precise Simulations of Multibunches in High Intensity Cyclotrons"
// and the paper
// "Matching of turn pattern measurements for cyclotrons using multiobjective optimization"
// (https://doi.org/10.1103/PhysRevAccelBeams.22.064602)
//
// 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 "TrimCoils/OpalTrimCoil.h"

#include "AbstractObjects/OpalData.h"
#include "Attributes/Attributes.h"
#include "TrimCoils/TrimCoilBFit.h"
#include "TrimCoils/TrimCoilMirrored.h"
#include "TrimCoils/TrimCoilPhaseFit.h"
#include "Utilities/OpalException.h"
#include "Utility/Inform.h"

#include <map>

extern Inform *gmsg;

// The attributes of class OpalTrimCoil.
namespace {
    enum {
        TYPE,       // The type of trim coil
        COEFNUM,
        COEFDENOM,
        COEFNUMPHI,
        COEFDENOMPHI,
        BMAX,
        PHIMIN,
        PHIMAX,
        RMIN,
        RMAX,
        SLPTC,
        SIZE
    };
}

OpalTrimCoil::OpalTrimCoil():
    Definition(SIZE, "TRIMCOIL",
               "The \"TRIMCOIL\" statement defines a trim coil."),
    trimcoil_m(nullptr) {
    itsAttr[TYPE] = Attributes::makePredefinedString
                    ("TYPE", "Specifies the type of trim coil.",
                     {"PSI-BFIELD", "PSI-PHASE", "PSI-BFIELD-MIRRORED"});

    itsAttr[COEFNUM] = Attributes::makeRealArray
                       ("COEFNUM", "Radial profile: list of polynomial coefficients for the numerator (not for PSI-BFIELD-MIRRORED)");

    itsAttr[COEFDENOM] = Attributes::makeRealArray
                         ("COEFDENOM", "Radial profile: list of polynomial coefficients for the denominator (not for PSI-BFIELD-MIRRORED)");

    itsAttr[COEFNUMPHI] = Attributes::makeRealArray
                          ("COEFNUMPHI", "Angular profile: list of polynomial coefficients for the numerator (not for PSI-BFIELD-MIRRORED)");

    itsAttr[COEFDENOMPHI] = Attributes::makeRealArray
                            ("COEFDENOMPHI", "Angular profile: list of polynomial coefficients for the denominator (not for PSI-BFIELD-MIRRORED)");

    itsAttr[BMAX] = Attributes::makeReal
                    ("BMAX", "Maximum magnetic field [T]");

    itsAttr[PHIMIN] = Attributes::makeReal
                      ("PHIMIN", "Minimal azimuth [deg] (default 0)", 0.0);

    itsAttr[PHIMAX] = Attributes::makeReal
                      ("PHIMAX", "Maximal azimuth [deg] (default 360)", 360.0);

    itsAttr[RMIN] = Attributes::makeReal
                    ("RMIN", "Minimum radius [mm]");

    itsAttr[RMAX] = Attributes::makeReal
                    ("RMAX", "Maximum radius [mm]");

    itsAttr[SLPTC] = Attributes::makeReal
                     ("SLPTC", "Slopes of the rising edge [1/mm] (for PSI-BFIELD-MIRRORED)");

    registerOwnership(AttributeHandler::STATEMENT);

    OpalTrimCoil* defTrimCoil = clone("UNNAMED_TRIMCOIL");
    defTrimCoil->builtin = true;
    try {
        defTrimCoil->update();
        OpalData::getInstance()->define(defTrimCoil);
    } catch(...) {
        delete defTrimCoil;
    }
}


OpalTrimCoil::OpalTrimCoil(const std::string& name, OpalTrimCoil* parent):
    Definition(name, parent),
    trimcoil_m(nullptr)
{}


OpalTrimCoil::~OpalTrimCoil() {
}


bool OpalTrimCoil::canReplaceBy(Object* object) {
    // Can replace only by another trim coil.
    return dynamic_cast<OpalTrimCoil*>(object) != nullptr;
}


OpalTrimCoil* OpalTrimCoil::clone(const std::string& name) {
    return new OpalTrimCoil(name, this);
}


void OpalTrimCoil::execute() {
    update();
}


OpalTrimCoil* OpalTrimCoil::find(const std::string& name) {
    OpalTrimCoil* trimcoil = dynamic_cast<OpalTrimCoil*>(OpalData::getInstance()->find(name));

    if (trimcoil == nullptr) {
        throw OpalException("OpalTrimCoil::find()", "OpalTrimCoil \"" + name + "\" not found.");
    }
    return trimcoil;
}


void OpalTrimCoil::update() {
    // Set default name.
    if (getOpalName().empty()) setOpalName("UNNAMED_TRIMCOIL");
}


OpalTrimCoil::TrimCoilType OpalTrimCoil::getTrimCoilType() const {
    static const std::map<std::string, TrimCoilType> stringTrimCoilType_s = {
        {"PSI-BFIELD",          TrimCoilType::BFIELD},
        {"PSI-PHASE",           TrimCoilType::PHASE},
        {"PSI-BFIELD-MIRRORED", TrimCoilType::BFIELDMIRRORED},
    };

    if (!itsAttr[TYPE]) {
        throw OpalException("OpalTrimCoil::getTrimCoilType",
                            "The attribute \"TYPE\" isn't set for the \"TRIMCOIL\" statement");
    } else {
        return stringTrimCoilType_s.at(Attributes::getString(itsAttr[TYPE]));
    }
}


void OpalTrimCoil::initOpalTrimCoil() {
    if (trimcoil_m != nullptr) return;

    TrimCoilType type = getTrimCoilType();

    double bmax   = Attributes::getReal(itsAttr[BMAX]);
    double phimin = Attributes::getReal(itsAttr[PHIMIN]);
    double phimax = Attributes::getReal(itsAttr[PHIMAX]);
    double rmin   = Attributes::getReal(itsAttr[RMIN]);
    double rmax   = Attributes::getReal(itsAttr[RMAX]);

    if (type == TrimCoilType::BFIELD || type == TrimCoilType::PHASE) {
        std::vector<double> coefnum      = Attributes::getRealArray(itsAttr[COEFNUM]);
        std::vector<double> coefdenom    = Attributes::getRealArray(itsAttr[COEFDENOM]);
        std::vector<double> coefnumphi   = Attributes::getRealArray(itsAttr[COEFNUMPHI]);
        std::vector<double> coefdenomphi = Attributes::getRealArray(itsAttr[COEFDENOMPHI]);
        if (type == TrimCoilType::BFIELD) {
            trimcoil_m = std::unique_ptr<TrimCoilBFit>     (new TrimCoilBFit    (bmax, rmin, rmax, coefnum, coefdenom, coefnumphi, coefdenomphi));
        } else { // TrimCoilType::PHASE
            trimcoil_m = std::unique_ptr<TrimCoilPhaseFit> (new TrimCoilPhaseFit(bmax, rmin, rmax, coefnum, coefdenom, coefnumphi, coefdenomphi));
        }
    } else if (type == TrimCoilType::BFIELDMIRRORED) {
        double slope = Attributes::getReal(itsAttr[SLPTC]);
        trimcoil_m = std::unique_ptr<TrimCoilMirrored>     (new TrimCoilMirrored(bmax, rmin, rmax, slope));
    }

    trimcoil_m->setAzimuth(phimin, phimax);

    *gmsg << level3 << *this << endl;
}


Inform& OpalTrimCoil::print(Inform& os) const {
    os << "* ******************************** T R I M C O I L ********************************\n"
       << "* TRIMCOIL       " << getOpalName() << '\n'
       << "* TYPE           " << Attributes::getString(itsAttr[TYPE]) << '\n';

    TrimCoilType type = getTrimCoilType();
    if (type == TrimCoilType::BFIELD || type == TrimCoilType::PHASE) {
        printPolynom(os,itsAttr[COEFNUM]);
        printPolynom(os,itsAttr[COEFDENOM]);
        printPolynom(os,itsAttr[COEFNUMPHI]);
        printPolynom(os,itsAttr[COEFDENOMPHI]);
    }
    os << "* BMAX           " << Attributes::getReal(itsAttr[BMAX])   << '\n'
       << "* RMIN           " << Attributes::getReal(itsAttr[RMIN])   << '\n'
       << "* RMAX           " << Attributes::getReal(itsAttr[RMAX])   << '\n'
       << "* PHIMIN         " << Attributes::getReal(itsAttr[PHIMIN]) << '\n'
       << "* PHIMAX         " << Attributes::getReal(itsAttr[PHIMAX]) << '\n';
    if (type == TrimCoilType::BFIELDMIRRORED) {
        os << "* SLPTC          " << Attributes::getReal(itsAttr[SLPTC]) << '\n';
    }
    os << "* *********************************************************************************" << endl;
    return os;
}

void OpalTrimCoil::printPolynom(Inform& os, const Attribute& attr) const {
    std::stringstream ss;
    std::vector<double> coef = Attributes::getRealArray(attr);
    for (std::size_t i = 0; i < coef.size(); ++i) {
        ss << ((i > 0) ? "+ " : "") << coef[i]
           << ((i > 0) ? (" * x^" + std::to_string(i)) : "") << ' ';
    }
    os << "* POLYNOM " << attr.getName() << "   " << ss.str() << '\n';
}