src/Classic/AbsBeamline/EndFieldModel/Tanh.cpp

@@ -30,6 +30,14 @@
 #include "gsl/gsl_sf_gamma.h"
 #include "gsl/gsl_sf_pow_int.h"
 
+#include <gsl/gsl_integration.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_sf_gamma.h>
+
+#include "AbsBeamline/MultipoleTFunctions/tanhDeriv.h"
 #include "AbsBeamline/EndFieldModel/Tanh.h"
 
 namespace endfieldmodel {
@@ -46,7 +54,6 @@ Tanh* 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;
@@ -63,14 +70,22 @@ double Tanh::getNegTanh(double x, int n) const {
   double t = 0;
   double lam_n = gsl_sf_pow_int(_lambda, n);
   double tanh_x = tanh((x-_x0)/_lambda);
-  for (size_t i = 0; i < _tdi[n].size(); i++)
+  for (size_t i = 0; i < _tdi[n].size(); i++) {
     t += 1./lam_n*static_cast<double>(_tdi[n][i][0])
             *gsl_sf_pow_int(tanh_x, _tdi[n][i][1]);
+  }
   return t;
 }
 
 double Tanh::function(double x, int n) const {
-  return (getTanh(x, n)-getNegTanh(x, n))/2.;
+  if (n > 10) {
+        // BUG this does not work and I do not understand for why not
+        double t = tanhderiv::integrate(x, _x0, _lambda, _lambda, n);
+        return t;
+  }
+  double tplus = getTanh(x, n);
+  double tminus = getNegTanh(x, n);
+  return (tplus-tminus)/2.;
 }
 
 void Tanh::setMaximumDerivative(size_t n) {
@@ -114,8 +129,6 @@ void Tanh::rescale(double scaleFactor) {
   _lambda *= scaleFactor;
 }
 
-
-
 std::ostream& Tanh::print(std::ostream& out) const {
     out << "Tanh model with centre length: " << _x0 
         << " end length: " << _lambda;

src/Classic/AbsBeamline/MultipoleT.cpp

@@ -72,6 +72,7 @@ MultipoleT::MultipoleT(const MultipoleT &right):
     rotation_m(right.rotation_m),
     variableRadius_m(right.variableRadius_m),
     boundingBoxLength_m(right.boundingBoxLength_m),
+    magnetStart_m(right.magnetStart_m),
     verticalApert_m(right.verticalApert_m),
     horizApert_m(right.horizApert_m),
     dummy() {
@@ -192,11 +193,11 @@ Vector_t MultipoleT::transformCoords(const Vector_t &R) {
         if (alpha != 0.0 && angle_m != 0.0) {
             X[0] = R[2] / sin(alpha) - radius;
             X[1] = R[1];
-            X[2] = radius * alpha;// + boundingBoxLength_m;
+            X[2] = radius * alpha - fringeField_l.Tanh::getX0()-magnetStart_m;
         } else {
             X[0] = R[0];
             X[1] = R[1];
-            X[2] = R[2];// + boundingBoxLength_m;
+            X[2] = R[2] - fringeField_l.Tanh::getX0()-magnetStart_m;
         }
     } else {
         // If radius is variable
@@ -208,7 +209,7 @@ Vector_t MultipoleT::transformCoords(const Vector_t &R) {
         std::vector<double> r = t.getTransformation();
         X[0] = r[0];
         X[1] = r[1];
-        X[2] = r[2];
+        X[2] = r[2] - fringeField_l.Tanh::getX0()-magnetStart_m;
     }
     return X;
 }
@@ -536,6 +537,16 @@ double MultipoleT::getFn(std::size_t n, double x, double s) {
 void MultipoleT::initialise() {
     planarArcGeometry_m.setElementLength(length_m);
     planarArcGeometry_m.setCurvature(angle_m / length_m);
+    if (verticalApert_m <= 0.0) {
+        throw OpalException("MultipoleT::initialise()",
+            "MultipoleT "+getName()+" had vertical aperture <= 0");
+    } else if (horizApert_m <= 0.0) {
+        throw OpalException("MultipoleT::initialise()",
+            "MultipoleT "+getName()+" had horizontal aperture <= 0");
+    } else if (boundingBoxLength_m <= 0.0) {
+        throw OpalException("MultipoleT::initialise()",
+            "MultipoleT "+getName()+" had bounding box <= 0");
+    }
 }
 
 void MultipoleT::initialise(PartBunchBase<double, 3>* bunch,

src/Classic/AbsBeamline/MultipoleT.h

 /*
  *  Copyright (c) 2017, Titus Dascalu
  *  Copyright (c) 2018, Martin Duy Tat
+ *  Copyright (c) 2023, 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:
@@ -245,6 +246,14 @@ public:
      */
     void setBoundingBoxLength(const double &boundingBoxLength);
 
+    /** Get distance to entrance of magnet */
+    double getMagnetStart() const;
+    /** Set distance to magnet start from placement pointer
+     *  \param magnetStart -> Distance to magnet entrance
+     */
+    void setMagnetStart(const double &magnetStart);
+
+
 private:
     MultipoleT operator=(const MultipoleT& rhs);
     // End fields
@@ -285,8 +294,12 @@ private:
     double rotation_m;
     /** Variable radius flag */
     bool variableRadius_m;
-    /** Distance between centre of magnet and entrance */
+    /** Full length of the bounding box */
     double boundingBoxLength_m;
+    /** Distance from the placement pointer to the start of the field */
+    double magnetStart_m = 0.0;
+    /** Distance from the placement pointer to the next placement pointer */
+    double magnetEnd_m = 0.0;
     /** Get field component methods
      */
     double getBx (const Vector_t &R);
@@ -434,5 +447,14 @@ inline
     void MultipoleT::setBoundingBoxLength(const double &boundingBoxLength) {
         boundingBoxLength_m = boundingBoxLength;
 }
+inline
+    double MultipoleT::getMagnetStart() const {
+        return magnetStart_m;
+}
+
+inline
+    void MultipoleT::setMagnetStart(const double &magnetStart) {
+        magnetStart_m = magnetStart;
+}
 
 #endif

src/Elements/OpalMultipoleT.cpp

@@ -35,7 +35,9 @@ OpalMultipoleT::OpalMultipoleT():
     OpalElement(SIZE, "MULTIPOLET",
                 "The \"MULTIPOLET\" element defines a combined function multipole.") {
     itsAttr[TP] = Attributes::makeRealArray
-        ("TP", "Transverse Profile derivatives in T m^(-k)");
+        ("TP",
+        "Array of multipolar field strengths b_k. The field generated in the "
+        "flat top is B = b_k x^k [T m^(-k)]");
 
     itsAttr[LFRINGE] = Attributes::makeReal
         ("LFRINGE", "The length of the left end field [m]");
@@ -70,11 +72,13 @@ OpalMultipoleT::OpalMultipoleT():
     itsAttr[BBLENGTH] = Attributes::makeReal
         ("BBLENGTH", "Distance between centre of magnet and entrance [m]");
 
+    itsAttr[MAGNET_START] = Attributes::makeReal
+        ("MAGNET_START", "Distance between the placement pointer and the magnet placement [m]");
+
     registerOwnership();
     setElement(new MultipoleT("MULTIPOLET"));
 }
 
-
 OpalMultipoleT::OpalMultipoleT(const std::string& name,
                                OpalMultipoleT* parent):
     OpalElement(name, parent) {
@@ -127,6 +131,7 @@ void OpalMultipoleT::update() {
     multT->setMaxXOrder(Attributes::getReal(itsAttr[MAXXORDER]));
     multT->setRotation(Attributes::getReal(itsAttr[ROTATION]));
     multT->setEntranceAngle(Attributes::getReal(itsAttr[EANGLE]));
+    multT->setMagnetStart(Attributes::getReal(itsAttr[MAGNET_START]));
 
     for (int comp = 0; comp < transSize; comp++) {
         double units = Units::T2kG * gsl_sf_pow_int(1.0, comp); // T m^-comp -> kG mm^-comp

src/Elements/OpalMultipoleT.h

@@ -37,7 +37,8 @@ public:
         ROTATION,        // Rotation angle about central axis for skew elements
         EANGLE,          // Entrance angle
         VARRADIUS,       // Variable radius flag
-        BBLENGTH,        // Distance between centre of magnet and entrance
+        BBLENGTH,        // Length within which field is calculated
+        MAGNET_START,    // Distance to magnet entrance
         SIZE             // size of the enum
     };
 

src/Elements/OpalScalingFFAMagnet.cpp

@@ -44,8 +44,8 @@ OpalScalingFFAMagnet::OpalScalingFFAMagnet() :
 
     itsAttr[END_FIELD_MODEL] = Attributes::makeString
         ("END_FIELD_MODEL",
-         "Names the end field model of the ring, giving dipole field along a line of "
-         "constant radius. If blank, uses the default 'END_LENGTH' and 'CENTRE_LENGTH' "
+         "Names the end field model of the magnet, giving the field magnitude along a line of "
+         "constant radius. If blank, uses the 'END_LENGTH' and 'CENTRE_LENGTH' "
          "parameters and a tanh model. If 'END_FIELD_MODEL' is not blank, Opal will seek "
          "an END_FIELD_MODEL corresponding to the name defined in this string.");
 

src/PyOpal/PyElements/PyMultipoleT.cpp

@@ -40,6 +40,7 @@ std::vector<PyOpalObjectNS::AttributeDef> PyOpalObjectNS::PyOpalObject<OpalMulti
     {"ROTATION", "rotation", "", PyOpalObjectNS::DOUBLE},
     {"VARRADIUS", "variable_radius", "", PyOpalObjectNS::BOOL},
     {"BBLENGTH", "bounding_box_length", "", PyOpalObjectNS::DOUBLE},
+    {"MAGNET_START", "magnet_start", "", PyOpalObjectNS::DOUBLE},
     {"L", "length", "", PyOpalObjectNS::DOUBLE},
     {"DELETEONTRANSVERSEEXIT", "delete_on_transverse_exit", "", PyOpalObjectNS::BOOL}
 };

src/PyOpal/PyPython/ffa_field_mapper.py

@@ -34,6 +34,10 @@ class FFAFieldMapper():
     Class to make field maps, intended for FFAs/ring geometries
     """
     def __init__(self):
+        # arbitrary height and time offsets for map generation
+        self.z_position = 0.0
+        self.time = 0.0
+
         # for cylindrical field map
         self.r_points = []
         self.phi_points = []
@@ -134,8 +138,8 @@ class FFAFieldMapper():
                 phi_grid.append(phi)
                 point = (radius*math.cos(math.radians(phi)),
                          radius*math.sin(math.radians(phi)),
-                         0,
-                         0)
+                         self.z_position,
+                         self.time)
                 value = pyopal.objects.field.get_field_value(*point)
                 bz_grid.append(value[3])
                 if self.verbose > 0:
@@ -211,7 +215,7 @@ class FFAFieldMapper():
         axes.set_xlim(xlim)
         axes.set_ylim(ylim)
 
-    def field_map_cartesian(self, axes = None):
+    def field_map_cartesian(self, axes = None, variable = "bz"):
         """
         Plot a field map in cartesian coordinates.
 
@@ -228,9 +232,9 @@ class FFAFieldMapper():
             for pos_y in self.y_points:
                 x_grid.append(pos_x)
                 y_grid.append(pos_y)
-                point = (pos_x, pos_y, 0, 0)
-                value = pyopal.objects.field.get_field_value(*point)
-                bz_grid.append(value[3])
+                point = (pos_x, pos_y, self.z_position, self.time)
+                value = self.get_value(variable, [point])[0]
+                bz_grid.append(value)
                 if self.verbose > 0:
                     print("Field value at point", point,
                           "is B:", value[1:4], "E:", value[4:])
@@ -244,9 +248,9 @@ class FFAFieldMapper():
                     cmin=min_by, cmax=max_by, cmap=self.cmap, vmin=-cmax, vmax=cmax)
         axes.set_xlabel("x [m]")
         axes.set_ylabel("y [m]")
-        axes.set_title("$B_{z}$ [T]")
+        axes.set_title(self.axis_labels[variable])
         figure.colorbar(hist[3])
-        fig_fname = os.path.join(self.plot_dir, "scaling_ffa_map_cart.png")
+        fig_fname = os.path.join(self.plot_dir, f"scaling_ffa_map_cart_{variable}.png")
         figure.savefig(fig_fname)
         print("Generated cartesian field map in", fig_fname)
         return figure
@@ -290,8 +294,8 @@ class FFAFieldMapper():
         for phi in self.phi_points:
             point = (radius*math.cos(math.radians(phi)),
                      radius*math.sin(math.radians(phi)),
-                     0,
-                     0)
+                     self.z_position,
+                     self.time)
             value = pyopal.objects.field.get_field_value(*point)
             bz_points.append(value[3])
 
@@ -338,6 +342,7 @@ class FFAFieldMapper():
         - var2: variable in the denominator (range 1 to 4)
         """
         pos_vec = [pos_x, pos_y, pos_z, time]
+        var1 = self.field_variables.index(var1)
         var2 = self.position_variables.index(var2)
         pos_vec[var2] += self.delta_x
         field_plus = pyopal.objects.field.get_field_value(*pos_vec)[var1]
@@ -378,6 +383,22 @@ class FFAFieldMapper():
         ]
         return curl_b
 
+    def get_value(self, variable, point_list):
+        value_list = []
+        if variable in self.field_variables:
+            index = self.field_variables.index(variable)
+            for point in point_list:
+                value_list.append(pyopal.objects.field.get_field_value(*point)[index])
+        elif variable in self.deriv_variables:
+            if variable == "div_b":
+                for point in point_list:
+                    value_list.append(self.get_div_b(point[0], point[1], point[2], point[3]))
+            if variable == "curl_b":
+                for point in point_list:
+                    curl_b = self.get_curl_b(point[0], point[1], point[2], point[3])
+                    value_list.append((curl_b[0]**2+curl_b[1]**2+curl_b[2]**2)**0.5)
+        return value_list
+
     default_radial_contour = {"radius":0.0,
                               "linestyle":"-",
                               "colour":"grey",
@@ -391,3 +412,7 @@ class FFAFieldMapper():
     }
     position_variables = ["x", "y", "z", "t"]
     field_variables = ["out_of_bounds", "bx", "by", "bz", "ex", "ey", "ez"]
+    deriv_variables = ["div_b", "curl_b"]
+    axis_labels = {"bx":"B$_{x}$ [T]", "by":"B$_{y}$ [T]", "bz":"B$_{z}$ [T]",
+                   "ex":"E$_{x}$ []", "ey":"E$_{y}$ []", "ez":"E$_{z}$ []",
+                   "div_b":"$\\nabla .$ B [T/m]", "curl_b":"$|\\nabla \\times B|$ [T/m]"}

tests/opal_src/PyOpal/PyElements/test_multipolet.py

@@ -11,11 +11,21 @@
 # along with OPAL.  If not, see <https://www.gnu.org/licenses/>.
 
 import math
-import unittest 
+import unittest
+import sys
+import os
+
+try:
+    import matplotlib
+    import matplotlib.pyplot
+    import pyopal.objects.ffa_field_mapper
+except ImportError:
+    pass
 
 import pyopal.objects.field
 import pyopal.objects.minimal_runner
 import pyopal.elements.multipolet
+import pyopal.objects.field
 
 class TestMultipoleTRunner(pyopal.objects.minimal_runner.MinimalRunner):
     """Test runner - I wanted to check the placement was okay"""
@@ -34,20 +44,23 @@ class TestMultipoleTRunner(pyopal.objects.minimal_runner.MinimalRunner):
     def make_multipolet(self, bz):
         """Make a default straight multipole"""
         multipole = pyopal.elements.multipolet.MultipoleT()
-        multipole.t_p = [bz] # dipole
-        multipole.left_fringe = 0.01
-        multipole.right_fringe = 0.01
-        multipole.length = self.length
-        multipole.horizontal_aperture = 0.1
-        multipole.vertical_aperture = 0.1
-        multipole.maximum_f_order = 1
-        multipole.entrance_angle = 0.0
-        multipole.maximum_x_order = 1
-        multipole.variable_radius = 0
-        multipole.rotation = 0
-        multipole.angle = 0.0
-        multipole.bounding_box_length = self.bb_length
-        multipole.delete_on_transverse_exit = False
+        multipole.set_attributes(
+            t_p = [bz, 0.0, 0.0, 0.0/(0.05**3)], # dipole
+            left_fringe = 0.01,
+            right_fringe = 0.01,
+            length = self.length,
+            horizontal_aperture = 0.1,
+            vertical_aperture = 0.1,
+            maximum_f_order = 6,
+            entrance_angle = 0.0,
+            maximum_x_order = 6,
+            variable_radius = 0,
+            rotation = 0,
+            angle = 0.0,
+            bounding_box_length = self.bb_length,
+            delete_on_transverse_exit = False,
+            magnet_start = 0.1,
+        )
         self.multipole = multipole
         return multipole
 
@@ -61,6 +74,7 @@ class TestMultipoleTRunner(pyopal.objects.minimal_runner.MinimalRunner):
 
     def postprocess(self):
         """Test placement - note that placement is from centre for this element"""
+        self.plot()
         self.test_field(-0.301, 0.0)
         self.test_field(-0.299, 8.317278708416809e-05)
         self.test_field(0.299, 8.317278708416809e-05)
@@ -68,6 +82,47 @@ class TestMultipoleTRunner(pyopal.objects.minimal_runner.MinimalRunner):
         self.test_field(-0.25, 0.75)
         self.test_field(0.25, 0.75)
 
+    def plot(self):
+        mapper = pyopal.objects.ffa_field_mapper.FFAFieldMapper()
+        mapper.x_points = [i*0.004+1.9 for i in range(50)]
+        mapper.y_points = [i*0.01 for i in range(80)]
+        mapper.z_point = 0.02
+        mapper.field_map_cartesian(variable = "bz")
+        mapper.field_map_cartesian(variable = "div_b")
+        mapper.field_map_cartesian(variable = "curl_b")
+
+    def _plot(self):
+        """Plot the field"""
+        if 'matplotlib' not in sys.modules:
+            return
+        x0 = 1.5
+        x1 = 2.5
+        y0 = 0.0
+        y1 = 1.0
+        n_x = 1000
+        n_y = 1000
+        x_list = []
+        y_list = []
+        bz_list = []
+        dx = (x1-x0)/n_x
+        dy = (y1-y0)/n_y
+        figure = matplotlib.pyplot.figure()
+        axes = figure.add_subplot(1, 1, 1)
+        for i in range(n_x-1):
+            for j in range(n_y-1):
+                x = x0+(i+0.5)*dx
+                y = y0+(j+0.5)*dy
+                field = pyopal.objects.field.get_field_value(x, y, 0, 0)
+                x_list.append(x)
+                y_list.append(y)
+                bz_list.append(field[3])
+        bin_x = [x0+(x1-x0)*i/n_x for i in range(n_x)]
+        bin_y = [y0+(y1-y0)*j/n_y for j in range(n_y)]
+        axes.hist2d(x_list, y_list, weights=bz_list, bins=[bin_x, bin_y])
+        fname = "test_multipolet_1.png"
+        figure.savefig(fname)
+        print(f"Saved figure in {os.getcwd()}/{fname}")
+
 class TestMultipoleT(unittest.TestCase):
     """
     Simple test on multipolet class
@@ -141,6 +196,7 @@ class TestMultipoleT(unittest.TestCase):
     def test_placement(self):
         """Check placement is okay"""
         runner = TestMultipoleTRunner()
+        runner.tmp_dir = "tmp"
         runner.execute_fork()
 
 if __name__ == "__main__":

tests/opal_src/PyOpal/PyObjects/test_track_run_scaling_ffa.py

@@ -53,7 +53,7 @@ class ScalingFFARunner(pyopal.objects.minimal_runner.MinimalRunner):
         self.probe_id = 1
         self.n_cells = 16
         self.bend_direction = 1 # set to -1 for ring centred on x,y = 8,0 m
-        self.dr = 1.0
+        self.dr = 0.5
         self.spiral_angle = math.pi/6
         ke = 3e-3 # 3 MeV
         self.momentum = ((ke+self.mass)**2-self.mass**2)**0.5
@@ -202,19 +202,20 @@ class ScalingFFARunner(pyopal.objects.minimal_runner.MinimalRunner):
         ]
 
 
-        phi_start = self.f_start/self.r0
-        phi_middle = (self.f_start+self.f_centre_length/2)/self.r0
-        phi_fringe_end = (self.f_start+self.f_centre_length)/self.r0
-        phi_end = self.f_end/self.r0
+        phi_start = math.degrees(self.f_start/self.r0)
+        phi_middle = math.degrees((self.f_start+self.f_centre_length/2)/self.r0)
+        phi_fringe_end = math.degrees((self.f_start+self.f_centre_length)/self.r0)
+        phi_end = math.degrees(self.f_end/self.r0)
+        spiral_deg = math.degrees(self.spiral_angle)
         mapper.spiral_contours = [
-            {"phi0":phi_start, "r0":self.r0, "spiral_angle":self.spiral_angle, "linestyle":"dashed", "colour":"blue", "label":"magnet_start"},
-            {"phi0":phi_middle, "r0":self.r0, "spiral_angle":self.spiral_angle, "linestyle":"dashed", "colour":"grey", "label":"magnet_start+\ncentre_length/2"},
-            {"phi0":phi_fringe_end, "r0":self.r0, "spiral_angle":self.spiral_angle, "linestyle":"dashed", "colour":"grey", "label":"magnet_start+\ncentre_length"},
-            {"phi0":phi_end, "r0":self.r0, "spiral_angle":self.spiral_angle, "linestyle":"dashed", "colour":"blue", "label":"magnet_end"},
+            {"phi0":phi_start, "r0":self.r0, "spiral_angle":spiral_deg, "linestyle":"dashed", "colour":"blue", "label":"magnet_start/r0"},
+            {"phi0":phi_middle, "r0":self.r0, "spiral_angle":spiral_deg, "linestyle":"dashed", "colour":"grey", "label":"(magnet_start+\ncentre_length/2)/r0"},
+            {"phi0":phi_fringe_end, "r0":self.r0, "spiral_angle":spiral_deg, "linestyle":"dashed", "colour":"grey", "label":"(magnet_start+\ncentre_length)/r0"},
+            {"phi0":phi_end, "r0":self.r0, "spiral_angle":spiral_deg, "linestyle":"dashed", "colour":"blue", "label":"magnet_end/r0"},
         ]
 
         mapper.r_points = [self.r0+i*0.001 for i in range(-1200, 1200+1, 10)]
-        mapper.phi_points = [i/8.0 for i in range(0, 180+1, 5)]
+        mapper.phi_points = [i/8.0/10 for i in range(0, 1800+1, 5)]
         mapper.field_map_cylindrical()
 
         mapper.x_points = [i*0.05 for i in range(-100, 100+1, 5)]