From 9594f6f6e8105489a43ac33ff9e9beac356e5c9a Mon Sep 17 00:00:00 2001 From: Jochem Snuverink <jochem.snuverink@psi.ch> Date: Fri, 8 Feb 2019 14:58:30 +0100 Subject: [PATCH 1/4] refactor --- src/Classic/TrimCoils/TrimCoil.cpp | 4 ++++ src/Classic/TrimCoils/TrimCoilBFit.cpp | 1 - src/Classic/TrimCoils/TrimCoilMirrored.cpp | 2 -- src/Classic/TrimCoils/TrimCoilPhaseFit.cpp | 1 - 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Classic/TrimCoils/TrimCoil.cpp b/src/Classic/TrimCoils/TrimCoil.cpp index 15b302a73..51be9ebf9 100644 --- a/src/Classic/TrimCoils/TrimCoil.cpp +++ b/src/Classic/TrimCoils/TrimCoil.cpp @@ -1,5 +1,7 @@ #include "TrimCoil.h" +#include <cmath> + TrimCoil::TrimCoil(double bmax, double rmin, double rmax) @@ -14,5 +16,7 @@ TrimCoil::TrimCoil(double bmax, void TrimCoil::applyField(const double r, const double z, double *br, double *bz) { + if (std::abs(bmax_m) < 1e-20) return; + doApplyField(r,z,br,bz); } \ No newline at end of file diff --git a/src/Classic/TrimCoils/TrimCoilBFit.cpp b/src/Classic/TrimCoils/TrimCoilBFit.cpp index dcfd4d11e..197e2dbca 100644 --- a/src/Classic/TrimCoils/TrimCoilBFit.cpp +++ b/src/Classic/TrimCoils/TrimCoilBFit.cpp @@ -18,7 +18,6 @@ TrimCoilBFit::TrimCoilBFit(double bmax, void TrimCoilBFit::doApplyField(const double r, const double z, double *br, double *bz) { - if (std::abs(bmax_m) < 1e-20) return; // check range if (r < rmin_m || r > rmax_m) return; diff --git a/src/Classic/TrimCoils/TrimCoilMirrored.cpp b/src/Classic/TrimCoils/TrimCoilMirrored.cpp index 2f9d3bdf5..c690980f3 100644 --- a/src/Classic/TrimCoils/TrimCoilMirrored.cpp +++ b/src/Classic/TrimCoils/TrimCoilMirrored.cpp @@ -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/TrimCoilPhaseFit.cpp b/src/Classic/TrimCoils/TrimCoilPhaseFit.cpp index 693f0fbb8..66c71d730 100644 --- a/src/Classic/TrimCoils/TrimCoilPhaseFit.cpp +++ b/src/Classic/TrimCoils/TrimCoilPhaseFit.cpp @@ -18,7 +18,6 @@ TrimCoilPhaseFit::TrimCoilPhaseFit(double bmax, void TrimCoilPhaseFit::doApplyField(const double r, const double z, double *br, double *bz) { - if (std::abs(bmax_m) < 1e-20) return; // check range if (r < rmin_m || r > rmax_m) return; -- GitLab From 3e07610a81a3d98557d8ec829cac6ba5bf0dbc7c Mon Sep 17 00:00:00 2001 From: Jochem Snuverink <jochem.snuverink@psi.ch> Date: Fri, 8 Feb 2019 17:34:32 +0100 Subject: [PATCH 2/4] add interface for azimuthal coils --- src/Classic/AbsBeamline/Cyclotron.cpp | 14 ++++---- src/Classic/AbsBeamline/Cyclotron.h | 7 ++-- src/Classic/TrimCoils/OpalTrimCoil.cpp | 37 ++++++++++++++-------- src/Classic/TrimCoils/TrimCoil.cpp | 20 ++++++++++-- src/Classic/TrimCoils/TrimCoil.h | 13 ++++++-- src/Classic/TrimCoils/TrimCoilBFit.cpp | 2 +- src/Classic/TrimCoils/TrimCoilBFit.h | 2 +- src/Classic/TrimCoils/TrimCoilMirrored.cpp | 2 +- src/Classic/TrimCoils/TrimCoilMirrored.h | 2 +- src/Classic/TrimCoils/TrimCoilPhaseFit.cpp | 2 +- src/Classic/TrimCoils/TrimCoilPhaseFit.h | 2 +- 11 files changed, 68 insertions(+), 35 deletions(-) diff --git a/src/Classic/AbsBeamline/Cyclotron.cpp b/src/Classic/AbsBeamline/Cyclotron.cpp index 41f1792e8..dfd70ac70 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 d27dd5ec6..7edebe04f 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/OpalTrimCoil.cpp b/src/Classic/TrimCoils/OpalTrimCoil.cpp index edff5b06f..cc9f38ebb 100644 --- a/src/Classic/TrimCoils/OpalTrimCoil.cpp +++ b/src/Classic/TrimCoils/OpalTrimCoil.cpp @@ -22,6 +22,8 @@ namespace { COEFNUM, // COEFDENOM, // BMAX, // + PHIMIN, + PHIMAX, RMIN, // RMAX, // SLPTC, @@ -45,15 +47,21 @@ OpalTrimCoil::OpalTrimCoil(): 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 +84,6 @@ OpalTrimCoil::OpalTrimCoil(const std::string &name, OpalTrimCoil *parent): OpalTrimCoil::~OpalTrimCoil() { - } @@ -114,13 +121,15 @@ 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]); @@ -136,7 +145,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,7 +155,7 @@ 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]); @@ -154,7 +165,7 @@ Inform& OpalTrimCoil::print(Inform &os) const { << ((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) { diff --git a/src/Classic/TrimCoils/TrimCoil.cpp b/src/Classic/TrimCoils/TrimCoil.cpp index 51be9ebf9..99de80115 100644 --- a/src/Classic/TrimCoils/TrimCoil.cpp +++ b/src/Classic/TrimCoils/TrimCoil.cpp @@ -2,6 +2,8 @@ #include <cmath> +#include "Physics/Physics.h" + TrimCoil::TrimCoil(double bmax, double rmin, double rmax) @@ -14,9 +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); +} - doApplyField(r,z,br,bz); -} \ No newline at end of file +void TrimCoil::setAzimuth(const double phimin, const double phimax) +{ + // 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 35647ba5e..254a868a4 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 197e2dbca..1c1321788 100644 --- a/src/Classic/TrimCoils/TrimCoilBFit.cpp +++ b/src/Classic/TrimCoils/TrimCoilBFit.cpp @@ -16,7 +16,7 @@ TrimCoilBFit::TrimCoilBFit(double bmax, coefdenom_m.push_back(1.0); } -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) { // check range if (r < rmin_m || r > rmax_m) return; diff --git a/src/Classic/TrimCoils/TrimCoilBFit.h b/src/Classic/TrimCoils/TrimCoilBFit.h index 128ea9815..cb760dcf5 100644 --- a/src/Classic/TrimCoils/TrimCoilBFit.h +++ b/src/Classic/TrimCoils/TrimCoilBFit.h @@ -27,7 +27,7 @@ private: 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/TrimCoilMirrored.cpp b/src/Classic/TrimCoils/TrimCoilMirrored.cpp index c690980f3..3dea80876 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 diff --git a/src/Classic/TrimCoils/TrimCoilMirrored.h b/src/Classic/TrimCoils/TrimCoilMirrored.h index 7787e9fd1..6ce955f4a 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 66c71d730..c81de98ba 100644 --- a/src/Classic/TrimCoils/TrimCoilPhaseFit.cpp +++ b/src/Classic/TrimCoils/TrimCoilPhaseFit.cpp @@ -16,7 +16,7 @@ TrimCoilPhaseFit::TrimCoilPhaseFit(double bmax, coefdenom_m.push_back(1.0); } -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) { // check range if (r < rmin_m || r > rmax_m) return; diff --git a/src/Classic/TrimCoils/TrimCoilPhaseFit.h b/src/Classic/TrimCoils/TrimCoilPhaseFit.h index 22e17e7c4..a5898354b 100644 --- a/src/Classic/TrimCoils/TrimCoilPhaseFit.h +++ b/src/Classic/TrimCoils/TrimCoilPhaseFit.h @@ -26,7 +26,7 @@ private: 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 -- GitLab From 8bc0bc9b673f0fd8acb5d546df84bdc7be4b8082 Mon Sep 17 00:00:00 2001 From: Jochem Snuverink <jochem.snuverink@psi.ch> Date: Thu, 14 Mar 2019 09:44:06 +0100 Subject: [PATCH 3/4] small logic fix and add phi to unit tests --- src/Classic/TrimCoils/TrimCoil.cpp | 2 +- .../classic_src/AbsBeamline/TrimCoilTest.cpp | 88 +++++++++++++------ 2 files changed, 62 insertions(+), 28 deletions(-) diff --git a/src/Classic/TrimCoils/TrimCoil.cpp b/src/Classic/TrimCoils/TrimCoil.cpp index 99de80115..c34225e4d 100644 --- a/src/Classic/TrimCoils/TrimCoil.cpp +++ b/src/Classic/TrimCoils/TrimCoil.cpp @@ -20,7 +20,7 @@ void TrimCoil::applyField(const double r, const double z, const double phi_rad, { 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; + (phimin_m > phimax_m && (phi_rad < phimin_m && phi_rad > phimax_m)) ) return; doApplyField(r,z,phi_rad,br,bz); } diff --git a/tests/classic_src/AbsBeamline/TrimCoilTest.cpp b/tests/classic_src/AbsBeamline/TrimCoilTest.cpp index 16b7cc6ed..6f3315f60 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,32 @@ TEST(TrimCoil, TrimCoilBFitZeros) double bmax = 0.0; double rmin = 1000; // mm double rmax = 2000; + double phi = 0.0; // rad 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.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); } +// TrimCoilBFit class TEST(TrimCoil, TrimCoilBFit) { OpalTestUtilities::SilenceTest silencer; @@ -47,26 +51,52 @@ 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}, {}); - double br = 1.0, bz = 1.0; - myTrimCoil.applyField(2.0, 2.0, &br, &bz); + const double brStart = 1.0, bzStart = 1.0; + double br = brStart, bz = bzStart; + myTrimCoil.applyField(2.0, 2.0, 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); + br = brStart, bz = bzStart; + myTrimCoil.applyField(2.0, 2.0, 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, 2.0, 3.2, &br, &bz); // outside range + EXPECT_NEAR(bz, bzStart, margin); + myTrimCoil.applyField(2.0, 2.0, 0.0, &br, &bz); // outside range + EXPECT_NEAR(bz, bzStart, margin); + myTrimCoil.applyField(2.0, 2.0, 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, 2.0, 1.0, &br, &bz); // outside range + EXPECT_NEAR(bz, bzStart, margin); + myTrimCoil.applyField(2.0, 2.0, 3.2, &br, &bz); // inside range + EXPECT_NEAR(bz, bzSolution, margin); + myTrimCoil.applyField(2.0, 2.0, 0.0, &br, &bz); // inside range + EXPECT_NEAR(bz, 2*bzSolution - bzStart, margin); } +// TrimCoilPhaseFit class TEST(TrimCoil, TrimCoilPhaseFit) { OpalTestUtilities::SilenceTest silencer; @@ -74,26 +104,28 @@ 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); + myTrimCoil.applyField(2.0, 2.0, 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.applyField(2.0, 2.0, phi, &br, &bz); + // 1 - 10*(-0.2) = 3; 1 + 2*10*2*(-0.2*2) / 5 = -2.2 + EXPECT_NEAR(bz, 3.0, margin); + EXPECT_NEAR(br,-2.2, margin); } +// TrimCoilMirrored class TEST(TrimCoil, TrimCoilMirrored) { OpalTestUtilities::SilenceTest silencer; @@ -101,13 +133,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 +} + -- GitLab From 84f4b2fd5353fd332f8d208b807c028d41335c9f Mon Sep 17 00:00:00 2001 From: Jochem Snuverink <jochem.snuverink@psi.ch> Date: Thu, 14 Mar 2019 12:00:10 +0100 Subject: [PATCH 4/4] refactor TrimCoilBFit and PhaseFit with common base class; add options to shape field in phi --- src/Classic/TrimCoils/CMakeLists.txt | 2 + src/Classic/TrimCoils/OpalTrimCoil.cpp | 53 ++++---- src/Classic/TrimCoils/OpalTrimCoil.h | 3 + src/Classic/TrimCoils/TrimCoilBFit.cpp | 52 ++------ src/Classic/TrimCoils/TrimCoilBFit.h | 11 +- src/Classic/TrimCoils/TrimCoilFit.cpp | 118 ++++++++++++++++++ src/Classic/TrimCoils/TrimCoilFit.h | 40 ++++++ src/Classic/TrimCoils/TrimCoilPhaseFit.cpp | 72 ++--------- src/Classic/TrimCoils/TrimCoilPhaseFit.h | 11 +- .../classic_src/AbsBeamline/TrimCoilTest.cpp | 79 ++++++++---- 10 files changed, 284 insertions(+), 157 deletions(-) create mode 100644 src/Classic/TrimCoils/TrimCoilFit.cpp create mode 100644 src/Classic/TrimCoils/TrimCoilFit.h diff --git a/src/Classic/TrimCoils/CMakeLists.txt b/src/Classic/TrimCoils/CMakeLists.txt index fb649d834..5baaf194c 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 cc9f38ebb..448ab4199 100644 --- a/src/Classic/TrimCoils/OpalTrimCoil.cpp +++ b/src/Classic/TrimCoils/OpalTrimCoil.cpp @@ -21,6 +21,8 @@ namespace { TYPE, // The type of trim coil COEFNUM, // COEFDENOM, // + COEFNUMPHI, + COEFDENOMPHI, BMAX, // PHIMIN, PHIMAX, @@ -39,10 +41,16 @@ 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."); @@ -131,12 +139,14 @@ void OpalTrimCoil::initOpalTrimCoil() { 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]); @@ -158,30 +168,29 @@ Inform& OpalTrimCoil::print(Inform &os) const { 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 9fe782e3d..26728c1c1 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/TrimCoilBFit.cpp b/src/Classic/TrimCoils/TrimCoilBFit.cpp index 1c1321788..4129879b7 100644 --- a/src/Classic/TrimCoils/TrimCoilBFit.cpp +++ b/src/Classic/TrimCoils/TrimCoilBFit.cpp @@ -1,57 +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, const double phi_rad, double *br, double *bz) { // 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 cb760dcf5..b734cdadd 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,23 +9,22 @@ /// 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, const double phi_rad, double *br, double *bz); }; diff --git a/src/Classic/TrimCoils/TrimCoilFit.cpp b/src/Classic/TrimCoils/TrimCoilFit.cpp new file mode 100644 index 000000000..6b671fef4 --- /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 000000000..be23d9be8 --- /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/TrimCoilPhaseFit.cpp b/src/Classic/TrimCoils/TrimCoilPhaseFit.cpp index c81de98ba..679618282 100644 --- a/src/Classic/TrimCoils/TrimCoilPhaseFit.cpp +++ b/src/Classic/TrimCoils/TrimCoilPhaseFit.cpp @@ -1,75 +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, const double phi_rad, double *br, double *bz) { // 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 a5898354b..dac12ed6a 100644 --- a/src/Classic/TrimCoils/TrimCoilPhaseFit.h +++ b/src/Classic/TrimCoils/TrimCoilPhaseFit.h @@ -1,30 +1,29 @@ #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, const double phi_rad, double *br, double *bz); }; diff --git a/tests/classic_src/AbsBeamline/TrimCoilTest.cpp b/tests/classic_src/AbsBeamline/TrimCoilTest.cpp index 6f3315f60..2030d212b 100644 --- a/tests/classic_src/AbsBeamline/TrimCoilTest.cpp +++ b/tests/classic_src/AbsBeamline/TrimCoilTest.cpp @@ -20,7 +20,7 @@ TEST(TrimCoil, TrimCoilBFitZeros) 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; @@ -30,7 +30,7 @@ TEST(TrimCoil, TrimCoilBFitZeros) EXPECT_NEAR(bz, one, margin); bmax = 1.0; - myTrimCoil = TrimCoilBFit(bmax, rmin, rmax, {}, {}); + myTrimCoil = TrimCoilBFit(bmax, rmin, rmax, {}, {}, {}, {}); myTrimCoil.applyField(rmin*mm2m - 1, 1, phi, &br, &bz); // not changed since r outside range @@ -41,6 +41,11 @@ TEST(TrimCoil, TrimCoilBFitZeros) // not changed since r outside range 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 @@ -54,20 +59,20 @@ TEST(TrimCoil, TrimCoilBFit) 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}, {}, {}, {}); - const double brStart = 1.0, bzStart = 1.0; + const double brStart = 1.0, bzStart = 1.0, zStart = 2.0; double br = brStart, bz = bzStart; - myTrimCoil.applyField(2.0, 2.0, phi, &br, &bz); + 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, 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}); + myTrimCoil = TrimCoilBFit(bmax, rmin, rmax, {4.0, 3.0}, {1.0, 2.0}, {}, {}); br = brStart, bz = bzStart; - myTrimCoil.applyField(2.0, 2.0, phi, &br, &bz); + 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) @@ -78,22 +83,36 @@ TEST(TrimCoil, TrimCoilBFit) // test phi angles myTrimCoil.setAzimuth(10,180); br = brStart, bz = bzStart; - myTrimCoil.applyField(2.0, 2.0, 3.2, &br, &bz); // outside range + myTrimCoil.applyField(2.0, zStart, 3.2, &br, &bz); // outside range EXPECT_NEAR(bz, bzStart, margin); - myTrimCoil.applyField(2.0, 2.0, 0.0, &br, &bz); // outside range + myTrimCoil.applyField(2.0, zStart, 0.0, &br, &bz); // outside range EXPECT_NEAR(bz, bzStart, margin); - myTrimCoil.applyField(2.0, 2.0, 1.0, &br, &bz); // inside range + 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, 2.0, 1.0, &br, &bz); // outside range + myTrimCoil.applyField(2.0, zStart, 1.0, &br, &bz); // outside range EXPECT_NEAR(bz, bzStart, margin); - myTrimCoil.applyField(2.0, 2.0, 3.2, &br, &bz); // inside range + myTrimCoil.applyField(2.0, zStart, 3.2, &br, &bz); // inside range EXPECT_NEAR(bz, bzSolution, margin); - myTrimCoil.applyField(2.0, 2.0, 0.0, &br, &bz); // inside range + 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 @@ -107,22 +126,36 @@ TEST(TrimCoil, TrimCoilPhaseFit) 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, phi, &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, 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}); + 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, phi, &br, &bz); - // 1 - 10*(-0.2) = 3; 1 + 2*10*2*(-0.2*2) / 5 = -2.2 - EXPECT_NEAR(bz, 3.0, margin); - EXPECT_NEAR(br,-2.2, margin); + 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 -- GitLab