diff --git a/src/Classic/AbsBeamline/Cyclotron.cpp b/src/Classic/AbsBeamline/Cyclotron.cpp
index 42b357e78cc4049f533944068244900bec99eeaa..344b194d3d01ed0a50e39c43aa955b08db4e0f06 100644
--- a/src/Classic/AbsBeamline/Cyclotron.cpp
+++ b/src/Classic/AbsBeamline/Cyclotron.cpp
@@ -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() {
diff --git a/src/Classic/AbsBeamline/Cyclotron.h b/src/Classic/AbsBeamline/Cyclotron.h
index d27dd5ec621e4d18fce25fdb0884d55d5ad5610d..7edebe04fef436995cf8b2f3b279b305fb703587 100644
--- a/src/Classic/AbsBeamline/Cyclotron.h
+++ b/src/Classic/AbsBeamline/Cyclotron.h
@@ -98,8 +98,6 @@ public:
Cyclotron();
Cyclotron(const Cyclotron &);
- void applyTrimCoil(const double r, const double z, double& br, double& bz);
-
virtual ~Cyclotron();
/// Apply visitor to Cyclotron.
@@ -216,7 +214,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:
diff --git a/src/Classic/TrimCoils/CMakeLists.txt b/src/Classic/TrimCoils/CMakeLists.txt
index fb649d834b9b3ac651babc09626023b8a2c57cd8..5baaf194c758c06bc1bec5bad643da9b4f42a0a2 100644
--- a/src/Classic/TrimCoils/CMakeLists.txt
+++ b/src/Classic/TrimCoils/CMakeLists.txt
@@ -1,6 +1,7 @@
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
diff --git a/src/Classic/TrimCoils/OpalTrimCoil.cpp b/src/Classic/TrimCoils/OpalTrimCoil.cpp
index edff5b06f2998c251ed234dd40ebb170b8f7c5e7..448ab41990f8f8bafede9d80bcafbdfd14c17751 100644
--- a/src/Classic/TrimCoils/OpalTrimCoil.cpp
+++ b/src/Classic/TrimCoils/OpalTrimCoil.cpp
@@ -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
diff --git a/src/Classic/TrimCoils/OpalTrimCoil.h b/src/Classic/TrimCoils/OpalTrimCoil.h
index 9fe782e3d8d623f39c6a9a355b8640c008131049..26728c1c13970a6aab6a8cb9cc1df318c08d2a85 100644
--- a/src/Classic/TrimCoils/OpalTrimCoil.h
+++ b/src/Classic/TrimCoils/OpalTrimCoil.h
@@ -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) {
diff --git a/src/Classic/TrimCoils/TrimCoil.cpp b/src/Classic/TrimCoils/TrimCoil.cpp
index 15b302a7370ed73bce21df52f7a3c2d62218798e..c34225e4d143107d7eeddc7bd990c71a9cc6d37a 100644
--- a/src/Classic/TrimCoils/TrimCoil.cpp
+++ b/src/Classic/TrimCoils/TrimCoil.cpp
@@ -1,5 +1,9 @@
#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;
+}
diff --git a/src/Classic/TrimCoils/TrimCoil.h b/src/Classic/TrimCoils/TrimCoil.h
index 35647ba5e22acedbed891437171c9c349d1abc89..254a868a4d9d493de3a16b3142c413da989cb7c8 100644
--- a/src/Classic/TrimCoils/TrimCoil.h
+++ b/src/Classic/TrimCoils/TrimCoil.h
@@ -1,6 +1,8 @@
#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
diff --git a/src/Classic/TrimCoils/TrimCoilBFit.cpp b/src/Classic/TrimCoils/TrimCoilBFit.cpp
index dcfd4d11eebe6afecfae1da74a20c25ef632b174..4129879b7a8dd503526e9896ea638899774d0ed2 100644
--- a/src/Classic/TrimCoils/TrimCoilBFit.cpp
+++ b/src/Classic/TrimCoils/TrimCoilBFit.cpp
@@ -1,58 +1,27 @@
#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
diff --git a/src/Classic/TrimCoils/TrimCoilBFit.h b/src/Classic/TrimCoils/TrimCoilBFit.h
index 128ea9815e241c1901138a8962c0210ccdbffd0a..b734cdadd83cc2d7a0512fcc3aa794ffdb435609 100644
--- a/src/Classic/TrimCoils/TrimCoilBFit.h
+++ b/src/Classic/TrimCoils/TrimCoilBFit.h
@@ -1,7 +1,7 @@
#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
diff --git a/src/Classic/TrimCoils/TrimCoilFit.cpp b/src/Classic/TrimCoils/TrimCoilFit.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6b671fef4f3f7e5b66c6b545b5893e9954a53b8f
--- /dev/null
+++ b/src/Classic/TrimCoils/TrimCoilFit.cpp
@@ -0,0 +1,118 @@
+#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
diff --git a/src/Classic/TrimCoils/TrimCoilFit.h b/src/Classic/TrimCoils/TrimCoilFit.h
new file mode 100644
index 0000000000000000000000000000000000000000..be23d9be868ca021b337fd3b8ab7157b3ab72ef9
--- /dev/null
+++ b/src/Classic/TrimCoils/TrimCoilFit.h
@@ -0,0 +1,40 @@
+#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;
+
+ /// rational function coefficients
+ std::vector<std::vector<double>> coefs;
+};
+
+#endif //TRIM_COILFIT_H
diff --git a/src/Classic/TrimCoils/TrimCoilMirrored.cpp b/src/Classic/TrimCoils/TrimCoilMirrored.cpp
index 2f9d3bdf59b75b7299d31e37f58f49a13612a809..3dea808760ded2e5f49455c59bd62147f0d3d78e 100644
--- a/src/Classic/TrimCoils/TrimCoilMirrored.cpp
+++ b/src/Classic/TrimCoils/TrimCoilMirrored.cpp
@@ -16,7 +16,7 @@ TrimCoilMirrored::TrimCoilMirrored(double bmax,
bslope_m = bslope / mm2m;
}
-void TrimCoilMirrored::doApplyField(const double r, const double z, double *br, double *bz)
+void TrimCoilMirrored::doApplyField(const double r, const double z, const double phi_rad, double *br, double *bz)
{
/// update bz and br with trim coil contributions
// for some discussion on the formulas see
@@ -24,8 +24,6 @@ void TrimCoilMirrored::doApplyField(const double r, const double z, double *br,
// https://gitlab.psi.ch/OPAL/src/issues/157
// https://gitlab.psi.ch/OPAL/src/issues/110
- if (std::abs(bmax_m) < 1e-20) return;
-
// unitless constants
const double Amax1 = 1;
const double Amax2 = 3;
diff --git a/src/Classic/TrimCoils/TrimCoilMirrored.h b/src/Classic/TrimCoils/TrimCoilMirrored.h
index 7787e9fd1cbbec7e9e5f0c424e023bda1a082e4d..6ce955f4a56fc9c8b993d422b8b839241b990765 100644
--- a/src/Classic/TrimCoils/TrimCoilMirrored.h
+++ b/src/Classic/TrimCoils/TrimCoilMirrored.h
@@ -22,7 +22,7 @@ private:
double bslope_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);
TrimCoilMirrored() = delete;
};
diff --git a/src/Classic/TrimCoils/TrimCoilPhaseFit.cpp b/src/Classic/TrimCoils/TrimCoilPhaseFit.cpp
index 693f0fbb8eccc711b8e3fa48d2eb4e85f88b260c..6796182823621f0530f46451fc4fc4aec518c09c 100644
--- a/src/Classic/TrimCoils/TrimCoilPhaseFit.cpp
+++ b/src/Classic/TrimCoils/TrimCoilPhaseFit.cpp
@@ -1,76 +1,29 @@
#include "TrimCoils/TrimCoilPhaseFit.h"
-#include <cmath>
-
TrimCoilPhaseFit::TrimCoilPhaseFit(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 TrimCoilPhaseFit::doApplyField(const double r, const double z, double *br, double *bz)
+void TrimCoilPhaseFit::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 d_num = 0.0; // derivative of numerator
- double d2_num = 0.0; // second derivative of numerator
- double powr = 1.0; // power of r
- unsigned int order = coefnum_m.size();
-
- // add constant and first term
- num += coefnum_m[0];
- if (order > 1) {
- num += coefnum_m[1] * r;
- d_num += coefnum_m[1];
- }
- for (std::size_t i = 2; i < coefnum_m.size(); ++i) {
- d2_num += coefnum_m[i] * powr * i * (i-1);
- powr *= r; // r^(i-1)
- d_num += coefnum_m[i] * powr * i;
- num += coefnum_m[i] * powr * r;
- }
-
- double denom = 0.0; // denominator
- double d_denom = 0.0; // derivative of denominator
- double d2_denom = 0.0; // derivative of denominator
- powr = 1.0; // power of r
- order = coefdenom_m.size();
-
- // add constant
- denom += coefdenom_m[0];
- if (order > 1) {
- denom += coefdenom_m[1] * r;
- d_denom += coefdenom_m[1];
- }
- for (std::size_t i = 2; i < coefdenom_m.size(); ++i) {
- d2_denom += coefdenom_m[i] * powr * i * (i-1);
- powr *= r;
- d_denom += coefdenom_m[i] * powr * i;
- denom += coefdenom_m[i] * powr * r;
- }
-
- double phase = num / denom;
-
- // derivative of phase with quotient rule (B - field)
- double d_phase = (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
- double d2_phase = (d2_num - 2*d_phase*d_denom - phase * d2_denom) / denom;
+ double phase=0.0, d_phase=0.0, d2_phase=0.0;
+ calculateRationalFunction(RADIUS, r, phase, d_phase, d2_phase);
+ double phi = 0.0, d_phi = 0.0, d2_phi=0.0;
+ calculateRationalFunction(PHI, phi_rad, phi, d_phi, d2_phi);
//std::cout << "r " << r << " dr " << dr << std::endl;
+ double derivative = d_phase * phi + phase * d_phi;
+ double der2 = d2_phase * phi + phase * d2_phi + 2*d_phi*d_phase;
- *bz += - bmax_m * d_phase;
- *br += - bmax_m * d2_phase * z;
+ *bz += - bmax_m * derivative;
+ *br += - bmax_m * der2 * z;
}
\ No newline at end of file
diff --git a/src/Classic/TrimCoils/TrimCoilPhaseFit.h b/src/Classic/TrimCoils/TrimCoilPhaseFit.h
index 22e17e7c4dcdd674e43843c5ba4e61b8583acfb2..dac12ed6aeed5064ee3001a97293cedd92f4900e 100644
--- a/src/Classic/TrimCoils/TrimCoilPhaseFit.h
+++ b/src/Classic/TrimCoils/TrimCoilPhaseFit.h
@@ -1,32 +1,31 @@
#ifndef TRIM_COILPHASEFIT_H
#define TRIM_COILPHASEFIT_H
-#include "TrimCoils/TrimCoil.h"
+#include "TrimCoils/TrimCoilFit.h"
#include <vector>
/// TrimCoilPhaseFit class
/// General rational function fit of the phase shift
-class TrimCoilPhaseFit : public TrimCoil {
+class TrimCoilPhaseFit : public TrimCoilFit {
public:
TrimCoilPhaseFit(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 ~TrimCoilPhaseFit() { };
private:
TrimCoilPhaseFit() = 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_COILPHASEFIT_H
diff --git a/tests/classic_src/AbsBeamline/TrimCoilTest.cpp b/tests/classic_src/AbsBeamline/TrimCoilTest.cpp
index 16b7cc6ed91cd92024fb271b01c2eccb3ff3538e..2030d212b5044f61d6222f6de4783a488212d8f8 100644
--- a/tests/classic_src/AbsBeamline/TrimCoilTest.cpp
+++ b/tests/classic_src/AbsBeamline/TrimCoilTest.cpp
@@ -6,6 +6,8 @@
#include "opal_test_utilities/SilenceTest.h"
+const double margin = 1e-7;
+
// Test TrimCoilBFit on zeros
TEST(TrimCoil, TrimCoilBFitZeros)
{
@@ -16,30 +18,37 @@ TEST(TrimCoil, TrimCoilBFitZeros)
double bmax = 0.0;
double rmin = 1000; // mm
double rmax = 2000;
+ double phi = 0.0; // rad
- TrimCoilBFit myTrimCoil(bmax, rmin, rmax, {}, {});
+ TrimCoilBFit myTrimCoil(bmax, rmin, rmax, {}, {}, {}, {});
const double one = 1.0;
double br = one, bz = one;
- myTrimCoil.applyField((rmin+rmax)*mm2m/2.0, 1, &br, &bz);
+ myTrimCoil.applyField((rmin+rmax)*mm2m/2.0, 1, 0, &br, &bz);
// not changed since bmax 0.0
- EXPECT_NEAR(br, one, 1e-7);
- EXPECT_NEAR(bz, one, 1e-7);
+ EXPECT_NEAR(br, one, margin);
+ EXPECT_NEAR(bz, one, margin);
bmax = 1.0;
- myTrimCoil = TrimCoilBFit(bmax, rmin, rmax, {}, {});
+ myTrimCoil = TrimCoilBFit(bmax, rmin, rmax, {}, {}, {}, {});
- myTrimCoil.applyField(rmin*mm2m - 1, 1, &br, &bz);
+ myTrimCoil.applyField(rmin*mm2m - 1, 1, phi, &br, &bz);
// not changed since r outside range
- EXPECT_NEAR(br, one, 1e-7);
- EXPECT_NEAR(bz, one, 1e-7);
+ EXPECT_NEAR(br, one, margin);
+ EXPECT_NEAR(bz, one, margin);
- myTrimCoil.applyField(rmax*mm2m + 1, 1, &br, &bz);
+ myTrimCoil.applyField(rmax*mm2m + 1, 1, phi, &br, &bz);
// not changed since r outside range
- EXPECT_NEAR(br, one, 1e-7);
- EXPECT_NEAR(bz, one, 1e-7);
+ EXPECT_NEAR(br, one, margin);
+ EXPECT_NEAR(bz, one, margin);
+
+ myTrimCoil.applyField(rmax*mm2m - 1, 1, phi, &br, &bz);
+ // default constant field
+ EXPECT_NEAR(br, one, margin);
+ EXPECT_NEAR(bz, 11.0, margin); // 1 + 10*1
}
+// TrimCoilBFit class
TEST(TrimCoil, TrimCoilBFit)
{
OpalTestUtilities::SilenceTest silencer;
@@ -47,26 +56,66 @@ TEST(TrimCoil, TrimCoilBFit)
double bmax = 1.0; // 1 T will be converted to 10 kG
double rmin = 0.0;
double rmax = 3000.0;
+ double phi = 0.0; // rad
// polynom 1 + 2*x + 3*x^2
- TrimCoilBFit myTrimCoil(bmax, rmin, rmax, {1.0, 2.0, 3.0}, {});
+ TrimCoilBFit myTrimCoil(bmax, rmin, rmax, {1.0, 2.0, 3.0}, {}, {}, {});
- double br = 1.0, bz = 1.0;
- myTrimCoil.applyField(2.0, 2.0, &br, &bz);
+ const double brStart = 1.0, bzStart = 1.0, zStart = 2.0;
+ double br = brStart, bz = bzStart;
+ myTrimCoil.applyField(2.0, zStart, phi, &br, &bz);
// 1 + 10*(1 + 2*2 + 3*2*2) = 171 ; 1 + 10*2*(2 + 6*2) = 281
- EXPECT_NEAR(bz, 171.0, 1e-7);
- EXPECT_NEAR(br, 281.0, 1e-7);
+ EXPECT_NEAR(bz, 171.0, margin);
+ EXPECT_NEAR(br, 281.0, margin);
// rational function (4 + 3*x) / (1 + 2*x)
- myTrimCoil = TrimCoilBFit(bmax, rmin, rmax, {4.0, 3.0}, {1.0, 2.0});
-
- br = 1.0, bz = 1.0;
- myTrimCoil.applyField(2.0, 2.0, &br, &bz);
- // 1 + 10*(4+3*2) / (1+2*2) = 21.0 ; 1 + 10*2*(-0.2) = -3
- EXPECT_NEAR(bz, 21.0, 1e-7);
- EXPECT_NEAR(br, -3.0, 1e-7);
+ myTrimCoil = TrimCoilBFit(bmax, rmin, rmax, {4.0, 3.0}, {1.0, 2.0}, {}, {});
+
+ br = brStart, bz = bzStart;
+ myTrimCoil.applyField(2.0, zStart, phi, &br, &bz);
+
+ double bzSolution = 21.0; // = 1 + 10*(4+3*2) / (1+2*2)
+ double brSolution = -3.0; // = 1 + 10*2*(-0.2)
+
+ EXPECT_NEAR(bz, bzSolution, margin);
+ EXPECT_NEAR(br, brSolution, margin);
+
+ // test phi angles
+ myTrimCoil.setAzimuth(10,180);
+ br = brStart, bz = bzStart;
+ myTrimCoil.applyField(2.0, zStart, 3.2, &br, &bz); // outside range
+ EXPECT_NEAR(bz, bzStart, margin);
+ myTrimCoil.applyField(2.0, zStart, 0.0, &br, &bz); // outside range
+ EXPECT_NEAR(bz, bzStart, margin);
+ myTrimCoil.applyField(2.0, zStart, 1.0, &br, &bz); // inside range
+ EXPECT_NEAR(bz, bzSolution, margin);
+
+ // test phi angles: first larger
+ myTrimCoil.setAzimuth(180,20);
+ br = brStart, bz = bzStart;
+ myTrimCoil.applyField(2.0, zStart, 1.0, &br, &bz); // outside range
+ EXPECT_NEAR(bz, bzStart, margin);
+ myTrimCoil.applyField(2.0, zStart, 3.2, &br, &bz); // inside range
+ EXPECT_NEAR(bz, bzSolution, margin);
+ myTrimCoil.applyField(2.0, zStart, 0.0, &br, &bz); // inside range
+ EXPECT_NEAR(bz, 2*bzSolution - bzStart, margin);
+
+ // Test Phi
+ // same rational function
+ myTrimCoil = TrimCoilBFit(bmax, rmin, rmax, {}, {}, {4.0, 3.0}, {1.0, 2.0});
+ br = brStart, bz = bzStart;
+ myTrimCoil.applyField(1.0, zStart, 2.0, &br, &bz);
+ EXPECT_NEAR(bz, bzSolution, margin);
+ EXPECT_NEAR(br, brSolution, margin);
+ // Both r and phi
+ myTrimCoil = TrimCoilBFit(bmax, rmin, rmax, {4.0, 3.0}, {1.0, 2.0}, {4.0, 3.0}, {1.0, 2.0});
+ br = brStart, bz = bzStart;
+ myTrimCoil.applyField(2.0, zStart, 2.0, &br, &bz);
+ EXPECT_NEAR(bz, 41.0, margin); // 1 + 10 * ((4+3*2) / (1+2*2)) **2
+ EXPECT_NEAR(br, -15.0, margin);
}
+// TrimCoilPhaseFit class
TEST(TrimCoil, TrimCoilPhaseFit)
{
OpalTestUtilities::SilenceTest silencer;
@@ -74,26 +123,42 @@ TEST(TrimCoil, TrimCoilPhaseFit)
double bmax = 1.0; // 1 T will be converted to 10 kG
double rmin = 0.0;
double rmax = 3000.0;
+ double phi = 0.0;
// polynom 1 + 2*x + 3*x^2
- TrimCoilPhaseFit myTrimCoil(bmax, rmin, rmax, {1.0, 2.0, 3.0}, {});
-
- double br = 1.0, bz = 1.0;
- myTrimCoil.applyField(2.0, 2.0, &br, &bz);
+ TrimCoilPhaseFit myTrimCoil(bmax, rmin, rmax, {1.0, 2.0, 3.0}, {}, {}, {});
+ const double brStart = 1.0, bzStart = 1.0, zStart = 2.0;
+ double br = brStart, bz = bzStart;
+ myTrimCoil.applyField(2.0, zStart, phi, &br, &bz);
// 1 - 10*(2 + 6*2) = -139; 1 - 2*10*6 = -119
- EXPECT_NEAR(bz, -139.0, 1e-7);
- EXPECT_NEAR(br, -119.0, 1e-7);
+ EXPECT_NEAR(bz, -139.0, margin);
+ EXPECT_NEAR(br, -119.0, margin);
// rational function (4 + 3*x) / (1 + 2*x)
- myTrimCoil = TrimCoilPhaseFit(bmax, rmin, rmax, {4.0, 3.0}, {1.0, 2.0});
-
- br = 1.0, bz = 1.0;
- myTrimCoil.applyField(2.0, 2.0, &br, &bz);
- // 1 - 10*(-0.2) = 3; 1 + 2*10*2*(-0.2*2) / 5 = -2.2
- EXPECT_NEAR(bz, 3.0, 1e-7);
- EXPECT_NEAR(br,-2.2, 1e-7);
+ myTrimCoil = TrimCoilPhaseFit(bmax, rmin, rmax, {4.0, 3.0}, {1.0, 2.0}, {}, {});
+
+ br = brStart, bz = bzStart;
+ myTrimCoil.applyField(2.0, zStart, phi, &br, &bz);
+ double bzSolution = 3.0; // = 1 - 10*(-0.2)
+ double brSolution = -2.2; // = 1 + 2*10*2*(-0.2*2) / 5
+ EXPECT_NEAR(bz, bzSolution, margin);
+ EXPECT_NEAR(br, brSolution, margin);
+ // Test Phi
+ // same rational function
+ myTrimCoil = TrimCoilPhaseFit(bmax, rmin, rmax, {}, {}, {4.0, 3.0}, {1.0, 2.0});
+ br = brStart, bz = bzStart;
+ myTrimCoil.applyField(1.0, zStart, 2.0, &br, &bz);
+ EXPECT_NEAR(bz, bzSolution, margin);
+ EXPECT_NEAR(br, brSolution, margin);
+ // Both r and phi
+ myTrimCoil = TrimCoilPhaseFit(bmax, rmin, rmax, {4.0, 3.0}, {1.0, 2.0}, {4.0, 3.0}, {1.0, 2.0});
+ br = brStart, bz = bzStart;
+ myTrimCoil.applyField(2.0, zStart, 2.0, &br, &bz);
+ EXPECT_NEAR(bz, 9.0, margin);
+ EXPECT_NEAR(br, -13.4, margin);
}
+// TrimCoilMirrored class
TEST(TrimCoil, TrimCoilMirrored)
{
OpalTestUtilities::SilenceTest silencer;
@@ -101,13 +166,15 @@ TEST(TrimCoil, TrimCoilMirrored)
double bmax = 1.0; // 1 T will be converted to 10 kG
double rmin = 0.0;
double rmax = 3000.0;
+ double phi = 0.0;
double bslope = 1./6.;
TrimCoilMirrored myTrimCoil(bmax, rmin, rmax, bslope);
double br = 1.0, bz = 1.0;
- myTrimCoil.applyField(2.0, 2.0, &br, &bz);
+ myTrimCoil.applyField(2.0, 2.0, phi, &br, &bz);
EXPECT_NEAR(bz,-6.1943868603626751, 1e-6);
EXPECT_NEAR(br, 1.0032755233321968, 1e-6);
-}
\ No newline at end of file
+}
+