Commit e0545fcc authored by ext-rogers_c's avatar ext-rogers_c
Browse files

Add scaling FFAG magnet code

parents 6c350658 cefeb59b
......@@ -65,7 +65,7 @@ set (BOOSTROOT $ENV{BOOST_DIR})
set(Boost_USE_STATIC_LIBS ON)
set(Boost_USE_MULTITHREADED OFF)
set(Boost_USE_STATIC_RUNTIME OFF)
find_package (Boost 1.55.0 REQUIRED COMPONENTS filesystem system iostreams regex)
find_package (Boost 1.55.0 REQUIRED COMPONENTS filesystem system iostreams regex chrono)
  • @ext-rogers_c: why was a dependency to boost::chrono added? I don't see any mention elsewhere in the commit.

  • @snuverink you are right, but it was a missing pre-existing dependency. I noticed because it caused my build to fail. Specifically

    opt-pilot/Util/Trace/Timestamp.h

    requires it, so I added the dependency explicitly.

  • Hmm, strange, the dependency is checked in opt-pilot/CMakeLists.txt so I had thought you should have got it from there.

  • Oh so it is. Not sure why it didn't pick up, I probably screwed something up. I reverted CMakeLists.txt in OPAL anyway.

Please register or sign in to reply
if (Boost_INCLUDE_DIRS)
message (STATUS "Found boost include dir: ${Boost_INCLUDE_DIR}")
message (STATUS "Found boost library dir: ${Boost_LIBRARY_DIR}")
......@@ -250,4 +250,4 @@ install (
FILES ${CMAKE_CURRENT_BINARY_DIR}/${PROJECT_NAME}Config_install.cmake
DESTINATION "${CMAKE_INSTALL_PREFIX}/lib/cmake/${PROJECT_NAME}"
RENAME ${PROJECT_NAME}Config.cmake
)
\ No newline at end of file
)
......@@ -45,6 +45,7 @@
#include "AbsBeamline/RFQuadrupole.h"
#include "AbsBeamline/SBend.h"
#include "AbsBeamline/SBend3D.h"
#include "AbsBeamline/ScalingFFAGMagnet.h"
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
......@@ -786,6 +787,16 @@ void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) {
"Need to define a RINGDEFINITION to use SBend3D element");
}
void ParallelCyclotronTracker::visitScalingFFAGMagnet(const ScalingFFAGMagnet &bend) {
*gmsg << "Adding ScalingFFAGMagnet" << endl;
if (opalRing_m != NULL) {
opalRing_m->appendElement(bend);
} else {
throw OpalException("ParallelCyclotronTracker::visitScalingFFAGMagnet",
"Need to define a RINGDEFINITION to use ScalingFFAGMagnet element");
}
}
void ParallelCyclotronTracker::visitVariableRFCavity(const VariableRFCavity &cav) {
*gmsg << "Adding Variable RF Cavity" << endl;
if (opalRing_m != NULL)
......
......@@ -139,6 +139,9 @@ public:
/// Apply the algorithm to a SBend3D.
virtual void visitSBend3D(const SBend3D &);
/// Apply the algorithm to a ScalingFFAGMagnet.
virtual void visitScalingFFAGMagnet(const ScalingFFAGMagnet &bend);
/// Apply the algorithm to a Separator.
virtual void visitSeparator(const Separator &);
......
......@@ -44,25 +44,24 @@ std::string(" field file, for checking that fields are read in correctly")+
std::string(" from disk. The fields are written out on a Cartesian grid.");
DumpFields::DumpFields() :
Action(10, "DUMPFIELDS", dumpfields_docstring.c_str()),
grid_m(NULL), filename_m("") {
Action(10, "DUMPFIELDS", dumpfields_docstring.c_str()) {
// would be nice if "steps" could be integer
itsAttr[0] = Attributes::makeReal
("X_START", "Start point in the grid in x [mm]");
("X_START", "Start point in the grid in x [m]");
itsAttr[1] = Attributes::makeReal
("DX", "Grid step size in x [mm]");
("DX", "Grid step size in x [m]");
itsAttr[2] = Attributes::makeReal
("X_STEPS", "Number of steps in x");
itsAttr[3] = Attributes::makeReal
("Y_START", "Start point in the grid in y [mm]");
("Y_START", "Start point in the grid in y [m]");
itsAttr[4] = Attributes::makeReal
("DY", "Grid step size in y [mm]");
("DY", "Grid step size in y [m]");
itsAttr[5] = Attributes::makeReal
("Y_STEPS", "Number of steps in y");
itsAttr[6] = Attributes::makeReal
("Z_START", "Start point in the grid in z [mm]");
("Z_START", "Start point in the grid in z [m]");
itsAttr[7] = Attributes::makeReal
("DZ", "Grid step size in z [mm]");
("DZ", "Grid step size in z [m]");
itsAttr[8] = Attributes::makeReal
("Z_STEPS", "Number of steps in z");
itsAttr[9] = Attributes::makeString
......@@ -71,13 +70,17 @@ DumpFields::DumpFields() :
registerOwnership(AttributeHandler::STATEMENT);
}
DumpFields::DumpFields(const std::string &name, DumpFields *parent):
Action(name, parent)
{}
DumpFields::~DumpFields() {
delete grid_m;
dumpsSet_m.erase(this);
}
DumpFields* DumpFields::clone(const std::string &name) {
DumpFields* dumper = new DumpFields();
DumpFields* dumper = new DumpFields(name, this);
if (grid_m != NULL) {
dumper->grid_m = grid_m->clone();
}
......@@ -162,9 +165,9 @@ void DumpFields::writeFieldThis(Component* field) {
}
// set precision
fout << grid_m->end().toInteger() << "\n";
fout << 1 << " x [mm]\n";
fout << 2 << " y [mm]\n";
fout << 3 << " z [mm]\n";
fout << 1 << " x [m]\n";
fout << 2 << " y [m]\n";
fout << 3 << " z [m]\n";
fout << 4 << " Bx [kGauss]\n";
fout << 5 << " By [kGauss]\n";
fout << 6 << " Bz [kGauss]\n";
......@@ -184,4 +187,4 @@ void DumpFields::writeFieldThis(Component* field) {
"Something went wrong during writing "+filename_m);
}
fout.close();
}
\ No newline at end of file
}
......@@ -64,6 +64,9 @@ class DumpFields : public Action {
/** Constructor */
DumpFields();
/** Constructor */
DumpFields(const std::string &name, DumpFields *parent);
/** Destructor deletes grid_m and if in the dumps set, take it out */
virtual ~DumpFields();
......@@ -104,8 +107,8 @@ class DumpFields : public Action {
virtual void buildGrid();
static void checkInt(double value, std::string name, double tolerance = 1e-9);
interpolation::ThreeDGrid* grid_m;
std::string filename_m;
interpolation::ThreeDGrid* grid_m = NULL;
std::string filename_m = "";
static std::unordered_set<DumpFields*> dumpsSet_m;
static std::string dumpfields_docstring;
......
......@@ -56,6 +56,7 @@ class TravelingWave;
class RFQuadrupole;
class SBend;
class SBend3D;
class ScalingFFAGMagnet;
class Cyclotron;
class Separator;
class Septum;
......@@ -183,6 +184,9 @@ public:
/// Apply the algorithm to a solenoid.
virtual void visitSolenoid(const Solenoid &) = 0;
/// Apply the algorithm to a solenoid.
virtual void visitScalingFFAGMagnet(const ScalingFFAGMagnet &) = 0;
/// Apply the algorithm to a source.
virtual void visitSource(const Source &) = 0;
......@@ -240,4 +244,4 @@ void BeamlineVisitor::visitRBend3D(const RBend3D &) {
}
#endif // CLASSIC_BeamlineVisitor_HH
\ No newline at end of file
#endif // CLASSIC_BeamlineVisitor_HH
add_subdirectory(EndFieldModel)
set (_SRCS
AlignWrapper.cpp
AttributeSet.cpp
......@@ -31,6 +33,7 @@ set (_SRCS
Ring.cpp
SBend.cpp
SBend3D.cpp
ScalingFFAGMagnet.cpp
Separator.cpp
Septum.cpp
Solenoid.cpp
......@@ -79,6 +82,7 @@ set (HDRS
Ring.h
SBend3D.h
SBend.h
ScalingFFAGMagnet.h
SectorFieldMapComponent.h
Separator.h
Septum.h
......@@ -90,4 +94,4 @@ set (HDRS
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")
set (_SRCS
EndFieldModel.cpp
Enge.cpp
Tanh.cpp
)
include_directories (
${CMAKE_CURRENT_SOURCE_DIR}
)
ADD_OPAL_SOURCES(${_SRCS})
set (HDRS
EndFieldModel.h
Enge.h
Tanh.h
)
install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/AbsBeamline/EndFieldModel/")
/*
* 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 <algorithm>
#include "AbsBeamline/EndFieldModel/EndFieldModel.h"
namespace endfieldmodel {
bool GreaterThan(std::vector<int> v1, std::vector<int> v2) {
size_t n1(v1.size()), n2(v2.size());
for (size_t i = 0; i < n1 && i < n2; ++i) {
if (v1[n1-1-i] > v2[n2-1-i]) return true;
if (v1[n1-1-i] < v2[n2-1-i]) return false;
}
return false;
}
std::vector< std::vector<int> > CompactVector(
std::vector< std::vector<int> > vec) {
// first sort the list
std::sort(vec.begin(), vec.end(), GreaterThan);
// now look for n = n+1
for (size_t j = 0; j < vec.size()-1; ++j) {
while (j < vec.size()-1 && IterableEquality(
vec[j].begin()+1, vec[j].end(),
vec[j+1].begin()+1, vec[j+1].end()) ) {
vec[j][0] += vec[j+1][0];
vec.erase(vec.begin()+j+1);
}
}
return vec;
}
}
/*
* 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.
*/
#ifndef ENDFIELDMODEL_ENDFIELDMODEL_H_
#define ENDFIELDMODEL_ENDFIELDMODEL_H_
#include <iostream>
#include <vector>
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;
};
std::vector< std::vector<int> > CompactVector
(std::vector< std::vector<int> > vec);
/// CompactVector helper function, used for sorting
bool GreaterThan(std::vector<int> v1, std::vector<int> v2);
/** Return a == b if a and b are same size and a[i] == b[i] for all i.
*
* The following operations must be defined for TEMP_ITER it:
* - ++it prefix increment operator
* - (*it) (that is unary *, i.e. dereference operator)
* - it1 != it2 not equals operator
* - (*it1) != (*it2) not equals operator of dereferenced object
*
* Call like e.g. \n
* std::vector<int> a,b;\n
* bool test_equal = IterableEquality(a.begin(), a.end(), b.begin(),
* b.end());\n
*
* Can give a segmentation fault if a.begin() is not between a.begin() and
* a.end() (inclusive)
*/
template <class TEMP_ITER>
bool IterableEquality(TEMP_ITER a_begin, TEMP_ITER a_end,
TEMP_ITER b_begin, TEMP_ITER b_end);
/** Return a == b if a and b are same size and a[i] == b[i] for all i.
*
* The following operations must be defined for TEMP_ITER it:
* - ++it prefix increment operator
* - (*it) (that is unary *, i.e. dereference operator)
* - it1 != it2 not equals operator
* - (*it1) != (*it2) not equals operator of dereferenced object
*
* Call like e.g. \n
* std::vector<int> a,b;\n
* bool test_equal = IterableEquality(a.begin(), a.end(), b.begin(),
* b.end());\n
*
* Can give a segmentation fault if a.begin() is not between a.begin() and
* a.end() (inclusive)
*/
template <class TEMP_ITER>
bool IterableEquality(TEMP_ITER a_begin, TEMP_ITER a_end,
TEMP_ITER b_begin, TEMP_ITER b_end);
template <class TEMP_CLASS>
bool IterableEquality(const TEMP_CLASS& a, const TEMP_CLASS& b) {
return IterableEquality(a.begin(), a.end(), b.begin(), b.end());
}
template <class TEMP_ITER>
bool IterableEquality(TEMP_ITER a_begin, TEMP_ITER a_end, TEMP_ITER b_begin,
TEMP_ITER b_end) {
TEMP_ITER a_it = a_begin;
TEMP_ITER b_it = b_begin;
while (a_it != a_end && b_it != b_end) {
if (*a_it != *b_it) return false;
++a_it;
++b_it;
}
if ( a_it != a_end || b_it != b_end ) return false;
return true;
}
}
#endif
/*
* 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 <math.h>
#include "gsl/gsl_sf_gamma.h"
#include "gsl/gsl_sf_pow_int.h"
#include "AbsBeamline/EndFieldModel/Enge.h"
namespace endfieldmodel {
// Use
// d^n E/dx^n = a_n1m1 F(n1) g(m1) + a_n2m1m2 F(n2) g(m1)g(m2)+...
// where
double Enge::GetEnge(double x, int n) const {
std::vector< std::vector<int> > qt = GetQIndex(n);
std::vector<double> g;
double e(0.);
for (size_t i = 0; i < qt.size(); ++i) {
double ei(qt[i][0]);
for (size_t j = 1; j < qt[i].size(); ++j) {
if (j > g.size()) g.push_back(GN(x, j-1));
ei *= gsl_sf_pow_int(g[j-1], qt[i][j]);
}
if (ei != ei) ei = 0; // div 0, usually g == 0 and index < 0
e += ei;
}
return e;
}
// h = a_0+a_1 (x/w)+a_2 (x/w)^2+a_3 (x/w)^3+...+a_m (x/w)^m
// h^(n) = d^nh/dx^n = sum^m_{i=n} a_i x^{i-n}/w^i i!/n!
double Enge::HN(double x, int n) const {
double hn = 0;
// optimise by precalculating factor
for (unsigned int i = n; i <_a.size(); i++)
hn += _a[i]/gsl_sf_pow_int(_lambda, i)*gsl_sf_pow_int(x, i-n)*
static_cast<double>(gsl_sf_fact(i))
/static_cast<double>(gsl_sf_fact(i-n));
return hn;
}
// g = 1+exp(h)
// g^(n) = d^ng/dx^n
double Enge::GN(double x, int n) const {
if (n == 0) return 1+exp(HN(x, 0)); // special case
std::vector<double> hn(n+1);
for (int i = 0; i <= n; i++) hn[i] = HN(x, i);
double exp_h0 = exp(hn[0]);
double gn = 0;
for (size_t i = 0; i < _h[n].size(); ++i) {
double gnj = _h[n][i][0]*exp_h0;
for (size_t j = 1; j < _h[n][i].size(); ++j)
gnj *= gsl_sf_pow_int(hn[j], _h[n][i][j]);
gn += gnj;
}
return gn;
}
// _q[i][j][k]; urk, 3d vector
// i indexes the derivative of f;
// j indexes the element in f derivative
// k indexes the derivative of g
// this will quickly become grotesque
std::vector< std::vector< std::vector<int> > > Enge::_q;
std::vector< std::vector< std::vector<int> > > Enge::_h;
void Enge::SetEngeDiffIndices(size_t n) {
if (_q.size() == 0) {
_q.push_back(std::vector< std::vector<int> >(1, std::vector<int>(3)) );
_q[0][0][0] = +1; // f_0 = 1*g^(-1)
_q[0][0][1] = -1;
_q[0][0][2] = 0;
}
for (size_t i = _q.size(); i < n+1; ++i) {
_q.push_back(std::vector< std::vector<int> >() );
for (size_t j = 0; j < _q[i-1].size(); ++j) {
size_t k_max = _q[i-1][j].size();
std::vector<int> new_vec(_q[i-1][j]);
// derivative of g^-n0 = -n0*g^(-n0-1)*g(1)
new_vec[0] *= new_vec[1]; // alpha *= g(0) power
new_vec[1] -= 1; // g(0) power -= 1
new_vec[2] += 1; // g(1) power += 1
_q[i].push_back(new_vec);
for (size_t k = 2; k < k_max; ++k) { // 0 is alpha; 1 is g(0)
// derivative of g(k)^nk = nk g(k+1) g(k)^(nk-1)
if (_q[i-1][j][k] > 0) {
std::vector<int> new_vec(_q[i-1][j]);
if ( k == k_max-1 ) new_vec.push_back(0); // need enough coefficients
new_vec[0] *= new_vec[k];
new_vec[k] -= 1;
new_vec[k+1] += 1;
_q[i].push_back(new_vec);
}
}
}
}
if (_h.size() == 0) {
// first one is special case (1+e^h dealt with explicitly)
_h.push_back(std::vector< std::vector<int> >());
// second is (1*e^h h'^1)
_h.push_back(std::vector< std::vector<int> >());
_h[1].push_back(std::vector<int>(2, 1));
}
for (size_t i = _h.size(); i < n+1; ++i) {
_h.push_back(std::vector< std::vector<int> >());
for (size_t j = 0; j < _h[i-1].size(); ++j) {
// d/dx k0 e^g g(1)^k1 ... g(n)^kn ... = k0 e^g g(1)^(k1+1) ... g(n)^kn
// + SUM_n k0 kn e^g ... g(n)^(kn-1) g(n-1)
std::vector<int> new_vec(_h[i-1][j]);
new_vec[1] += 1;
_h[i].push_back(new_vec);
for (size_t k = 1; k <_h[i-1][j].size(); ++k) {
if (_h[i-1][j][k] > 0) {
std::vector<int> new_vec(_h[i-1][j]);
if (k == _h[i-1][j].size()-1) new_vec.push_back(0);
new_vec[0] *= new_vec[k];
new_vec[k] -= 1;
new_vec[k+1] += 1;
_h[i].push_back(new_vec);
}
}
}
_h[i] = CompactVector(_h[i]);
}
}
}
/*
* 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.
*/
#ifndef ENDFIELDMODEL_ENGE_H_
#define ENDFIELDMODEL_ENGE_H_
#include <iostream>
#include <vector>
#include "AbsBeamline/EndFieldModel/EndFieldModel.h"
namespace endfieldmodel {
class Enge : public EndFieldModel {
public:
/** Default constructor */
Enge() : _a(), _lambda(0.) {SetEngeDiffIndices(10);}
/** Builds Enge function with parameters a_0, a_1, ..., lambda and x0.
*
* max_index is the maximum derivative that will be used in calcula