Commit 67bd3e4d authored by ext-rogers_c's avatar ext-rogers_c
Browse files

Implement ScalingFFAGMagnet as a element for tracking

parent 7bbcdd0f
......@@ -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"
......@@ -751,6 +752,15 @@ 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)
......@@ -2219,7 +2229,6 @@ bool ParallelCyclotronTracker::rk4(double x[], const double &t, const double &ta
double xtemp[PSdim];
// Evaluate f1 = f(x,t).
bool outOfBound = derivate(x, t, deriv1 , Pindex);
if (outOfBound) {
return false;
......
......@@ -126,6 +126,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,7 +56,7 @@ class TravelingWave;
class RFQuadrupole;
class SBend;
class SBend3D;
class SpiralSector;
class ScalingFFAGMagnet;
class Cyclotron;
class Separator;
class Septum;
......@@ -185,7 +185,7 @@ public:
virtual void visitSolenoid(const Solenoid &) = 0;
/// Apply the algorithm to a solenoid.
virtual void visitSpiralSector(const SpiralSector &) = 0;
virtual void visitScalingFFAGMagnet(const ScalingFFAGMagnet &) = 0;
/// Apply the algorithm to a source.
virtual void visitSource(const Source &) = 0;
......
......@@ -33,11 +33,11 @@ set (_SRCS
Ring.cpp
SBend.cpp
SBend3D.cpp
ScalingFFAGMagnet.cpp
Separator.cpp
Septum.cpp
Solenoid.cpp
Source.cpp
SpiralSector.cpp
Stripper.cpp
TravelingWave.cpp
VariableRFCavity.cpp
......@@ -82,13 +82,13 @@ set (HDRS
Ring.h
SBend3D.h
SBend.h
ScalingFFAGMagnet.h
SectorFieldMapComponent.h
Separator.h
Septum.h
Solenoid.h
Source.h
SpecificElementVisitor.h
SpiralSector.h
Stripper.h
TravelingWave.h
VariableRFCavity.h
......
......@@ -95,7 +95,8 @@ class Offset : public Component {
* - theta_in angle between the previous element and the displacement
* vector
* - theta_out angle between the displacement vector and the next element
* - displacement length of the displacement vector
* - displacement length of the displacement vector in the theta_in
* direction
*/
static Offset localCylindricalOffset(std::string name,
double theta_in,
......
......@@ -127,7 +127,6 @@ bool Ring::apply(const Vector_t &R, const Vector_t &P,
B += (scale_m * B_temp);
E += (scale_m * E_temp);
}
// std::cerr << "Ring::apply " << sections.size() << " Pos: " << R << " B: " << B << std::endl;
return outOfBounds;
}
......
......@@ -27,127 +27,121 @@
#include <cmath>
#include "AbsBeamline/SpiralSector.h"
#include "AbsBeamline/ScalingFFAGMagnet.h"
#include "Algorithms/PartBunch.h"
#include "AbsBeamline/BeamlineVisitor.h"
SpiralSector::SpiralSector(const std::string &name)
ScalingFFAGMagnet::ScalingFFAGMagnet(const std::string &name)
: Component(name),
planarArcGeometry_m(1., 1.), dummy() {
setElType(isDrift);
}
SpiralSector::SpiralSector(const SpiralSector &right)
ScalingFFAGMagnet::ScalingFFAGMagnet(const ScalingFFAGMagnet &right)
: Component(right),
planarArcGeometry_m(right.planarArcGeometry_m), dummy() {
planarArcGeometry_m(right.planarArcGeometry_m),
dummy(), maxOrder_m(right.maxOrder_m), tanDelta_m(right.tanDelta_m),
k_m(right.k_m), r0_m(right.r0_m), Bz_m(right.Bz_m),
rMin_m(right.rMin_m), rMax_m(right.rMax_m), phiStart_m(right.phiStart_m),
phiEnd_m(right.phiEnd_m), azimuthalExtent_m(right.azimuthalExtent_m),
centre_m(right.centre_m), endField_m(right.endField_m->clone()),
dfCoefficients_m(right.dfCoefficients_m) {
RefPartBunch_m = right.RefPartBunch_m;
setElType(isDrift);
Bz_m = right.Bz_m;
r0_m = right.r0_m;
}
SpiralSector::~SpiralSector() {
ScalingFFAGMagnet::~ScalingFFAGMagnet() {
}
ElementBase* SpiralSector::clone() const {
return new SpiralSector(*this);
ElementBase* ScalingFFAGMagnet::clone() const {
ScalingFFAGMagnet* magnet = new ScalingFFAGMagnet(*this);
magnet->initialise();
return magnet;
}
EMField &SpiralSector::getField() {
EMField &ScalingFFAGMagnet::getField() {
return dummy;
}
const EMField &SpiralSector::getField() const {
const EMField &ScalingFFAGMagnet::getField() const {
return dummy;
}
bool SpiralSector::apply(const size_t &i, const double &t,
bool ScalingFFAGMagnet::apply(const size_t &i, const double &t,
Vector_t &E, Vector_t &B) {
return apply(RefPartBunch_m->R[i], RefPartBunch_m->P[i], t, E, B);
}
void SpiralSector::initialise(PartBunch *bunch, double &startField, double &endField) {
RefPartBunch_m = bunch;
void ScalingFFAGMagnet::initialise() {
calculateDfCoefficients();
planarArcGeometry_m.setElementLength(r0_m*phiEnd_m); // length = phi r
planarArcGeometry_m.setCurvature(1./r0_m);
}
void ScalingFFAGMagnet::initialise(PartBunch *bunch, double &startField, double &endField) {
RefPartBunch_m = bunch;
initialise();
}
void SpiralSector::finalise() {
void ScalingFFAGMagnet::finalise() {
RefPartBunch_m = NULL;
}
bool SpiralSector::bends() const {
bool ScalingFFAGMagnet::bends() const {
return true;
}
BGeometryBase& SpiralSector::getGeometry() {
BGeometryBase& ScalingFFAGMagnet::getGeometry() {
return planarArcGeometry_m;
}
const BGeometryBase& SpiralSector::getGeometry() const {
const BGeometryBase& ScalingFFAGMagnet::getGeometry() const {
return planarArcGeometry_m;
}
void SpiralSector::accept(BeamlineVisitor& visitor) const {
visitor.visitSpiralSector(*this);
void ScalingFFAGMagnet::accept(BeamlineVisitor& visitor) const {
visitor.visitScalingFFAGMagnet(*this);
}
bool SpiralSector::apply(const Vector_t &R, const Vector_t &P,
const double &t, Vector_t &E, Vector_t &B) {
bool ScalingFFAGMagnet::getFieldValue(const Vector_t &R, Vector_t &B) const {
Vector_t pos = R - centre_m;
double r = sqrt(pos[0]*pos[0]+pos[2]*pos[2]);
double phi = atan2(pos[2], pos[0]);
double phi = atan2(pos[0], pos[2]);
Vector_t posCyl(r, pos[1], phi);
Vector_t bCyl(0., 0., 0.);
applyCylindrical(posCyl, P, t, E, bCyl);
bool outOfBounds = getFieldValueCylindrical(posCyl, 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];
}
}
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 outOfBounds;
}
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 ScalingFFAGMagnet::getFieldValueCylindrical(const Vector_t &pos, Vector_t &B) const {
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 z = pos[1];
double phi = pos[2];
if (r < rMin_m || r > rMax_m) {
return true;
}
double normRadius = r/r0_m;
double g = tanDelta_m*log(normRadius);
double phiSpiral = phi - g;
double phiSpiral = phi - g + tanDelta_m - phiStart_m;
double h = pow(normRadius, k_m)*Bz_m;
if (phiSpiral < -azimuthalExtent_m || phiSpiral > azimuthalExtent_m) {
return true;
}
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] << " ";
fringeDerivatives[i] = endField_m->function(phiSpiral, i); // d^i_phi f
}
std::cerr << std::endl;
for (size_t n = 0; n < dfCoefficients_m.size(); n += 2) {
double f2n = 0;
Vector_t deltaB;
......@@ -158,17 +152,47 @@ bool SpiralSector::applyCylindrical(const Vector_t &pos, const Vector_t &P,
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;
deltaB[1] = f2n*h*pow(z/r, n); // Bz = sum(f_2n * h * (z/r)^2n
deltaB[2] = f2nplus1*h*pow(z/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(z/r, n+1);
B += deltaB;
}
std::cerr << B << std::endl;
return true;
return false;
}
bool ScalingFFAGMagnet::apply(const Vector_t &R, const Vector_t &P,
const double &t, Vector_t &E, Vector_t &B) {
return getFieldValue(R, B);
}
void ScalingFFAGMagnet::calculateDfCoefficients() {
dfCoefficients_m = std::vector<std::vector<double> >(maxOrder_m);
if (maxOrder_m == 0) {
return;
}
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]/(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-n)*(k_m-n)/(n+1)*dfCoefficients_m[n][i]/(n+2);
}
for(size_t i = 0; i < dfCoefficients_m[n+1].size(); ++i) { // i indexes the derivative
dfCoefficients_m[n+2][i] += 2*(k_m-n)*tanDelta_m*dfCoefficients_m[n+1][i]/(n+2);
dfCoefficients_m[n+2][i+1] -= (1+tanDelta_m*tanDelta_m)*dfCoefficients_m[n+1][i]/(n+2);
}
}
}
void SpiralSector::setEndField(endfieldmodel::EndFieldModel* endField) {
void ScalingFFAGMagnet::setEndField(endfieldmodel::EndFieldModel* endField) {
if (endField_m != NULL) {
delete endField_m;
}
......
......@@ -30,22 +30,22 @@
#include "AbsBeamline/EndFieldModel/EndFieldModel.h"
#include "AbsBeamline/Component.h"
#ifndef ABSBEAMLINE_SPIRALSECTOR_H
#define ABSBEAMLINE_SPIRALSECTOR_H
#ifndef ABSBEAMLINE_ScalingFFAGMagnet_H
#define ABSBEAMLINE_ScalingFFAGMagnet_H
/** Sector bending magnet with an FFAG-style field index and spiral end shape
*/
class SpiralSector : public Component {
class ScalingFFAGMagnet : public Component {
public:
/** Construct a new SpiralSector
/** Construct a new ScalingFFAGMagnet
*
* \param name User-defined name of the SpiralSector
* \param name User-defined name of the ScalingFFAGMagnet
*/
explicit SpiralSector(const std::string &name);
explicit ScalingFFAGMagnet(const std::string &name);
/** Destructor - deletes map */
~SpiralSector();
~ScalingFFAGMagnet();
/** Inheritable copy constructor */
ElementBase* clone() const;
......@@ -65,28 +65,33 @@ class SpiralSector : public Component {
*
* \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 t not used
* \param E not used
* \param B calculated magnetic field
* \returns true if particle is outside the field map, else false
*/
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 cartesian coordinates
*
* \param R position in the local coordinate system of the bend, in
* cartesian coordinates defined like (x, y, z)
* \param B calculated magnetic field defined like (Bx, By, Bz)
* \returns true if particle is outside the field map, else false
*/
bool getFieldValue(const Vector_t &R, Vector_t &B) const;
/** 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);
bool getFieldValueCylindrical(const Vector_t &R, Vector_t &B) const;
/** Initialise the SpiralSector
/** Initialise the ScalingFFAGMagnet
*
* \param bunch the global bunch object
* \param startField not used
......@@ -94,10 +99,17 @@ class SpiralSector : public Component {
*/
void initialise(PartBunch *bunch, double &startField, double &endField);
/** Finalise the SpiralSector - sets bunch to NULL */
/** Initialise the ScalingFFAGMagnet
*
* Sets up the field expansion and the geometry; call after changing any
* field parameters
*/
void initialise();
/** Finalise the ScalingFFAGMagnet - sets bunch to NULL */
void finalise();
/** Return true - SpiralSector always bends the reference particle */
/** Return true - ScalingFFAGMagnet always bends the reference particle */
inline bool bends() const;
/** Not implemented */
......@@ -150,24 +162,70 @@ class SpiralSector : public Component {
/** Get the fringe field
*
* Returns the fringe field model; SpiralSector retains ownership of the
* Returns the fringe field model; ScalingFFAGMagnet retains ownership of the
* returned memory.
*/
endfieldmodel::EndFieldModel* getEndField() const {return endField_m;}
/** Set the fringe field
*
* - endField: the new fringe field; SpiralSector takes ownership of the
* - endField: the new fringe field; ScalingFFAGMagnet takes ownership of the
* memory associated with endField.
*/
void setEndField(endfieldmodel::EndFieldModel* endField);
/** Get the maximum pole modelled in the off-midplane expansion; 0 is dipole; 2 is quadrupole; etc */
/** Get the maximum power of y modelled in the off-midplane expansion;
*/
size_t getMaxOrder() const {return maxOrder_m;}
/** Set the maximum pole modelled in the off-midplane expansion; 0 is dipole; 2 is quadrupole; etc */
/** Set the maximum power of y modelled in the off-midplane expansion;
*/
void setMaxOrder(size_t maxOrder) {maxOrder_m = maxOrder;}
/** Get the offset of the magnet centre from the start
*/
double getPhiStart() const {return phiStart_m;}
/** Set the offset of the magnet centre from the start
*/
void setPhiStart(double phiStart) {phiStart_m = phiStart;}
/** Get the offset of the magnet end from the start
*/
double getPhiEnd() const {return phiEnd_m;}
/** Set the offset of the magnet end from the start
*/
void setPhiEnd(double phiEnd) {phiEnd_m = phiEnd;}
/** Get the maximum radius
*/
double getRMin() const {return rMin_m;}
/** Set the maximum radius
*/
void setRMin(double rMin) {rMin_m = rMin;}
/** Get the maximum radius
*/
double getRMax() const {return rMax_m;}
/** Set the maximum radius
*/
void setRMax(double rMax) {rMax_m = rMax;}
/** Get the maximum azimuthal displacement from \psi=0
*/
double getAzimuthalExtent() const {return azimuthalExtent_m;}