OpalCyclotron.cpp

Go to the documentation of this file.

//
// Class OpalCyclotron
//   The OpalCyclotron element.
//
// Copyright (c) 200x - 2020, Paul Scherrer Institut, 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 Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "Elements/OpalCyclotron.h"

#include "AbstractObjects/Attribute.h"
#include "Attributes/Attributes.h"
#include "BeamlineCore/CyclotronRep.h"
#include "Physics/Units.h"
#include "Structure/BoundaryGeometry.h"
#include "TrimCoils/OpalTrimCoil.h"
#include "TrimCoils/TrimCoil.h"
#include "Utilities/OpalException.h"

#include <numeric>


OpalCyclotron::OpalCyclotron():
    OpalElement(SIZE, "CYCLOTRON",
                "The \"CYCLOTRON\" defines a cyclotron"),
    obgeo_m(nullptr)  {
    itsAttr[CYHARMON] = Attributes::makeReal
        ("CYHARMON", "The harmonic number of the cyclotron");

    itsAttr[SYMMETRY] = Attributes::makeReal
        ("SYMMETRY", "Defines how the field is stored");

    itsAttr[RINIT] = Attributes::makeReal
        ("RINIT", "Initial radius of the reference particle [mm]");

    itsAttr[PRINIT] = Attributes::makeReal
        ("PRINIT", "Initial radial momentum of the reference particle, pr=beta_r*gamma");

    itsAttr[PHIINIT] = Attributes::makeReal
        ("PHIINIT", "Initial azimuth of the reference particle [deg]");

    itsAttr[ZINIT] = Attributes::makeReal
        ("ZINIT", "Initial z-coordinate of the reference particle [mm] (default=0 mm)", 0.0);

    itsAttr[PZINIT] = Attributes::makeReal
        ("PZINIT", "Initial vertical momentum of the reference particle, pz=beta_z*gamma (default=0)", 0.0);

    itsAttr[FMAPFN] = Attributes::makeString
        ("FMAPFN", "Filename for the B fieldmap");

    itsAttr[BSCALE] = Attributes::makeReal
         ("BSCALE", "Scale factor for the B-field (default=1)", 1.0);

    itsAttr[RFFREQ] = Attributes::makeRealArray
        ("RFFREQ", "RF Frequency(ies) [MHz]");

    itsAttr[ESCALE] = Attributes::makeRealArray
        ("ESCALE", "Scale factor for the RF field(s)");

    itsAttr[SUPERPOSE] = Attributes::makeBoolArray
        ("SUPERPOSE", "If TRUE, all of the electric field maps are superposed, only used when TYPE=BANDRF");

    itsAttr[RFMAPFN] = Attributes::makeStringArray
        ("RFMAPFN", "Filename(s) for the RF fieldmap(s)");

    itsAttr[RFFCFN] = Attributes::makeStringArray
        ("RFFCFN", "Filename(s) for the RF Frequency Coefficients");

    itsAttr[RFVCFN] = Attributes::makeStringArray
        ("RFVCFN", "Filename(s) for the RF Voltage Coefficients");

    itsAttr[RFPHI] = Attributes::makeRealArray
        ("RFPHI", "Initial phase(s) of the electric field map(s) [rad]");

    itsAttr[MINZ] = Attributes::makeReal
        ("MINZ","Minimal vertical extent of the machine [mm] (default=-10000 mm)",-10000.0);

    itsAttr[MAXZ] = Attributes::makeReal
        ("MAXZ","Maximal vertical extent of the machine [mm] (default=10000 mm)",10000.0);

    itsAttr[MINR] = Attributes::makeReal
        ("MINR","Minimal radial extent of the machine [mm] (default=0 mm)", 0.0);

    itsAttr[MAXR] = Attributes::makeReal
        ("MAXR","Maximal radial extent of the machine [mm] (default=10000 mm)", 10000.0);

    itsAttr[GEOMETRY] = Attributes::makeString
        ("GEOMETRY", "Boundary Geometry for the cyclotron");

    itsAttr[FMLOWE] = Attributes::makeReal
        ("FMLOWE", "Minimal Energy [GeV] the fieldmap can accept. Used in GAUSSMATCHED", -1.0);

    itsAttr[FMHIGHE] = Attributes::makeReal
        ("FMHIGHE", "Maximal Energy [GeV] the fieldmap can accept. Used in GAUSSMATCHED", -1.0);

    itsAttr[SPIRAL] = Attributes::makeBool
        ("SPIRAL", "Flag whether or not this is a spiral inflector simulation (default=FALSE)", false);

    itsAttr[TRIMCOILTHRESHOLD] = Attributes::makeReal
        ("TRIMCOILTHRESHOLD", "Minimum magnetic field [T] for which trim coils are applied (default=0.0)", 0.0);

    itsAttr[TRIMCOIL] = Attributes::makeStringArray
        ("TRIMCOIL", "List of trim coils");

    registerOwnership();

    setElement(new CyclotronRep("CYCLOTRON"));
}

OpalCyclotron::OpalCyclotron(const std::string& name, OpalCyclotron* parent):
    OpalElement(name, parent),
    obgeo_m(nullptr) {
    setElement(new CyclotronRep(name));
}


OpalCyclotron::~OpalCyclotron()
{}


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


void OpalCyclotron::update() {

    CyclotronRep* cycl =
        dynamic_cast<CyclotronRep*>(getElement());

    std::string fmapfm = Attributes::getString(itsAttr[FMAPFN]);
    std::string type   = Attributes::getString(itsAttr[TYPE]);

    double harmnum  = Attributes::getReal(itsAttr[CYHARMON]);
    double symmetry = Attributes::getReal(itsAttr[SYMMETRY]);
    double bscale   = Attributes::getReal(itsAttr[BSCALE]);

    double rinit    = Units::mm2m * Attributes::getReal(itsAttr[RINIT]);
    double phiinit  = Units::deg2rad * Attributes::getReal(itsAttr[PHIINIT]);
    double zinit    = Units::mm2m * Attributes::getReal(itsAttr[ZINIT]);
    double prinit   = Attributes::getReal(itsAttr[PRINIT]);
    double pzinit   = Attributes::getReal(itsAttr[PZINIT]);

    double minz = Units::mm2m * Attributes::getReal(itsAttr[MINZ]);
    double maxz = Units::mm2m * Attributes::getReal(itsAttr[MAXZ]);
    double minr = Units::mm2m * Attributes::getReal(itsAttr[MINR]);
    double maxr = Units::mm2m * Attributes::getReal(itsAttr[MAXR]);

    double fmLowE  = Units::GeV2MeV * Attributes::getReal(itsAttr[FMLOWE]);
    double fmHighE = Units::GeV2MeV * Attributes::getReal(itsAttr[FMHIGHE]);

    bool spiral_flag = Attributes::getBool(itsAttr[SPIRAL]);
    double trimCoilThreshold = Units::T2kG * Attributes::getReal(itsAttr[TRIMCOILTHRESHOLD]);

    cycl->setFieldMapFN(fmapfm);
    cycl->setSymmetry(symmetry);
    cycl->setBScale(bscale);

    cycl->setRinit(rinit);
    cycl->setPRinit(prinit);
    cycl->setPHIinit(phiinit);
    cycl->setZinit(zinit);
    cycl->setPZinit(pzinit);

    cycl->setCyclotronType(type);
    cycl->setCyclHarm(harmnum);

    cycl->setMinR(minr);
    cycl->setMaxR(maxr);
    cycl->setMinZ(minz);
    cycl->setMaxZ(maxz);

    cycl->setFMLowE(fmLowE);
    cycl->setFMHighE(fmHighE);

    cycl->setSpiralFlag(spiral_flag);
    cycl->setTrimCoilThreshold(trimCoilThreshold);

    cycl->setOutputFN(Attributes::getString(itsAttr[OUTFN]));

    std::vector<std::string> fm_str     = Attributes::getStringArray(itsAttr[RFMAPFN]);
    std::vector<std::string> rffcfn_str = Attributes::getStringArray(itsAttr[RFFCFN]);
    std::vector<std::string> rfvcfn_str = Attributes::getStringArray(itsAttr[RFVCFN]);
    std::vector<double> scale_str       = Attributes::getRealArray(itsAttr[ESCALE]);
    std::vector<double> phi_str         = Attributes::getRealArray(itsAttr[RFPHI]);
    std::vector<double> rff_str         = Attributes::getRealArray(itsAttr[RFFREQ]);
    std::vector<bool> superpose_str     = Attributes::getBoolArray(itsAttr[SUPERPOSE]);
    std::vector<std::string> trimcoil   = Attributes::getStringArray(itsAttr[TRIMCOIL]);

    if ( !trimcoil.empty() ) {

        std::vector<TrimCoil* > trimcoils;

        for (std::vector<std::string>::const_iterator tit = trimcoil.begin();
             tit != trimcoil.end(); ++tit)
        {
            OpalTrimCoil *tc = OpalTrimCoil::find(*tit);

            if ( tc ) {
                tc->initOpalTrimCoil();
                trimcoils.push_back(tc->trimcoil_m.get());
            }
        }
        cycl->setTrimCoils(trimcoils);
    }

    cycl->setRfPhi(phi_str);
    cycl->setEScale(scale_str);
    cycl->setRfFieldMapFN(fm_str);
    cycl->setRFFCoeffFN(rffcfn_str);
    cycl->setRFVCoeffFN(rfvcfn_str);
    cycl->setRfFrequ(rff_str);
    cycl->setSuperpose(superpose_str);

    if (itsAttr[GEOMETRY] && obgeo_m == nullptr) {
      obgeo_m = BoundaryGeometry::find(Attributes::getString(itsAttr[GEOMETRY]));
      if (obgeo_m) {
          cycl->setBoundaryGeometry(obgeo_m);
      }
    }

    // Transmit "unknown" attributes.
    OpalElement::updateUnknown(cycl);
}