Resolve "Trimcoil implementation bug"

src/Classic/TrimCoils/OpalTrimCoil.cpp

@@ -76,10 +76,10 @@ OpalTrimCoil::OpalTrimCoil():
                          ("BMAX", "Maximum magnetic field in Tesla.");
 
     itsAttr[PHIMIN]    = Attributes::makeReal
-                         ("PHIMIN", "Minimal azimuth [deg] (default 0)");
+                         ("PHIMIN", "Minimal azimuth [deg] (default 0)", 0.0);
 
     itsAttr[PHIMAX]    = Attributes::makeReal
-                         ("PHIMAX", "Maximal azimuth [deg] (default 360)");
+                         ("PHIMAX", "Maximal azimuth [deg] (default 360)", 360.0);
 
     itsAttr[RMIN]      = Attributes::makeReal
                          ("RMIN", "Minimum radius [mm].");
@@ -191,10 +191,11 @@ Inform& OpalTrimCoil::print(Inform &os) const {
         printPolynom(os,itsAttr[COEFDENOMPHI]);
     }
 
-    os << "* BMAX           " << Attributes::getReal(itsAttr[BMAX]) << '\n'
-       << "* RMIN           " << Attributes::getReal(itsAttr[RMIN]) << '\n'
-       << "* RMAX           " << Attributes::getReal(itsAttr[RMAX]) << '\n';
-
+    os << "* BMAX           " << Attributes::getReal(itsAttr[BMAX])   << '\n'
+       << "* RMIN           " << Attributes::getReal(itsAttr[RMIN])   << '\n'
+       << "* RMAX           " << Attributes::getReal(itsAttr[RMAX])   << '\n'
+       << "* PHIMIN         " << Attributes::getReal(itsAttr[PHIMIN]) << '\n'
+       << "* PHIMAX         " << Attributes::getReal(itsAttr[PHIMAX]) << '\n';
     if (Attributes::getString(itsAttr[TYPE]) == "PSI-BFIELD-MIRRORED") {
         os << "* SLPTC          " << Attributes::getReal(itsAttr[SLPTC]) << '\n';
     }

src/Classic/TrimCoils/TrimCoil.cpp

@@ -27,34 +27,32 @@
 #include <cmath>
 
 #include "Physics/Units.h"
+#include "Utilities/Util.h"
 
 TrimCoil::TrimCoil(double bmax,
                    double rmin,
                    double rmax)
 {
     // convert to m
-    const double mm2m = 0.001;
-    rmin_m = rmin * mm2m;
-    rmax_m = rmax * mm2m;
+    rmin_m = rmin * Units::mm2m;
+    rmax_m = rmax * Units::mm2m;
     // convert to kG
-    bmax_m = bmax * 10.0;
+    bmax_m = bmax * Units::T2kG;
 }
 
 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);
+    double phi = Util::angle_0to2pi(phi_rad);
+    // check if phi is inside [phimin_m, phimax_m]
+    if (phimin_m == phimax_m || Util::angleBetweenAngles(phi, phimin_m, phimax_m))
+        doApplyField(r,z,phi_rad,br,bz);
 }
 
 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 * Units::deg2rad;
-    phimax_m = phimax * Units::deg2rad;
+    // phi convert to rad
+    phimin_m = Util::angle_0to2pi(phimin * Units::deg2rad);
+    phimax_m = Util::angle_0to2pi(phimax * Units::deg2rad);
 }
\ No newline at end of file

src/Classic/TrimCoils/TrimCoil.h

@@ -35,7 +35,7 @@ public:
     /// 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, const double phi_rad, double *br, double *bz);
-    /// Set azimuthal range
+    /// Set azimuthal range where trim coil acts: [phimin, phimax] (also when phimin > phimax)
     void setAzimuth(const double phimin, const double phimax);
     virtual ~TrimCoil() { };
 

src/Classic/Utilities/Util.h

@@ -21,6 +21,8 @@
 #include "Algorithms/Vektor.h"
 #include "Algorithms/Quaternion.h"
 
+#include "Physics/Physics.h"
+
 #include <algorithm>
 #include <cmath>
 #include <cstring>
@@ -194,6 +196,19 @@ namespace Util {
 
     Vector_t getTaitBryantAngles(Quaternion rotation, const std::string& elementName = "");
 
+    /// convert angle (in rad) to [0,2pi) range, from https://stackoverflow.com/a/29721295
+    inline double angle_0to2pi(double angle) {
+        // converts angle to range [-2*pi, 2*pi)
+        angle = std::fmod(angle, Physics::two_pi);
+        if (angle >= 0.0) return angle;
+        else              return angle + Physics::two_pi;
+    }
+    /// check if angle (in rad and in range [0,2pi]) is within [min, max]
+    inline bool angleBetweenAngles(const double angle, const double min, const double max) {
+        if (min <= max) return (angle >= min && angle <= max);
+        else            return (angle >= min || angle <= max);
+    }
+
     std::string toUpper(const std::string& str);
 
     std::string boolToUpperString(const bool& b);

tests/classic_src/AbsBeamline/TrimCoilTest.cpp

@@ -4,6 +4,9 @@
 #include "TrimCoils/TrimCoilPhaseFit.h"
 #include "TrimCoils/TrimCoilMirrored.h"
 
+#include "Physics/Units.h"
+#include "Utilities/Util.h"
+
 #include "opal_test_utilities/SilenceTest.h"
 
 const double margin = 1e-7;
@@ -13,8 +16,6 @@ TEST(TrimCoil, TrimCoilBFitZeros)
 {
     OpalTestUtilities::SilenceTest silencer;
 
-    const double mm2m = 0.001;
-
     double bmax = 0.0;
     double rmin = 1000; // mm
     double rmax = 2000;
@@ -24,7 +25,7 @@ TEST(TrimCoil, TrimCoilBFitZeros)
 
     const double one = 1.0;
     double br = one, bz = one;
-    myTrimCoil.applyField((rmin+rmax)*mm2m/2.0, 1, 0, &br, &bz);
+    myTrimCoil.applyField((rmin + rmax) * Units::mm2m / 2.0, 1, 0, &br, &bz);
     // not changed since bmax 0.0
     EXPECT_NEAR(br, one, margin);
     EXPECT_NEAR(bz, one, margin);
@@ -32,17 +33,17 @@ TEST(TrimCoil, TrimCoilBFitZeros)
     bmax = 1.0;
     myTrimCoil = TrimCoilBFit(bmax, rmin, rmax, {}, {}, {}, {});
 
-    myTrimCoil.applyField(rmin*mm2m - 1, 1, phi, &br, &bz);
+    myTrimCoil.applyField(rmin * Units::mm2m - 1, 1, phi, &br, &bz);
     // 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);
+    myTrimCoil.applyField(rmax * Units::mm2m + 1, 1, phi, &br, &bz);
     // 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);
+    myTrimCoil.applyField(rmax * Units::mm2m - 1, 1, phi, &br, &bz);
     // default constant field
     EXPECT_NEAR(br, one, margin);
     EXPECT_NEAR(bz, 11.0, margin); // 1 + 10*1
@@ -80,8 +81,19 @@ TEST(TrimCoil, TrimCoilBFit)
     EXPECT_NEAR(bz, bzSolution, margin);
     EXPECT_NEAR(br, brSolution, margin);
 
+    phi = Util::angle_0to2pi(1.0);
+    EXPECT_NEAR(phi, 1, margin);
+    phi = Util::angle_0to2pi(-1.0);
+    EXPECT_NEAR(phi, -1+Physics::two_pi, margin);
+
     // test phi angles
     myTrimCoil.setAzimuth(10,180);
+    bool valid = Util::angleBetweenAngles(0.0, 10*Units::deg2rad, 180*Units::deg2rad);
+    EXPECT_FALSE(valid);
+    valid = Util::angleBetweenAngles(1.0, 10*Units::deg2rad, 180*Units::deg2rad);
+    EXPECT_TRUE(valid);
+    valid = Util::angleBetweenAngles(3.2, 10*Units::deg2rad, 180*Units::deg2rad);
+    EXPECT_FALSE(valid);
     br = brStart, bz = bzStart;
     myTrimCoil.applyField(2.0, zStart, 3.2, &br, &bz); // outside range
     EXPECT_NEAR(bz, bzStart, margin);
@@ -92,6 +104,12 @@ TEST(TrimCoil, TrimCoilBFit)
 
     // test phi angles: first larger
     myTrimCoil.setAzimuth(180,20);
+    valid = Util::angleBetweenAngles(0.0, 180*Units::deg2rad, 20*Units::deg2rad);
+    EXPECT_TRUE(valid);
+    valid = Util::angleBetweenAngles(1.0, 180*Units::deg2rad, 20*Units::deg2rad);
+    EXPECT_FALSE(valid);
+    valid = Util::angleBetweenAngles(3.2, 180*Units::deg2rad, 20*Units::deg2rad);
+    EXPECT_TRUE(valid);
     br = brStart, bz = bzStart;
     myTrimCoil.applyField(2.0, zStart, 1.0, &br, &bz); // outside range
     EXPECT_NEAR(bz, bzStart, margin);
@@ -100,7 +118,24 @@ TEST(TrimCoil, TrimCoilBFit)
     myTrimCoil.applyField(2.0, zStart, 0.0, &br, &bz); //  inside range
     EXPECT_NEAR(bz, 2*bzSolution - bzStart, margin);
 
-    // Test Phi
+    // test phi angles: negative
+    myTrimCoil.setAzimuth(10,-180);
+    double pi = Util::angle_0to2pi(-180*Units::deg2rad);
+    valid = Util::angleBetweenAngles(0.0, 10*Units::deg2rad, pi);
+    EXPECT_FALSE(valid);
+    valid = Util::angleBetweenAngles(1.0, 10*Units::deg2rad, pi);
+    EXPECT_TRUE(valid);
+    valid = Util::angleBetweenAngles(3.2, 10*Units::deg2rad, pi);
+    EXPECT_FALSE(valid);
+    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 function
     // same rational function
     myTrimCoil = TrimCoilBFit(bmax, rmin, rmax, {}, {}, {4.0, 3.0}, {1.0, 2.0});
     br = brStart, bz = bzStart;