Commit 11ebb0e6 authored by ext-rogers_c's avatar ext-rogers_c
Browse files

Generic geometry elements for rings and 3D field maps in ring geometry

parent 1d28d3c6
......@@ -58,6 +58,9 @@ classic/5.0/src/AbsBeamline/RFQuadrupole.cpp -text
classic/5.0/src/AbsBeamline/RFQuadrupole.h -text
classic/5.0/src/AbsBeamline/SBend.cpp -text
classic/5.0/src/AbsBeamline/SBend.h -text
classic/5.0/src/AbsBeamline/SBend3D.cpp -text
classic/5.0/src/AbsBeamline/SBend3D.h -text
classic/5.0/src/AbsBeamline/SectorFieldMapComponent.h -text
classic/5.0/src/AbsBeamline/Separator.cpp -text
classic/5.0/src/AbsBeamline/Separator.h -text
classic/5.0/src/AbsBeamline/Septum.cpp -text
......@@ -347,6 +350,25 @@ classic/5.0/src/Fields/Fieldmap.icc -text
classic/5.0/src/Fields/NullField.cpp -text
classic/5.0/src/Fields/NullField.h -text
classic/5.0/src/Fields/OscillatingField.h -text
classic/5.0/src/Fields/SectorMagneticFieldMap.cpp -text
classic/5.0/src/Fields/SectorMagneticFieldMap.h -text
classic/5.0/src/Fields/SectorMagneticFieldMap/CMakeLists.txt -text
classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo1d.cpp -text
classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo1d.h -text
classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo3d.cpp -text
classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo3d.h -text
classic/5.0/src/Fields/SectorMagneticFieldMap/Makefile -text
classic/5.0/src/Fields/SectorMagneticFieldMap/Mesh-inl.icc -text
classic/5.0/src/Fields/SectorMagneticFieldMap/Mesh.cpp -text
classic/5.0/src/Fields/SectorMagneticFieldMap/Mesh.h -text
classic/5.0/src/Fields/SectorMagneticFieldMap/SectorField.cpp -text
classic/5.0/src/Fields/SectorMagneticFieldMap/SectorField.h -text
classic/5.0/src/Fields/SectorMagneticFieldMap/ThreeDGrid.cpp -text
classic/5.0/src/Fields/SectorMagneticFieldMap/ThreeDGrid.h -text
classic/5.0/src/Fields/SectorMagneticFieldMap/TriLinearInterpolator.cpp -text
classic/5.0/src/Fields/SectorMagneticFieldMap/TriLinearInterpolator.h -text
classic/5.0/src/Fields/SectorMagneticFieldMap/VectorMap.cpp -text
classic/5.0/src/Fields/SectorMagneticFieldMap/VectorMap.h -text
classic/5.0/src/Fields/StaticElectricField.cpp -text
classic/5.0/src/Fields/StaticElectricField.h -text
classic/5.0/src/Fields/StaticMagneticField.cpp -text
......@@ -741,8 +763,15 @@ src/Elements/OpalRBend.cpp -text
src/Elements/OpalRBend.h -text
src/Elements/OpalRCollimator.cpp -text
src/Elements/OpalRCollimator.h -text
src/Elements/OpalRing.cpp -text
src/Elements/OpalRing.h -text
src/Elements/OpalRingDefinition.cpp -text
src/Elements/OpalRingDefinition.h -text
src/Elements/OpalRingTest.h -text
src/Elements/OpalSBend.cpp -text
src/Elements/OpalSBend.h -text
src/Elements/OpalSBend3D.cpp -text
src/Elements/OpalSBend3D.h -text
src/Elements/OpalSRot.cpp -text
src/Elements/OpalSRot.h -text
src/Elements/OpalSeparator.cpp -text
......@@ -1025,6 +1054,8 @@ src/Utilities/OpalField.cpp -text
src/Utilities/OpalField.h -text
src/Utilities/OpalFilter.cpp -text
src/Utilities/OpalFilter.h -text
src/Utilities/OpalRingSection.cpp -text
src/Utilities/OpalRingSection.h -text
src/Utilities/OpalSection.cpp -text
src/Utilities/OpalSection.h -text
src/Utilities/Options.h -text
......
......@@ -53,6 +53,7 @@ class RFCavity;
class TravelingWave;
class RFQuadrupole;
class SBend;
class SBend3D;
class Cyclotron;
class Separator;
class Septum;
......@@ -60,6 +61,7 @@ class Solenoid;
class ParallelPlate;
class CyclotronValley;
class Stripper;
class OpalRing;
// Integrators.
class Integrator;
class MapIntegrator;
......@@ -117,6 +119,9 @@ public:
/// Apply the algorithm to a drift space.
virtual void visitDrift(const Drift &) = 0;
/// Apply the algorithm to an OpalRing
virtual void visitOpalRing(const OpalRing &) = 0;
/// Apply the algorithm to a cyclotron.
virtual void visitCyclotron(const Cyclotron &) = 0;
......@@ -153,6 +158,9 @@ public:
/// Apply the algorithm to a sector bend.
virtual void visitSBend(const SBend &) = 0;
/// Apply the algorithm to a Sector Bend with 3D field map.
virtual void visitSBend3D(const SBend3D &) = 0;
/// Apply the algorithm to an electrostatic separator.
virtual void visitSeparator(const Separator &) = 0;
......
......@@ -24,6 +24,7 @@ set (_SRCS
TravelingWave.cpp
RFQuadrupole.cpp
SBend.cpp
SBend3D.cpp
Separator.cpp
Septum.cpp
Solenoid.cpp
......
......@@ -270,7 +270,7 @@ bool Cyclotron::apply(const size_t &id, const double &t, Vector_t &E, Vector_t &
if (zpos > maxz_m || zpos < minz_m || rpos > maxr_m || rpos < minr_m){
flagNeedUpdate = true;
Inform gmsgALL("OPAL ", INFORM_ALL_NODES);
Inform gmsgALL("OPAL ", INFORM_ALL_NODES);
gmsgALL<<getName() << ": particle "<< id <<" out of the global aperture of cyclotron!"<< endl;
} else{
......
......@@ -106,7 +106,6 @@ struct BPositions {
// This class defines the abstract interface for a Cyclotron.
class Cyclotron: public Component {
public:
/// Constructor with given name.
......
/*
* Copyright (c) 2012, 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 "AbsBeamline/SBend3D.h"
#include "Algorithms/PartBunch.h"
#include "AbsBeamline/BeamlineVisitor.h"
SBend3D::SBend3D(const std::string &name)
: Component(name), map_m(NULL),
planarArcGeometry_m(1., 1.), fieldUnits_m(1.), lengthUnits_m(1.),
dummy() {
}
SBend3D::SBend3D(const SBend3D &right)
: Component(right), map_m(NULL),
planarArcGeometry_m(right.planarArcGeometry_m),
fieldUnits_m(right.fieldUnits_m), lengthUnits_m(right.lengthUnits_m),
dummy() {
RefPartBunch_m = right.RefPartBunch_m;
if (right.map_m != NULL)
map_m = new SectorMagneticFieldMap(*right.map_m);
}
SBend3D::~SBend3D() {
delete map_m;
}
ElementBase* SBend3D::clone() const {
return new SBend3D(*this);
}
EMField &SBend3D::getField() {
return dummy;
}
const EMField &SBend3D::getField() const {
return dummy;
}
bool SBend3D::apply(const size_t &i, const double &t, double E[], double B[]) {
Vector_t Ev(0, 0, 0), Bv(0, 0, 0);
if (apply(i, t, Ev, Bv))
return true;
E[0] = Ev(0);
E[1] = Ev(1);
E[2] = Ev(2);
B[0] = Bv(0);
B[1] = Bv(1);
B[2] = Bv(2);
return false;
}
bool SBend3D::apply(const size_t &i, const double &t,
Vector_t &E, Vector_t &B) {
return apply(RefPartBunch_m->R[i], RefPartBunch_m->get_centroid(), t,
E, B);
}
bool SBend3D::apply(const Vector_t &R, const Vector_t &centroid,
const double &t, Vector_t &E, Vector_t &B) {
return map_m->getFieldstrength(R, E, B);
}
void SBend3D::initialise(PartBunch *bunch, double &startField, double &endField,
const double &scaleFactor) {
RefPartBunch_m = bunch;
}
void SBend3D::finalise() {
RefPartBunch_m = NULL;
}
bool SBend3D::bends() const {
return true;
}
Geometry& SBend3D::getGeometry() {
return planarArcGeometry_m;
}
const Geometry& SBend3D::getGeometry() const {
return planarArcGeometry_m;
}
void SBend3D::setFieldMapFileName(std::string name) {
if (map_m != NULL) {
delete map_m;
map_m = NULL;
}
if (name != "") {
map_m = new SectorMagneticFieldMap
(name, "Dipole", lengthUnits_m, fieldUnits_m);
double r_curv = (map_m->getPolarBoundingBoxMax()[0]+
map_m->getPolarBoundingBoxMin()[0])/2.;
double delta_phi = (map_m->getPolarBoundingBoxMax()[2]-
map_m->getPolarBoundingBoxMin()[2]);
planarArcGeometry_m.setElementLength(r_curv*delta_phi);
planarArcGeometry_m.setCurvature(1./r_curv);
// std::cerr << "RCURV " << r_curv << " DPHI " << delta_phi << " DELTA X: "
// << planarArcGeometry_m.getTotalTransform().getVector()(0) << " Y: "
// << planarArcGeometry_m.getTotalTransform().getVector()(1) << " Z: "
// << planarArcGeometry_m.getTotalTransform().getVector()(2) << std::endl;
}
}
void SBend3D::accept(BeamlineVisitor& visitor) const {
visitor.visitSBend3D(*this);
}
/*
* Copyright (c) 2012, 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 "Fields/BMultipoleField.h"
#include "Fields/SectorMagneticFieldMap.h"
#include "BeamlineGeometry/PlanarArcGeometry.h"
#include "AbsBeamline/Component.h"
#ifndef ABSBEAMLINE_SBEND3D_H
#define ABSBEAMLINE_SBEND3D_H
/** Sector bending magnet from a 3D field map
*
* SBend3D provides interface between Component type and routines for linear
* interpolation from a 3D field map in a sector geometry.
*
* \param planarArcGeometry_m Cell geometry is a PlanarArcGeometry with radius
* of curvature equal to the radius at the middle of the field map and
* angle given by the field map opening angle.
* \param map_m the actual field map class (with associated interpolation and
* IO routines).
* \param field_m a dummy field - never used, only put in to obey inheritance
* requirements
*/
class SBend3D : public Component {
public:
/** Construct a new SBend3D
*
* \param name User-defined name of the SBend3D
*/
explicit SBend3D(const std::string &name);
/** Copy constructor */
SBend3D(const SBend3D &right);
/** Destructor - deletes map */
~SBend3D();
/** Inheritable copy constructor */
ElementBase* clone() const;
/** Calculate the field at the position of the ith particle
*
* \param i index of the particle event; field is calculated at this
* position
* \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
* \returns true if particle is outside the field map
*/
bool apply(const size_t &i, const double &t, double E[], double B[]);
/** Calculate the field at the position of the ith particle
*
* \param i index of the particle event; field is calculated at this
* position
* \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
* \returns true if particle is outside the field map
*/
bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
/** Calculate the field at some arbitrary position
*
* \param R position in the local coordinate system of the bend
* \param centroid 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
* \returns true if particle is outside the field map, else false
*/
bool apply(const Vector_t &R, const Vector_t &centroid, const double &t,
Vector_t &E, Vector_t &B);
/** Initialise the SBend3D
*
* \param bunch the global bunch object
* \param startField not used
* \param endField not used
* \param scaleFactor not used
*/
void initialise(PartBunch *bunch, double &startField, double &endField,
const double &scaleFactor);
/** Finalise the SBend3D - sets bunch to NULL */
void finalise();
/** Return true - SBend3D always bends the reference particle */
inline bool bends() const;
/** Not implemented */
void getDimensions(double &zBegin, double &zEnd) const {}
/** Return the cell geometry */
Geometry& getGeometry();
/** Return the cell geometry */
const Geometry& getGeometry() const;
/** Return a dummy (0.) field value (what is this for?) */
EMField &getField();
/** Return a dummy (0.) field value (what is this for?) */
const EMField &getField() const;
/** Accept a beamline visitor */
void accept(BeamlineVisitor& visitor) const;
/** Set the field map file
*
* @param name the file name of the field map
*
* This generates a new field map object with the given name. If name == ""
* then sets the field map to NULL.
*/
void setFieldMapFileName(std::string name);
/** Get the file name of the field map */
inline std::string getFieldMapFileName() const;
/** Get the scale factor */
double getLengthUnits() const {return lengthUnits_m;}
/** Set the scale factor */
void setLengthUnits(double lengthUnits) {lengthUnits_m = lengthUnits;}
/** Get the scale factor */
double getFieldUnits() const {return fieldUnits_m;}
/** Set the scale factor */
void setFieldUnits(double fieldUnits) {fieldUnits_m = fieldUnits;}
private:
SectorMagneticFieldMap* map_m;
// Geometry
PlanarArcGeometry planarArcGeometry_m;
// Units used for field maps
double fieldUnits_m;
double lengthUnits_m;
// Not implemented (I think this is something to do with error fields?)
BMultipoleField dummy;
};
std::string SBend3D::getFieldMapFileName() const {
return map_m->getFieldMapFileName();
}
#endif
class SectorFieldMapComponent: public SBend {
public:
/// Constructor with given name.
SectorFieldMapComponent(const std::string &name);
SectorFieldMapComponent();
SectorFieldMapComponent(const SectorFieldMapComponent &right);
~SectorFieldMapComponent();
/// Return field.
// The representation of the electro-magnetic field of the component
// (version for non-constant object).
virtual EMField &getField() = 0;
/// Return field.
// The representation of the electro-magnetic field of the component
// (version for constant object).
virtual const EMField &getField() const = 0;
virtual bool apply(const size_t &i, const double &t, double E[], double B[]) = 0;
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) = 0;
virtual bool apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B) = 0;
virtual void initialise(PartBunch *bunch, double &startField, double &endField, const double &scaleFactor) = 0;
virtual void finalise() = 0;
virtual bool bends() const = 0;
virtual void getDimensions(double &zBegin, double &zEnd) const = 0;
private:
};
......@@ -40,6 +40,7 @@
#include "AbsBeamline/TravelingWave.h"
#include "AbsBeamline/RFQuadrupole.h"
#include "AbsBeamline/SBend.h"
#include "AbsBeamline/SBend3D.h"
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
......@@ -59,6 +60,8 @@
#include "ComponentWrappers/SBendWrapper.h"
#include "ComponentWrappers/CyclotronWrapper.h"
#include "Elements/OpalRing.h"
// Class DefaultVisitor
// ------------------------------------------------------------------------
......@@ -138,6 +141,11 @@ void DefaultVisitor::visitMultipole(const Multipole &mult) {
}
void DefaultVisitor::visitOpalRing(const OpalRing &ring) {
applyDefault(ring);
}
void DefaultVisitor::visitPatch(const Patch &patch) {
applyDefault(patch);
}
......@@ -171,6 +179,11 @@ void DefaultVisitor::visitSBend(const SBend &bend) {
}
void DefaultVisitor::visitSBend3D(const SBend3D &bend) {
applyDefault(bend);
}
void DefaultVisitor::visitSeparator(const Separator &sep) {
applyDefault(sep);
}
......
......@@ -51,7 +51,6 @@ public:
/// Apply the algorithm to the top-level beamline.
virtual void execute();
/// Apply the algorithm to a beam-beam.
virtual void visitBeamBeam(const BeamBeam &);
......@@ -64,6 +63,9 @@ public:
/// Apply the algorithm to an cyclotron
virtual void visitCyclotron(const Cyclotron &);
/// Apply the algorithm to an opal ring..
virtual void visitOpalRing(const OpalRing &);
/// Apply the algorithm to a corrector.
virtual void visitCorrector(const Corrector &);
......@@ -109,6 +111,9 @@ public:
/// Apply the algorithm to a sector bend.
virtual void visitSBend(const SBend &);
/// Apply the algorithm to a sector bend.
virtual void visitSBend3D(const SBend3D &);
/// Apply the algorithm to a separator.
virtual void visitSeparator(const Separator &);
......
......@@ -45,6 +45,7 @@ add_subdirectory(Beamlines)
add_subdirectory(Channels)
add_subdirectory(ComponentWrappers)
add_subdirectory(Construction)
add_subdirectory(Fields/SectorMagneticFieldMap)
add_subdirectory(Fields)
add_subdirectory(Filters)
add_subdirectory(FixedAlgebra)
......
......@@ -36,10 +36,12 @@ set (_SRCS
FM1DProfile1.cpp
FM1DProfile2.cpp
FMDummy.cpp
SectorMagneticFieldMap.cpp
)
include_directories (
SectorMagneticFieldMap
${CMAKE_CURRENT_SOURCE_DIR}
)
add_cl_sources(${_SRCS})
\ No newline at end of file
add_cl_sources(${_SRCS})
/*
* Copyright (c) 2012, 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 <algorithm>
#include <vector>
#include <iostream>
#include <fstream>
#include <string>
#include "Utilities/LogicalError.h"
#include "Fields/SectorMagneticFieldMap/Mesh.h"
#include "Fields/SectorMagneticFieldMap/ThreeDGrid.h"
#include "Fields/SectorMagneticFieldMap/Interpolator3dGridTo3d.h"
#include "Fields/SectorMagneticFieldMap.h"
// allow a fairly generous phi tolerance - we don't care about phi much and
// calculation can be flaky due to ascii truncation of double and conversions