Forked from
OPAL / OPAL-X / src
188 commits behind the upstream repository.
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
Ring.cpp 14.92 KiB
/*
* Copyright (c) 2012-2014, 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/Ring.h"
#include <cmath>
#include <limits>
#include <sstream>
#include "Utility/Inform.h" // ippl
#include "AbsBeamline/BeamlineVisitor.h"
#include "Fields/EMField.h"
#include "PartBunch/PartBunch.hpp"
#include "Physics/Physics.h"
#include "Structure/LossDataSink.h"
// fairly generous tolerance here... maybe too generous? Maybe should be
// user parameter?
const double Ring::lengthTolerance_m = 1e-2;
const double Ring::angleTolerance_m = 1e-2;
Ring::Ring(std::string ring)
: Component(ring),
planarArcGeometry_m(1, 1),
refPartBunch_m(nullptr),
lossDS_m(nullptr),
beamRInit_m(0.),
beamPRInit_m(0.),
beamPhiInit_m(0.),
latticeRInit_m(0.),
latticePhiInit_m(0.),
latticeThetaInit_m(0.),
isLocked_m(false),
isClosed_m(true),
symmetry_m(0),
cyclHarm_m(0),
rfFreq_m(0),
phiStep_m(Physics::pi / 100.),
ringSections_m(),
section_list_m() {
setRefPartBunch(nullptr);
}
void Ring::accept(BeamlineVisitor& visitor) const {
visitor.visitRing(*this);}
Ring::Ring(const Ring& ring)
: Component(ring.getName()),
planarArcGeometry_m(ring.planarArcGeometry_m),
lossDS_m(nullptr),
beamRInit_m(ring.beamRInit_m),
beamPRInit_m(ring.beamPRInit_m),
beamPhiInit_m(ring.beamPhiInit_m),
latticeRInit_m(ring.latticeRInit_m),
latticePhiInit_m(ring.latticePhiInit_m),
latticeThetaInit_m(ring.latticeThetaInit_m),
willDoAperture_m(ring.willDoAperture_m),
minR2_m(ring.minR2_m),
maxR2_m(ring.maxR2_m),
isLocked_m(ring.isLocked_m),
isClosed_m(ring.isClosed_m),
symmetry_m(ring.symmetry_m),
cyclHarm_m(ring.cyclHarm_m),
phiStep_m(ring.phiStep_m),
ringSections_m(),
section_list_m(ring.section_list_m.size()) {
setRefPartBunch(ring.refPartBunch_m);
if (ring.lossDS_m != nullptr)
throw GeneralClassicException(
"Ring::Ring(const Ring&)",
"Can't copy construct LossDataSink so copy constructor fails");
for (size_t i = 0; i < section_list_m.size(); ++i) {
section_list_m[i] = new RingSection(*ring.section_list_m[i]);
}
buildRingSections();
}
Ring::~Ring() {
for (size_t i = 0; i < section_list_m.size(); ++i)
delete section_list_m[i];
}
bool Ring::apply(
const size_t& id, const double& t, Vector_t<double, 3>& E, Vector_t<double, 3>& B) {
std::shared_ptr<ParticleContainer_t> pc = RefPartBunch_m->getParticleContainer();
auto Rview = pc->R.getView();
auto Pview = pc->P.getView();
const Vector_t<double, 3> R = Rview(id);
const Vector_t<double, 3> P = Pview(id);
bool flagNeedUpdate = apply(R, P, t, E, B);
if (flagNeedUpdate) {
Inform gmsgALL("OPAL ", INFORM_ALL_NODES);
auto Qview = pc->Q.getView();
auto Mview = pc->M.getView();
auto Binview = pc->Bin.getView();
const double Q = Qview(id);
const double M = Mview(id);
// const short int Bin = Binview(id);
gmsgALL << getName() << ": particle " << id << " at " << R
<< " m out of the field map boundary" << endl;
lossDS_m->addParticle(OpalParticle(id, R * Vector_t<double, 3>(1000.0), P, t, Q, M));
Binview(id) = -1;
}
return flagNeedUpdate;
}
bool Ring::apply( const Vector_t<double, 3>& R, const Vector_t<double, 3>& /*P*/, const double& t,
Vector_t<double, 3>& E, Vector_t<double, 3>& B) {
B = Vector_t<double, 3>({0.0, 0.0, 0.0});
E = Vector_t<double, 3>({0.0, 0.0, 0.0});
std::vector<RingSection*> sections =
getSectionsAt(R); // I think this doesn't actually use R -DW
bool outOfBounds = true;
// assume field maps don't set B, E to 0...
if (willDoAperture_m) {
double rSquared = R[0] * R[0] + R[1] * R[1];
if (rSquared < minR2_m || rSquared > maxR2_m) {
return true;
}
}
for (size_t i = 0; i < sections.size(); ++i) {
Vector_t<double, 3> B_temp({0.0, 0.0, 0.0});
Vector_t<double, 3> E_temp({0.0, 0.0, 0.0});
// Super-TEMP! cyclotron tracker now uses m internally, have to change to mm here to match
// old field limits -DW
outOfBounds &= sections[i]->getFieldValue(
R * Vector_t<double, 3>(1000.0),
refPartBunch_m->get_centroid() * Vector_t<double, 3>(1000.0), t, E_temp, B_temp);
B += (scale_m * B_temp);
E += (scale_m * E_temp);
}
return outOfBounds;
}
void Ring::setLossDataSink(LossDataSink* sink) {
// lossDS_m is a borrowed pointer - we never delete it
lossDS_m = sink;
}
void Ring::getDimensions(double& /*zBegin*/, double& /*zEnd*/) const {
throw GeneralClassicException("Ring::getDimensions", "Cannot get s-dimension of a ring");
}
void Ring::initialise(PartBunch_t* bunch) {
online_m = true;
setRefPartBunch(bunch);
setLossDataSink(new LossDataSink(getName(), false));
}
void Ring::initialise(PartBunch_t* bunch, double& /*startField*/, double& /*endField*/) {
initialise(bunch);
}
void Ring::finalise() {
lossDS_m->save();
online_m = false;
setLossDataSink(nullptr);
}
void Ring::setRefPartBunch(PartBunch_t* bunch) {
RefPartBunch_m = bunch; // inherited from Component
refPartBunch_m = bunch; // private data (obeys style guide)
}
std::vector<RingSection*> Ring::getSectionsAt(const Vector_t<double, 3>& /*r*/) {
return section_list_m;
}
Rotation3D Ring::getRotationStartToEnd(Euclid3D delta) const {
// rotation from start to end in local coordinates
// Euclid3D/Rotation3D doesnt have a "getAngle" method so we use fairly
// obscure technique to extract it.
Vector3D rotationTest(1., 0., 0.);
rotationTest = delta.getRotation() * rotationTest; double deltaAngle = std::atan2(rotationTest(2), rotationTest(0));
Rotation3D elementRotation = Rotation3D::ZRotation(deltaAngle);
return elementRotation;
}
void Ring::checkMidplane(Euclid3D delta) const {
if (std::abs(delta.getVector()(2)) > lengthTolerance_m
|| std::abs(delta.getRotation().getAxis()(0)) > angleTolerance_m
|| std::abs(delta.getRotation().getAxis()(1)) > angleTolerance_m) {
throw GeneralClassicException(
"Ring::checkMidplane",
std::string("Placement of elements out of the ") + "midplane is not supported by Ring");
}
}
void Ring::rotateToCyclCoordinates(Euclid3D& delta) const {
Vector3D v = delta.getVector();
Vector3D r = delta.getRotation().getAxis();
// Euclid3D euclid(v(0), -v(2), v(1), r(0), -r(2), r(1));
Euclid3D euclid(v(2), v(0), -v(1), r(2), r(0), -r(1));
delta = euclid;
}
Vector_t<double, 3> Ring::getNextPosition() const {
if (!section_list_m.empty()) {
return section_list_m.back()->getEndPosition();
}
return Vector_t<double, 3>(
{latticeRInit_m * std::sin(latticePhiInit_m), latticeRInit_m * std::cos(latticePhiInit_m),
0.});
}
Vector_t<double, 3> Ring::getNextNormal() const {
if (!section_list_m.empty()) {
return section_list_m.back()->getEndNormal();
}
return Vector_t<double, 3>(
{std::cos(latticePhiInit_m + latticeThetaInit_m),
-std::sin(latticePhiInit_m + latticeThetaInit_m), 0.});
}
void Ring::appendElement(const Component& element) {
if (isLocked_m) {
throw GeneralClassicException(
"Ring::appendElement",
"Attempt to append element " + element.getName() + " when ring is locked");
}
// delta is transform from start of bend to end with x, z as horizontal
// I failed to get Rotation3D to work so use rotations written by hand.
// Probably an error in my call to Rotation3D.
Euclid3D delta = element.getGeometry().getTotalTransform();
rotateToCyclCoordinates(delta);
checkMidplane(delta);
RingSection* section = new RingSection();
Vector_t<double, 3> startPos = getNextPosition();
Vector_t<double, 3> startNorm = getNextNormal();
section->setComponent(dynamic_cast<Component*>(element.clone()));
section->setStartPosition(startPos);
section->setStartNormal(startNorm);
double startF = std::atan2(startNorm(1), startNorm(0));
Vector_t<double, 3> endPos =
Vector_t<double, 3>(
{+delta.getVector()(0) * std::cos(startF) - delta.getVector()(1) * std::sin(startF),
+delta.getVector()(0) * std::sin(startF) + delta.getVector()(1) * std::cos(startF), 0})
+ startPos;
section->setEndPosition(endPos);
double endF = delta.getRotation().getAxis()(2); //+
// atan2(delta.getVector()(1), delta.getVector()(0));
Vector_t<double, 3> endNorm = Vector_t<double, 3>(
{+startNorm(0) * std::cos(endF) + startNorm(1) * std::sin(endF),
-startNorm(0) * std::sin(endF) + startNorm(1) * std::cos(endF), 0});
section->setEndNormal(endNorm);
section->setComponentPosition(startPos);
section->setComponentOrientation(Vector_t<double, 3>({0, 0, startF}));
section_list_m.push_back(section);
double dphi = atan2(startNorm(0), startNorm(1));
Inform msg("OPAL");
msg << "* Added " << element.getName() << " to Ring" << endl;
msg << "* Start position (" << section->getStartPosition()(0) << ", "
<< section->getStartPosition()(1) << ", " << section->getStartPosition()(2) << ") normal ("
<< section->getStartNormal()(0) << ", " << section->getStartNormal()(1) << ", "
<< section->getStartNormal()(2) << "), phi " << dphi << endl;
msg << "* End position (" << section->getEndPosition()(0) << ", "
<< section->getEndPosition()(1) << ", " << section->getEndPosition()(2) << ") normal ("
<< section->getEndNormal()(0) << ", " << section->getEndNormal()(1) << ", "
<< section->getEndNormal()(2) << ")" << endl;
}
void Ring::lockRing() {
Inform msg("OPAL ");
if (isLocked_m) {
throw GeneralClassicException(
"Ring::lockRing", "Attempt to lock ring when it is already locked");
}
// check for any elements at all
size_t sectionListSize = section_list_m.size();
if (sectionListSize == 0) {
throw GeneralClassicException("Ring::lockRing", "Failed to find any elements in Ring");
}
// Apply symmetry properties; I think it is fastest to just duplicate
// elements rather than try to do some angle manipulations when we do field
// lookups because we do field lookups in Cartesian coordinates in general.
msg << "Applying symmetry to Ring" << endl;
for (int i = 1; i < symmetry_m; ++i) {
for (size_t j = 0; j < sectionListSize; ++j) {
appendElement(*section_list_m[j]->getComponent());
}
}
// resetAzimuths();
// Check that the ring is closed within tolerance and close exactly
if (isClosed_m)
checkAndClose();
buildRingSections();
isLocked_m = true;
for (size_t i = 0; i < section_list_m.size(); i++) {
}
}
void Ring::resetAzimuths() {
for (size_t i = 0; i < section_list_m.size(); ++i) {
Vector_t<double, 3> startPos = section_list_m[i]->getEndPosition();
Vector_t<double, 3> startDir = section_list_m[i]->getEndNormal();
Vector_t<double, 3> endPos = section_list_m[i]->getEndPosition();
Vector_t<double, 3> endDir = section_list_m[i]->getEndNormal();
if (!section_list_m[i]->isOnOrPastStartPlane(endPos)) {
section_list_m[i]->setEndPosition(startPos);
section_list_m[i]->setEndNormal(startDir);
section_list_m[i]->setStartPosition(endPos);
section_list_m[i]->setStartNormal(endDir);
}
}
}
void Ring::checkAndClose() {
Vector_t<double, 3> first_pos = section_list_m[0]->getStartPosition();
Vector_t<double, 3> first_norm = section_list_m[0]->getStartNormal();
Vector_t<double, 3> last_pos = section_list_m.back()->getEndPosition();
Vector_t<double, 3> last_norm = section_list_m.back()->getEndNormal();
for (int i = 0; i < 3; ++i) {
if (std::abs(first_pos(i) - last_pos(i)) > lengthTolerance_m
|| std::abs(first_norm(i) - last_norm(i)) > angleTolerance_m)
throw GeneralClassicException("Ring::lockRing", "Ring is not closed");
}
section_list_m.back()->setEndPosition(first_pos);
section_list_m.back()->setEndNormal(first_norm);
}
void Ring::buildRingSections() {
size_t nSections = 2. * Physics::pi / phiStep_m + 1;
ringSections_m = std::vector<std::vector<RingSection*> >(nSections);
for (size_t i = 0; i < ringSections_m.size(); ++i) {
double phi0 = i * phiStep_m;
double phi1 = (i + 1) * phiStep_m;
// std::cerr << phi0 << " " << phi1 << std::endl;
for (size_t j = 0; j < section_list_m.size(); ++j) {
if (section_list_m[j]->doesOverlap(phi0, phi1))
ringSections_m[i].push_back(section_list_m[j]);
}
}
}
RingSection* Ring::getLastSectionPlaced() const {
if (section_list_m.empty()) {
return nullptr;
}
return section_list_m.back();
}
bool Ring::sectionCompare(RingSection const* const sec1, RingSection const* const sec2) {
return sec1 > sec2;
}
void Ring::setRingAperture(double minR, double maxR) {
if (minR < 0 || maxR < 0) {
throw GeneralClassicException(
"Ring::setRingAperture", "Could not parse negative or undefined aperture limit");
}
willDoAperture_m = true;
minR2_m = minR * minR;
maxR2_m = maxR * maxR;
}