Maybe the maths is right and implemented; but Maxwell still not obeyed; more work needed

src/Classic/AbsBeamline/CMakeLists.txt

+add_subdirectory(EndFieldModel)
+
 set (_SRCS
   AlignWrapper.cpp
   AttributeSet.cpp
@@ -35,6 +37,7 @@ set (_SRCS
   Septum.cpp
   Solenoid.cpp
   Source.cpp
+  SpiralSector.cpp
   Stripper.cpp
   TravelingWave.cpp
   VariableRFCavity.cpp
@@ -85,9 +88,10 @@ set (HDRS
     Solenoid.h
     Source.h
     SpecificElementVisitor.h
+    SpiralSector.h
     Stripper.h
     TravelingWave.h
     VariableRFCavity.h
 )
 
-install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/AbsBeamline")
\ No newline at end of file
+install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/AbsBeamline")

src/Classic/AbsBeamline/EndFieldModel/CMakeLists.txt

 set (_SRCS
-  Tanh.cpp 
+  EndFieldModel.cpp
   Enge.cpp
   Sin.cpp
-  EndFieldModel.cpp
+  Tanh.cpp 
   )
 
 include_directories (
@@ -12,10 +12,10 @@ include_directories (
 ADD_OPAL_SOURCES(${_SRCS})
 
 set (HDRS
-  Tanh.h
+  EndFieldModel.h
   Enge.h
   Sin.h
-  EndFieldModel.h
+  Tanh.h
 )
 
 install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/AbsBeamline/EndFieldModel/")

src/Classic/AbsBeamline/EndFieldModel/EndFieldModel.h

@@ -25,7 +25,8 @@
  *  POSSIBILITY OF SUCH DAMAGE.
  */
 
-#ifndef ENDFIELDMODEL_H_
+#ifndef ENDFIELDMODEL_ENDFIELDMODEL_H_
+#define ENDFIELDMODEL_ENDFIELDMODEL_H_
 
 #include <iostream>
 #include <vector>
@@ -35,9 +36,9 @@ namespace endfieldmodel {
 class EndFieldModel {
  public:
   virtual ~EndFieldModel() {;}
-  virtual std::ostream& Print(std::ostream& out) const = 0;
-  virtual double Function(double x, int n) const = 0;
-  virtual EndFieldModel* Clone() const = 0;
+  virtual std::ostream& print(std::ostream& out) const = 0;
+  virtual double function(double x, int n) const = 0;
+  virtual EndFieldModel* clone() const = 0;
 };
 
 std::vector< std::vector<int> > CompactVector

src/Classic/AbsBeamline/EndFieldModel/Enge.cpp

@@ -36,7 +36,7 @@ namespace endfieldmodel {
 
 // Use
 // d^n E/dx^n = a_n1m1 F(n1) g(m1) + a_n2m1m2 F(n2) g(m1)g(m2)+...
-// where
+// where 
 double Enge::GetEnge(double x, int n) const {
   std::vector< std::vector<int> > qt = GetQIndex(n);
   std::vector<double> g;

src/Classic/AbsBeamline/EndFieldModel/Enge.h

@@ -26,6 +26,7 @@
  */
 
 #ifndef ENDFIELDMODEL_ENGE_H_
+#define ENDFIELDMODEL_ENGE_H_
 
 #include <iostream>
 #include <vector>

src/Classic/AbsBeamline/EndFieldModel/Tanh.cpp

@@ -37,10 +37,17 @@ namespace endfieldmodel {
 std::vector< std::vector< std::vector<int> > > Tanh::_tdi;
 
 Tanh::Tanh(double x0, double lambda, int max_index) : _x0(x0), _lambda(lambda) {
-  SetTanhDiffIndices(max_index);
+  setTanhDiffIndices(max_index);
 }
 
-double Tanh::GetTanh(double x, int n) const {
+Tanh::~Tanh() {}
+
+EndFieldModel* Tanh::clone() const {
+    return new Tanh(*this);
+}
+
+
+double Tanh::getTanh(double x, int n) const {
   if (n == 0) return tanh((x+_x0)/_lambda);
   double t = 0;
   double lam_n = gsl_sf_pow_int(_lambda, n);
@@ -51,7 +58,7 @@ double Tanh::GetTanh(double x, int n) const {
   return t;
 }
 
-double Tanh::GetNegTanh(double x, int n) const {
+double Tanh::getNegTanh(double x, int n) const {
   if (n == 0) return tanh((x-_x0)/_lambda);
   double t = 0;
   double lam_n = gsl_sf_pow_int(_lambda, n);
@@ -62,11 +69,11 @@ double Tanh::GetNegTanh(double x, int n) const {
   return t;
 }
 
-double Tanh::GetDoubleTanh(double x, int n) const {
-  return (GetTanh(x, n)-GetNegTanh(x, n))/2.;
+double Tanh::function(double x, int n) const {
+  return (getTanh(x, n)-getNegTanh(x, n))/2.;
 }
 
-void Tanh::SetTanhDiffIndices(size_t n) {
+void Tanh::setTanhDiffIndices(size_t n) {
   if (_tdi.size() == 0) {
     _tdi.push_back(std::vector< std::vector<int> >(1, std::vector<int>(2)));
     _tdi[0][0][0] = 1;  // 1*tanh(x) - third index is redundant
@@ -92,8 +99,8 @@ void Tanh::SetTanhDiffIndices(size_t n) {
   }
 }
 
-std::vector< std::vector<int> > Tanh::GetTanhDiffIndices(size_t n) {
-  SetTanhDiffIndices(n);
+std::vector< std::vector<int> > Tanh::getTanhDiffIndices(size_t n) {
+  setTanhDiffIndices(n);
   return _tdi[n];
 }
 

src/Classic/AbsBeamline/EndFieldModel/Tanh.h

@@ -26,6 +26,7 @@
  */
 
 #ifndef ENDFIELDMODEL_TANH_H_
+#define ENDFIELDMODEL_TANH_H_
 
 #include <iostream>
 #include <vector>
@@ -55,46 +56,52 @@ class Tanh : public EndFieldModel {
     Tanh(double x0, double lambda, int max_index);
 
     /** Default constructor (initialises x0 and lambda to 0) */
-    Tanh() : _x0(0.), _lambda(0.) {SetTanhDiffIndices(10);}
+    Tanh() : _x0(0.), _lambda(0.) {setTanhDiffIndices(10);}
+
+    /** Copy constructor */
+    Tanh(const Tanh& rhs) : _x0(rhs._x0), _lambda(rhs._lambda) {}
 
     /** Destructor (no mallocs so does nothing) */
-    ~Tanh() {}
+    ~Tanh();
 
-    /** Returns the value of DoubleTanh or its \f$n^{th}\f$ derivative. */
-    Tanh* Clone() const;
+    /** Inherited copy constructor. */
+    EndFieldModel* clone() const;
 
     /** Double Tanh is given by\n
      *  \f$d(x) = \f$
      */
-    double GetDoubleTanh(double x, int n) const;
+    double function(double x, int n) const;
 
     /** Returns the value of tanh((x+x0)/lambda) or its \f$n^{th}\f$ derivative. */
-    double GetTanh(double x, int n) const;
+    double getTanh(double x, int n) const;
 
     /** Returns the value of tanh((x-x0)/lambda) or its \f$n^{th}\f$ derivative. */
-    double GetNegTanh(double x, int n) const;
+    double getNegTanh(double x, int n) const;
 
     /** Get all the tanh differential indices \f$I_{pq}\f$.
      *
      *  Returns vector of vector of ints where p indexes the differential and
      * q indexes the tanh power - so
      */
-    static std::vector< std::vector<int> > GetTanhDiffIndices(size_t n);
+    static std::vector< std::vector<int> > getTanhDiffIndices(size_t n);
 
     /** Set the value of tanh differential indices to nth order differentials. */
-    static void SetTanhDiffIndices(size_t n);
+    static void setTanhDiffIndices(size_t n);
 
     /** Return lambda (end length) */
-    inline double GetLambda() const {return _lambda;}
+    inline double getLambda() const {return _lambda;}
 
     /** Return x0 (flat top length) */
-    inline double GetX0() const {return _x0;}
+    inline double getX0() const {return _x0;}
 
     /** Set lambda (end length) */
-    inline void   SetLambda(double lambda) {_lambda = lambda;}
+    inline void setLambda(double lambda) {_lambda = lambda;}
 
     /** Set x0 (flat top length) */
-    inline void   SetX0(double x0)     {_x0 = x0;}
+    inline void setX0(double x0)     {_x0 = x0;}
+
+    /** Bug - does nothing */
+    std::ostream& print(std::ostream& out) const {return out;}
   private:
     double _x0, _lambda;
 

src/Classic/AbsBeamline/SpiralSector.cpp

@@ -25,7 +25,7 @@
  *  POSSIBILITY OF SUCH DAMAGE.
  */
 
-#include <math.h>
+#include <cmath>
 
 #include "AbsBeamline/SpiralSector.h"
 #include "Algorithms/PartBunch.h"
@@ -51,7 +51,6 @@ ElementBase* SpiralSector::clone() const {
     return new SpiralSector(*this);
 }
 
-
 EMField &SpiralSector::getField() {
     return dummy;
 }
@@ -67,6 +66,7 @@ bool SpiralSector::apply(const size_t &i, const double &t,
 
 void SpiralSector::initialise(PartBunch *bunch, double &startField, double &endField) {
     RefPartBunch_m = bunch;
+    calculateDfCoefficients();
 }
 
 void SpiralSector::finalise() {
@@ -93,11 +93,86 @@ void SpiralSector::accept(BeamlineVisitor& visitor) const {
 bool SpiralSector::apply(const Vector_t &R, const Vector_t &P,
                     const double &t, Vector_t &E, Vector_t &B) {
     Vector_t pos = R - centre_m;
-    double radius = sqrt(pos[0]*pos[0]+pos[2]*pos[2]);
+    double r = sqrt(pos[0]*pos[0]+pos[2]*pos[2]);
     double phi = atan2(pos[2], pos[0]);
-    double fringe = getFringe(phi);
-    B[1] = fringe*Bz_m*pow(radius/r0_m, k_m);
+    Vector_t posCyl(r, pos[1], phi);
+    Vector_t bCyl(0., 0., 0.);
+    applyCylindrical(posCyl, P, t, E, bCyl);
+    // this is cartesian coordinates
+    B[1] = bCyl[1];
+    B[0] = bCyl[0]*cos(phi) + bCyl[2]*sin(phi);
+    B[2] = bCyl[2]*cos(phi) + bCyl[0]*sin(phi);
+    return true;
+
+}
+
+void SpiralSector::calculateDfCoefficients() {
+    dfCoefficients_m = std::vector<std::vector<double> >(maxOrder_m);
+    dfCoefficients_m[0] = std::vector<double>(1, 1.); // f_0 = 1.*0th derivative
+    for (size_t n = 0; n+1 < maxOrder_m; n += 2) { // n indexes the power in z
+        dfCoefficients_m[n+1] = std::vector<double>(dfCoefficients_m[n].size()+1, 0);
+        for (size_t i = 0; i < dfCoefficients_m[n].size(); ++i) { // i indexes the derivative
+            dfCoefficients_m[n+1][i+1] = dfCoefficients_m[n][i]/(2*n+1);
+        }
+        if (n+2 == maxOrder_m) {
+            break;
+        }
+        dfCoefficients_m[n+2] = std::vector<double>(dfCoefficients_m[n].size()+2, 0);
+        for(size_t i = 0; i < dfCoefficients_m[n].size(); ++i) { // i indexes the derivative
+            dfCoefficients_m[n+2][i] = -(k_m-2*n)*(k_m-2*n)/(n+1)*dfCoefficients_m[n][i];
+        }
+        for(size_t i = 0; i < dfCoefficients_m[n+1].size(); ++i) { // i indexes the derivative
+            dfCoefficients_m[n+2][i] += 2*(k_m-2*n)*tanDelta_m*dfCoefficients_m[n+1][i];
+            dfCoefficients_m[n+2][i+1] -= (1+tanDelta_m*tanDelta_m)*dfCoefficients_m[n+1][i];
+        }
+    }
+
+}
 
+BUG - IT COMPILES AND RUNS; I BELIEVE THE MATHS IS RIGHT; BUT MAXWELL DOES NOT IMPROVE!
+
+STRATEGY... Try Testing M2a-c; validate each of the field elements in turn; then consider M1 again
+
+bool SpiralSector::applyCylindrical(const Vector_t &pos, const Vector_t &P,
+                    const double &t, Vector_t &E, Vector_t &B) {
+    double r = pos[0];
+    double phi = pos[2];
+    double normRadius = r/r0_m;
+    double g = tanDelta_m*log(normRadius);
+    double phiSpiral = phi - g;
+    double h = pow(normRadius, k_m)*Bz_m;
+    std::vector<double> fringeDerivatives(maxOrder_m, 0.);
+    std::cerr << "(r, y, phi) " << pos << " h " << h << " psi " << phiSpiral << std::endl;
+    for (size_t i = 0; i < fringeDerivatives.size(); ++i) {
+        fringeDerivatives[i] = endField_m->function(phiSpiral, i);
+        std::cerr << fringeDerivatives[i] << " ";
+    }
+    std::cerr << std::endl;
+    for (size_t n = 0; n < dfCoefficients_m.size(); n += 2) {
+        double f2n = 0;
+        Vector_t deltaB;
+        for (size_t i = 0; i < dfCoefficients_m[n].size(); ++i) {
+            f2n += dfCoefficients_m[n][i]*fringeDerivatives[i];
+        }
+        double f2nplus1 = 0;
+        for (size_t i = 0; i < dfCoefficients_m[n+1].size() && n+1 < dfCoefficients_m.size(); ++i) {
+            f2nplus1 += dfCoefficients_m[n+1][i]*fringeDerivatives[i];
+        }
+        deltaB[1] = f2n*h*pow(pos[1]/r, n); // Bz = sum(f_2n * h * (z/r)^2n
+        deltaB[2] = f2nplus1*h*pow(pos[1]/r, n+1); // Bphi = sum(f_2n+1 * h * (z/r)^2n+1
+        deltaB[0] = (f2n*(k_m-n)/(n+1) - tanDelta_m*f2nplus1)*h*pow(pos[1]/r, n+1);
+        std::cerr << "    " << n << " " << deltaB << std::endl;
+        B += deltaB;
+    }
+    std::cerr << B << std::endl;
     return true;
 }
 
+void SpiralSector::setEndField(endfieldmodel::EndFieldModel* endField) {
+    if (endField_m != NULL) {
+        delete endField_m;
+    }
+
+    endField_m = endField;
+}
+

src/Classic/AbsBeamline/SpiralSector.h

@@ -44,9 +44,6 @@ class SpiralSector : public Component {
      */
     explicit SpiralSector(const std::string &name);
 
-    /** Copy constructor */
-    SpiralSector(const SpiralSector &right);
-
     /** Destructor - deletes map */
     ~SpiralSector();
 
@@ -67,6 +64,7 @@ class SpiralSector : public Component {
     /** Calculate the field at some arbitrary position
      *
      *  \param R position in the local coordinate system of the bend
+     *  \param P not used
      *  \param t time at which the field is to be calculated
      *  \param E calculated electric field - always 0 (no E-field)
      *  \param B calculated magnetic field
@@ -75,6 +73,19 @@ class SpiralSector : public Component {
     bool apply(const Vector_t &R, const Vector_t &P, const double &t,
                Vector_t &E, Vector_t &B);
 
+    /** Calculate the field at some arbitrary position in cylindrical coordinates
+     *
+     *  \param R position in the local coordinate system of the bend, in
+     *           cylindrical polar coordinates defined like (r, y, phi)
+     *  \param P not used
+     *  \param t not used (field is time independent)
+     *  \param E not used (no E-field)
+     *  \param B calculated magnetic field defined like (Br, By, Bphi)
+     *  \returns true if particle is outside the field map, else false
+     */
+    bool applyCylindrical(const Vector_t &R, const Vector_t &P,
+                    const double &t, Vector_t &E, Vector_t &B);
+
      /** Initialise the SpiralSector
       *
       *  \param bunch the global bunch object
@@ -107,6 +118,36 @@ class SpiralSector : public Component {
     /** Accept a beamline visitor */
     void accept(BeamlineVisitor& visitor) const;
 
+    /** Get tan delta - delta is the spiral angle */
+    Vector_t getTanDelta() const {return tanDelta_m;}
+
+    /** Set tan delta - delta is the spiral angle */
+    void setTanDelta(double tanDelta) {tanDelta_m = tanDelta;}
+
+    /** Get the field index k */
+    double getFieldIndex() const {return k_m;}
+
+    /** Set the field index k */
+    void setFieldIndex(double k) {k_m = k;}
+
+    /** Get the dipole constant B_0 */
+    double getDipoleConstant() const {return Bz_m;}
+
+    /** Set the dipole constant B_0 */
+    void setDipoleConstant(double Bz) {Bz_m = Bz;}
+
+    /** Get the radius constant R_0 */
+    double getR0() const {return r0_m;}
+
+    /** Set the radius constant R_0 */
+    void setR0(double r0) {r0_m = r0;}
+
+    /** Get the centre of the sector */
+    Vector_t getCentre() const {return centre_m;}
+
+    /** Set the centre of the sector */
+    void setCentre(Vector_t centre) {centre_m = centre;}
+
     /** Get the fringe field
      *
      *  Returns the fringe field model; SpiralSector retains ownership of the
@@ -121,27 +162,40 @@ class SpiralSector : public Component {
       */
     void setEndField(endfieldmodel::EndFieldModel* endField);
 
-    /** Get the centre of the sector */
-    Vector_t getCentre() const {return centre_m;}
+    /** Get the maximum pole modelled in the off-midplane expansion; 0 is dipole; 2 is quadrupole; etc */
+    size_t getMaxOrder() const {return maxOrder_m;}
 
-    /** Set the centre of the sector */
-    void setCentre(Vector_t centre) const {centre_m = centre;}
+    /** Set the maximum pole modelled in the off-midplane expansion; 0 is dipole; 2 is quadrupole; etc */
+    void setMaxOrder(size_t maxOrder) {maxOrder_m = maxOrder;}
+
+    /** Calculate the df coefficients, ready for field generation
+     *
+     *  Must be called following any update to the the field parameters, in
+     *  order for correct field to be calculated.
+     */
+    void calculateDfCoefficients();
+
+    /** Return the calculated df coefficients */
+    std::vector<std::vector<double> > getDfCoefficients() {return dfCoefficients_m;}
 
   private:
-    // get the fringe field model at azimuthal angle phi
-    double getFringe(double phi);
+    /** Copy constructor */
+    SpiralSector(const SpiralSector &right);
 
+    SpiralSector& operator=(const SpiralSector& rhs);
     PlanarArcGeometry planarArcGeometry_m;
     BMultipoleField dummy;
 
-    double theta1_m;
-    double theta2_m;
+    size_t maxOrder_m = 0;
+    double tanDelta_m = 0.;
     double k_m;
     double Bz_m;
     double r0_m;
     Vector_t centre_m;
     endfieldmodel::EndFieldModel* endField_m = NULL;
+    const double fp_tolerance = 1e-18;
+    std::vector<std::vector<double> > dfCoefficients_m;
 };
 
-
 #endif
+

tests/classic_src/AbsBeamline/SpiralSectorTest.cpp

+/*
+ *  Copyright (c) 2017, Chris Rogers
+ *  All rights reserved.
+ *  Redistribution and use in source and binary forms, with or without
+ *  modification, are permitted provided that the following conditions are met:
+ *  1. Redistributions of source code must retain the above copyright notice,
+ *     this list of conditions and the following disclaimer.
+ *  2. Redistributions in binary form must reproduce the above copyright notice,
+ *     this list of conditions and the following disclaimer in the documentation
+ *     and/or other materials provided with the distribution.
+ *  3. Neither the name of STFC nor the names of its contributors may be used to
+ *     endorse or promote products derived from this software without specific
+ *     prior written permission.
+ *
+ *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ *  POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include <cmath>
+#include <fstream>
+
+#include "gtest/gtest.h"
+
+
+#include "Classic/AbsBeamline/EndFieldModel/Tanh.h"
+#include "Classic/AbsBeamline/SpiralSector.h"
+
+class SpiralSectorTest : public ::testing::Test { 
+public: 
+    SpiralSectorTest() : sector_m(NULL), fout_m("spiralsectortest.out") { 
+    }
+ 
+    void SetUp( ) { 
+        sector_m = new SpiralSector("test");
+        // characteristic length is R*dphi => 0.6545 m
+        endfieldmodel::Tanh* tanh = new endfieldmodel::Tanh(M_PI/12., M_PI/96., 10);
+        sector_m->setEndField(tanh);
+        sector_m->setTanDelta(0.); //tan(M_PI/3.));
+        sector_m->setR0(10.);
+        sector_m->setDipoleConstant(1.);
+        sector_m->setFieldIndex(6);
+        sector_m->setMaxOrder(4);
+        sector_m->calculateDfCoefficients();
+    }
+ 
+    void TearDown( ) { 
+        delete sector_m;
+        sector_m = NULL;
+    }
+
+    ~SpiralSectorTest() {
+    }
+
+    SpiralSector* sector_m;
+    std::ofstream fout_m;
+
+    double get_div_B(Vector_t pos, Vector_t delta) {
+        Vector_t mom(0., 0., 0.);
+        Vector_t E(0., 0., 0.);
+        Vector_t B(0., 0., 0.);
+        double t = 0;
+        Vector_t div_elements(0., 0., 0.);
+        // dB_x/dx ~ (B_x(x+dx) - B_x(x-dx))/2/dx
+        for (size_t i = 0; i < 3; ++i) {
+            pos[i] += delta[i];
+            sector_m->apply(pos, mom, t, E, B);
+            div_elements[i] += B[i];
+            pos[i] -= 2.*delta[i];
+            sector_m->apply(pos, mom, t, E, B);
+            div_elements[i] -= B[i];
+            pos[i] += delta[i];
+            div_elements[i] /= delta[i]*2.;
+        }
+        std::cerr << "\n\n";
+        sector_m->apply(pos, mom, t, E, B);
+        std::cerr << "Div elements " << div_elements << std::endl;
+        return div_elements[0] + div_elements[1] + div_elements[2];
+    }
+
+    double get_div_B_cylindrical(Vector_t posCyl, Vector_t delta) {
+        Vector_t mom(0., 0., 0.);
+        Vector_t E(0., 0., 0.);
+        Vector_t B(0., 0., 0.);
+        double t = 0;
+        Vector_t div_elements(0., 0., 0.);
+        // dB_x/dx ~ (B_x(x+dx) - B_x(x-dx))/2/dx
+        for (size_t i = 0; i < 3; ++i) {
+            posCyl[i] += delta[i];
+            sector_m->applyCylindrical(posCyl, mom, t, E, B);
+            div_elements[i] += B[i];
+            posCyl[i] -= 2.*delta[i];
+            sector_m->applyCylindrical(posCyl, mom, t, E, B);
+            div_elements[i] -= B[i];
+            posCyl[i] += delta[i];
+            div_elements[i] /= delta[i]*2.;
+        }
+        sector_m->applyCylindrical(posCyl, mom, t, E, B);
+        div_elements[0] += B[0]/posCyl[0];
+        div_elements[1] /= posCyl[0];
+        return div_elements[0] + div_elements[1] + div_elements[2];
+    }
+
+private:
+
+};
+
+Vector_t get_curl_B(Vector_t pos, Vector_t delta) {
+    return Vector_t();
+}
+
+TEST_F(SpiralSectorTest, ApplyTest) {
+    std::vector<std::vector<double> > coeff = sector_m->getDfCoefficients();
+    std::cerr << "Coefficients" << std::endl;
+    for (size_t n = 0; n < coeff.size(); ++n) {
+        for (size_t i = 0; i < coeff[n].size(); ++i) {
+            std::cerr << coeff[n][i] << " ";
+        }
+        std::cerr << std::endl;
+    }
+    std::cerr << std::endl;
+    for (double y = 0.01; y < 0.055; y += 0.1 /*0.01*/) {
+        for (double r = 10.; r < 12.1; r += 10. /*0.5*/) {
+                for (double phi = M_PI/12./*-M_PI/3.*/; phi < M_PI/3.; phi += M_PI/*/600.*/) {
+                double x = r*cos(phi);
+                double z = r*sin(phi);
+                Vector_t pos(r, y, phi);
+                Vector_t posCart(x, y, z);
+                Vector_t mom(0., 0., 0.);
+                Vector_t E(0., 0., 0.);
+                Vector_t B(0., 0., 0.);
+                Vector_t BCart(0., 0., 0.);
+                double t = 0;
+                
+                sector_m->applyCylindrical(pos, mom, t, E, B);
+                //sector_m->apply(posCart, mom, t, E, BCart);
+                double div_B_cyl = 0.; //get_div_B_cylindrical(pos, Vector_t(1e-3, 1e-3, 1e-3));
+                double div_B_cart = get_div_B(posCart, Vector_t(1e-3, 1e-3, 1e-3));
+                fout_m << x << "    " << y << "    " << z << "    " << r << "    " << phi << "    " << B[0] << "    " << B[1] << "    "
+                       << B[2] << "    " << div_B_cyl << "    " << BCart[0] << "    " << BCart[1] << "    "
+                       << BCart[2] << "     " << div_B_cart << std::endl;
+            }
+        }
+        fout_m << std::endl;
+    }
+}
+
+