Commit 071be009 authored by snuverink_j's avatar snuverink_j
Browse files

Merge branch 'master' of gitlab.psi.ch:OPAL/src

parents f8eb4330 f8636c1e
......@@ -87,21 +87,21 @@ Cyclotron::~Cyclotron() {
}
void Cyclotron::applyTrimCoil_m(const double r, const double z, double *br, double *bz) {
void Cyclotron::applyTrimCoil_m(const double r, const double z, const double tet_rad, double *br, double *bz) {
for (auto trimcoil : trimcoils_m) {
trimcoil->applyField(r,z,br,bz);
trimcoil->applyField(r,tet_rad,z,br,bz);
}
}
void Cyclotron::applyTrimCoil(const double r, const double z, double& br, double& bz) {
void Cyclotron::applyTrimCoil(const double r, const double z, const double tet_rad, double& br, double& bz) {
//Changed from > to >= to include case where bz == 0 and trimCoilThreshold_m == 0 -DW
if (std::abs(bz) >= trimCoilThreshold_m) {
applyTrimCoil_m(r, z, &br, &bz);
applyTrimCoil_m(r, z, tet_rad, &br, &bz);
}
else {
// make sure to have a smooth transition
double tmp_bz = 0.0;
applyTrimCoil_m(r, z, &br, &tmp_bz);
applyTrimCoil_m(r, z, tet_rad, &br, &tmp_bz);
bz += tmp_bz * std::abs(bz) / trimCoilThreshold_m;
}
}
......@@ -411,7 +411,7 @@ bool Cyclotron::apply(const Vector_t &R, const Vector_t &P, const double &t, Vec
/* Bz */
double bz = - bzint;
this->applyTrimCoil(rad, R[2], br, bz);
this->applyTrimCoil(rad, R[2], tet_rad, br, bz);
/* Br Btheta -> Bx By */
B[0] = br * cos(tet_rad) - bt * sin(tet_rad);
......@@ -525,7 +525,7 @@ void Cyclotron::apply(const double& rad, const double& z,
const double& tet_rad, double& br,
double& bt, double& bz) {
this->interpolate(rad, tet_rad, br, bt, bz);
this->applyTrimCoil(rad, z, br, bz);
this->applyTrimCoil(rad, z, tet_rad, br, bz);
}
void Cyclotron::finalise() {
......
......@@ -97,8 +97,6 @@ public:
Cyclotron();
Cyclotron(const Cyclotron &);
void applyTrimCoil(const double r, const double z, double& br, double& bz);
virtual ~Cyclotron();
/// Apply visitor to Cyclotron.
......@@ -215,7 +213,10 @@ public:
void read(const int &fieldflag, const double &scaleFactor);
private:
void applyTrimCoil_m(const double r, const double z, double *br, double *bz);
/// Apply trim coils (calculate field contributions) with smooth field transition
void applyTrimCoil (const double r, const double z, const double tet_rad, double& br, double& bz);
/// Apply trim coils (calculate field contributions)
void applyTrimCoil_m(const double r, const double z, const double tet_rad, double *br, double *bz);
protected:
......
set (_SRCS
OpalTrimCoil.cpp
TrimCoil.cpp
TrimCoilFit.cpp
TrimCoilBFit.cpp
TrimCoilPhaseFit.cpp
TrimCoilMirrored.cpp
......@@ -15,6 +16,7 @@ ADD_OPAL_SOURCES(${_SRCS})
set (HDRS
OpalTrimCoil.h
TrimCoil.h
TrimCoilFit.h
TrimCoilBFit.h
TrimCoilPhaseFit.h
TrimCoilMirrored.h
......
......@@ -21,7 +21,11 @@ namespace {
TYPE, // The type of trim coil
COEFNUM, //
COEFDENOM, //
COEFNUMPHI,
COEFDENOMPHI,
BMAX, //
PHIMIN,
PHIMAX,
RMIN, //
RMAX, //
SLPTC,
......@@ -37,23 +41,35 @@ OpalTrimCoil::OpalTrimCoil():
("TYPE", "Specifies the type of trim coil: PSI-BFIELD, PSI-PHASE, PSI-BFIELD-MIRRORED");
itsAttr[COEFNUM] = Attributes::makeRealArray
("COEFNUM", "List of polynomial coefficients for the numerator");
("COEFNUM", "Radial profile: list of polynomial coefficients for the numerator (not for PSI-BFIELD-MIRRORED)");
itsAttr[COEFDENOM] = Attributes::makeRealArray
("COEFDENOM", "List of polynomial coefficients for the denominator");
("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 in Tesla.");
itsAttr[PHIMIN] = Attributes::makeReal
("PHIMIN", "Minimal azimuth [deg] (default 0)");
itsAttr[PHIMAX] = Attributes::makeReal
("PHIMAX", "Maximal azimuth [deg] (default 360)");
itsAttr[RMIN] = Attributes::makeReal
("RMIN", "Minimum radius in millimeters.");
("RMIN", "Minimum radius [mm].");
itsAttr[RMAX] = Attributes::makeReal
("RMAX", "Maximum radius in millimeters.");
("RMAX", "Maximum radius [mm].");
itsAttr[SLPTC] = Attributes::makeReal
("SLPTC", "Slopes of the rising edge [1/mm] (for PSI-BFIELD-MIRRORED)");
registerOwnership(AttributeHandler::STATEMENT);
......@@ -76,7 +92,6 @@ OpalTrimCoil::OpalTrimCoil(const std::string &name, OpalTrimCoil *parent):
OpalTrimCoil::~OpalTrimCoil() {
}
......@@ -114,20 +129,24 @@ void OpalTrimCoil::update() {
void OpalTrimCoil::initOpalTrimCoil() {
if (trimcoil_m != nullptr) return;
std::string type = Util::toUpper(Attributes::getString(itsAttr[TYPE]));
double bmax = Attributes::getReal(itsAttr[BMAX]);
double rmin = Attributes::getReal(itsAttr[RMIN]);
double rmax = Attributes::getReal(itsAttr[RMAX]);
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 == "PSI-BFIELD" || type == "PSI-PHASE") {
std::vector<double> coefnum = Attributes::getRealArray(itsAttr[COEFNUM]);
std::vector<double> coefdenom = Attributes::getRealArray(itsAttr[COEFDENOM]);
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 == "PSI-BFIELD")
trimcoil_m = std::unique_ptr<TrimCoilBFit> (new TrimCoilBFit (bmax, rmin, rmax, coefnum, coefdenom));
trimcoil_m = std::unique_ptr<TrimCoilBFit> (new TrimCoilBFit (bmax, rmin, rmax, coefnum, coefdenom, coefnumphi, coefdenomphi));
else // type == "PSI-PHASE"
trimcoil_m = std::unique_ptr<TrimCoilPhaseFit> (new TrimCoilPhaseFit(bmax, rmin, rmax, coefnum, coefdenom));
trimcoil_m = std::unique_ptr<TrimCoilPhaseFit> (new TrimCoilPhaseFit(bmax, rmin, rmax, coefnum, coefdenom, coefnumphi, coefdenomphi));
} else if (type == "PSI-BFIELD-MIRRORED") {
double slope = Attributes::getReal(itsAttr[SLPTC]);
......@@ -136,7 +155,9 @@ void OpalTrimCoil::initOpalTrimCoil() {
throw OpalException("OpalTrimCoil::initOpalTrimCoil",
type + " is not a valid trim coil type");
}
trimcoil_m->setAzimuth(phimin, phimax);
*gmsg << level3 << *this << endl;
}
......@@ -144,33 +165,32 @@ 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';
std::string type = Util::toUpper(Attributes::getString(itsAttr[TYPE]));
if (type == "PSI-BFIELD" || type == "PSI-PHASE") {
std::vector<double> coefnum = Attributes::getRealArray(itsAttr[COEFNUM]);
std::stringstream ss;
for (std::size_t i = 0; i < coefnum.size(); ++i) {
ss << ((i > 0) ? "+ " : "") << coefnum[i]
<< ((i > 0) ? (" * x^" + std::to_string(i)) : "") << ' ';
}
os << "* POLYNOM NUM " << ss.str() << '\n';
std::vector<double> coefdenom = Attributes::getRealArray(itsAttr[COEFDENOM]);
ss.str("");
for (std::size_t i = 0; i < coefdenom.size(); ++i) {
ss << ((i > 0) ? "+ " : "") << coefdenom[i]
<< ((i > 0) ? (" * x^" + std::to_string(i)) : "") << ' ';
}
os << "* POLYNOM DENOM " << ss.str() << '\n';
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';
if (Util::toUpper(Attributes::getString(itsAttr[TYPE])) == "PSI-BFIELD-MIRRORED") {
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';
}
\ No newline at end of file
......@@ -5,6 +5,7 @@
#include <memory>
#include "AbstractObjects/Definition.h"
class Attribute;
class TrimCoil;
// Class OpalTrimCoil
......@@ -56,6 +57,8 @@ private:
/// Private copy constructor, called by clone
OpalTrimCoil(const std::string &name, OpalTrimCoil *parent);
/// Helper method for printing
void printPolynom(Inform& os, const Attribute& attr) const;
};
inline Inform &operator<<(Inform &os, const OpalTrimCoil &b) {
......
#include "TrimCoil.h"
#include <cmath>
#include "Physics/Physics.h"
TrimCoil::TrimCoil(double bmax,
double rmin,
double rmax)
......@@ -12,7 +16,21 @@ TrimCoil::TrimCoil(double bmax,
bmax_m = bmax * 10.0;
}
void TrimCoil::applyField(const double r, const double z, double *br, double *bz)
void TrimCoil::applyField(const double r, const double z, const double phi_rad, double *br, double *bz)
{
if (std::abs(bmax_m) < 1e-20) return;
if ((phimin_m <= phimax_m && (phi_rad < phimin_m || phi_rad > phimax_m)) ||
(phimin_m > phimax_m && (phi_rad < phimin_m && phi_rad > phimax_m)) ) return;
doApplyField(r,z,phi_rad,br,bz);
}
void TrimCoil::setAzimuth(const double phimin, const double phimax)
{
doApplyField(r,z,br,bz);
}
\ No newline at end of file
// phi convert to rad in [0,two pi]
if (phimin_m < 0) phimin_m += 360;
if (phimax_m < 0) phimax_m += 360;
phimin_m = phimin * Physics::deg2rad;
phimax_m = phimax * Physics::deg2rad;
}
#ifndef TRIM_COIL_H
#define TRIM_COIL_H
#include "Physics/Physics.h"
/// Abstract TrimCoil class
class TrimCoil {
......@@ -10,8 +12,9 @@ public:
TrimCoil(double bmax, double rmin, double rmax);
/// Apply the trim coil at position r and z to Bfields br and bz
/// Calls virtual doApplyField
void applyField(const double r, const double z, double *br, double *bz);
void applyField(const double r, const double z, const double phi_rad, double *br, double *bz);
/// Set azimuthal range
void setAzimuth(const double phimin, const double phimax);
virtual ~TrimCoil() { };
protected:
......@@ -22,11 +25,15 @@ protected:
double rmin_m;
/// Maximum radius (m)
double rmax_m;
/// Minimal azimuth (rad)
double phimin_m = 0.0;
/// Maximal azimuth (rad)
double phimax_m = Physics::two_pi;
private:
/// virtual implementation of applyField
virtual void doApplyField(const double r, const double z, double *br, double *bz) = 0;
virtual void doApplyField(const double r, const double z, const double phi_rad, double *br, double *bz) = 0;
};
#endif //TRIM_COIL_H
#include "TrimCoils/TrimCoilBFit.h"
#include <cmath>
TrimCoilBFit::TrimCoilBFit(double bmax,
double rmin,
double rmax,
const std::vector<double>& coefnum,
const std::vector<double>& coefdenom):
TrimCoil(bmax, rmin, rmax),
coefnum_m(coefnum),
coefdenom_m(coefdenom)
{
// normal polynom if no denominator coefficients (denominator = 1)
if (coefdenom_m.empty())
coefdenom_m.push_back(1.0);
}
const std::vector<double>& coefdenom,
const std::vector<double>& coefnumphi,
const std::vector<double>& coefdenomphi):
TrimCoilFit(bmax, rmin, rmax, coefnum, coefdenom, coefnumphi, coefdenomphi)
{}
void TrimCoilBFit::doApplyField(const double r, const double z, double *br, double *bz)
void TrimCoilBFit::doApplyField(const double r, const double z, const double phi_rad, double *br, double *bz)
{
if (std::abs(bmax_m) < 1e-20) return;
// check range
if (r < rmin_m || r > rmax_m) return;
double num = 0.0; // numerator
double dnum_dr = 0.0; // derivative of numerator
double powr = 1.0; // power of r
// add constant
num += coefnum_m[0];
for (std::size_t i = 1; i < coefnum_m.size(); ++i) {
dnum_dr += coefnum_m[i] * powr * i;
powr *= r;
num += coefnum_m[i] * powr;
}
double denom = 0.0; // denominator
double ddenom_dr = 0.0; // derivative of denominator
powr = 1.0; // power of r
// add constant
denom += coefdenom_m[0];
for (std::size_t i = 1; i < coefdenom_m.size(); ++i) {
ddenom_dr += coefdenom_m[i] * powr * i;
powr *= r;
denom += coefdenom_m[i] * powr;
}
double btr = num / denom;
// derivative of dr with quotient rule
double dr = (dnum_dr * denom - ddenom_dr * num) / (denom*denom);
double btr = 0.0, dr = 0.0;
calculateRationalFunction(RADIUS, r, btr, dr);
double phi = 0.0, dphi = 0.0;
calculateRationalFunction(PHI, phi_rad, phi, dphi);
//std::cout << "r " << r << " dr " << dr << std::endl;
*bz += bmax_m * btr;
*br += bmax_m * dr * z;
*bz += bmax_m * btr * phi;
*br += bmax_m * (dr * phi + btr*dphi) * z;
}
\ No newline at end of file
#ifndef TRIM_COILBFIT_H
#define TRIM_COILBFIT_H
#include "TrimCoils/TrimCoil.h"
#include "TrimCoils/TrimCoilFit.h"
#include <vector>
......@@ -9,25 +9,24 @@
/// General rational function fit
/// https://gitlab.psi.ch/OPAL/src/issues/157
class TrimCoilBFit : public TrimCoil {
class TrimCoilBFit : public TrimCoilFit {
public:
TrimCoilBFit(double bmax,
double rmin,
double rmax,
const std::vector<double>& coefnum,
const std::vector<double>& coefdenom);
const std::vector<double>& coefdenom,
const std::vector<double>& coefnumphi,
const std::vector<double>& coefdenomphi);
virtual ~TrimCoilBFit() { };
private:
TrimCoilBFit() = delete;
std::vector<double> coefnum_m;
std::vector<double> coefdenom_m;
/// @copydoc TrimCoil::doApplyField
virtual void doApplyField(const double r, const double z, double *br, double *bz);
virtual void doApplyField(const double r, const double z, const double phi_rad, double *br, double *bz);
};
#endif //TRIM_COILBFIT_H
#include "TrimCoils/TrimCoilFit.h"
#include <cmath>
TrimCoilFit::TrimCoilFit(double bmax,
double rmin,
double rmax,
const std::vector<double>& coefnum,
const std::vector<double>& coefdenom,
const std::vector<double>& coefnumphi,
const std::vector<double>& coefdenomphi):
TrimCoil(bmax, rmin, rmax)
{
coefs.resize(4);
coefs[NUM] = coefnum;
coefs[DENOM] = coefdenom;
coefs[NUMPHI] = coefnumphi;
coefs[DENOMPHI] = coefdenomphi;
// normal polynom if no denominator coefficients (denominator = 1)
if (coefs[DENOM].empty())
coefs[DENOM].push_back(1.0);
if (coefs[DENOMPHI].empty())
coefs[DENOMPHI].push_back(1.0);
// default constant nominator
if (coefs[NUM].empty())
coefs[NUM].push_back(1.0);
if (coefs[NUMPHI].empty())
coefs[NUMPHI].push_back(1.0);
}
void TrimCoilFit::calculateRationalFunction(FunctionType type, double val, double& quot, double& der_quot) const
{
double num = 0.0; // numerator
double dnum = 0.0; // derivative of numerator
double powval = 1.0; // power of value
const std::vector<double>& coefnum = coefs[type];
const std::vector<double>& coefdenom = coefs[type+1];
// add constant
num += coefnum[0];
for (std::size_t i = 1; i < coefnum.size(); ++i) {
dnum += coefnum[i] * powval * i;
powval *= val;
num += coefnum[i] * powval;
}
double denom = 0.0; // denominator
double ddenom = 0.0; // derivative of denominator
powval = 1.0; // power of value
// add constant
denom += coefdenom[0];
for (std::size_t i = 1; i < coefdenom.size(); ++i) {
ddenom += coefdenom[i] * powval * i;
powval *= val;
denom += coefdenom[i] * powval;
}
quot = num / denom;
// derivative with quotient rule
der_quot = (dnum * denom - ddenom * num) / (denom*denom);
}
void TrimCoilFit::calculateRationalFunction(FunctionType type, double val, double& quot, double& der_quot, double& der2_quot) const
{
double num = 0.0; // numerator
double d_num = 0.0; // derivative of numerator
double d2_num = 0.0; // second derivative of numerator
double powval = 1.0; // power of r
const std::vector<double>& coefnum = coefs[type];
const std::vector<double>& coefdenom = coefs[type+1];
unsigned int order = coefnum.size();
// add constant and first term
num += coefnum[0];
if (order > 1) {
num += coefnum[1] * val;
d_num += coefnum[1];
}
for (std::size_t i = 2; i < coefnum.size(); ++i) {
d2_num += coefnum[i] * powval * i * (i-1);
powval *= val; // r^(i-1)
d_num += coefnum[i] * powval * i;
num += coefnum[i] * powval * val;
}
double denom = 0.0; // denominator
double d_denom = 0.0; // derivative of denominator
double d2_denom = 0.0; // derivative of denominator
powval = 1.0; // power of r
order = coefdenom.size();
// add constant
denom += coefdenom[0];
if (order > 1) {
denom += coefdenom[1] * val;
d_denom += coefdenom[1];
}
for (std::size_t i = 2; i < coefdenom.size(); ++i) {
d2_denom += coefdenom[i] * powval * i * (i-1);
powval *= val;
d_denom += coefdenom[i] * powval * i;
denom += coefdenom[i] * powval * val;
}
quot = num / denom;
// derivative of phase with quotient rule (B - field)
der_quot = (d_num * denom - d_denom * num) / (denom*denom);
// second derivitive of phase (dB/dr) with quotient rule
// (d2_num - 2*(num/denom)' * d_denom - (num/denom) * d2_denom) / denom
der2_quot = (d2_num - 2*der_quot*d_denom - quot * d2_denom) / denom;
}
\ No newline at end of file
#ifndef TRIM_COILFIT_H
#define TRIM_COILFIT_H
#include "TrimCoils/TrimCoil.h"
#include <vector>
/// Abstract TrimCoilFit class
/// General rational function fit
class TrimCoilFit : public TrimCoil {
public:
TrimCoilFit(double bmax,
double rmin,
double rmax,
const std::vector<double>& coefnum,
const std::vector<double>& coefdenom,
const std::vector<double>& coefnumphi,
const std::vector<double>& coefdenomphi);
virtual ~TrimCoilFit() {};
protected:
enum PolynomType {NUM, DENOM, NUMPHI, DENOMPHI};
enum FunctionType {RADIUS=0, PHI=2};
/// calculate rational function and its first derivative
void calculateRationalFunction(FunctionType, double value, double& quot, double& der_quot) const;
/// calculate rational function and its first and second derivative
void calculateRationalFunction(FunctionType, double value, double& quot, double& der_quot, double& der2_quot) const;
private:
TrimCoilFit() = delete;
<