Resolve "Specify trim coil range azimuthally"

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() {

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:
     

src/Classic/TrimCoils/CMakeLists.txt

 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

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

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) {

src/Classic/TrimCoils/TrimCoil.cpp

 #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;
+}

src/Classic/TrimCoils/TrimCoil.h

 #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

src/Classic/TrimCoils/TrimCoilBFit.cpp

 #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

src/Classic/TrimCoils/TrimCoilBFit.h

 #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

src/Classic/TrimCoils/TrimCoilFit.cpp

+#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

src/Classic/TrimCoils/TrimCoilFit.h

+#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

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;

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;
 };

src/Classic/TrimCoils/TrimCoilPhaseFit.cpp

 #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

src/Classic/TrimCoils/TrimCoilPhaseFit.h

 #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

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
+}
+