Commit af9160f9 authored by gsell's avatar gsell

Merge branch 'master' of gitorious.psi.ch:opal/src

parents b5757442 c2fc3641
// ------------------------------------------------------------------------
// $RCSfile: Bend.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Definitions for class: Bend
// Defines the abstract interface for a sector bend magnet.
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:31 $
// $Author: fci $
//
// ------------------------------------------------------------------------
#include "AbsBeamline/Bend.h"
#include "Algorithms/PartPusher.h"
#include "Algorithms/PartBunch.h"
#include "AbsBeamline/BeamlineVisitor.h"
#include "Utilities/Options.h"
#include "Fields/Fieldmap.h"
#include "AbstractObjects/OpalData.h"
#include <iostream>
#include <fstream>
extern Inform *gmsg;
// Class Bend
// ------------------------------------------------------------------------
Bend::Bend():
Component(),
messageHeader_m(" * "),
pusher_m(),
fileName_m(""),
fieldmap_m(NULL),
fast_m(false),
angle_m(0.0),
aperture_m(0.0),
designEnergy_m(0.0),
designRadius_m(0.0),
fieldAmplitude_m(0.0),
bX_m(0.0),
bY_m(0.0),
angleGreaterThanPi_m(false),
entranceAngle_m(0.0),
exitAngle_m(0.0),
fieldIndex_m(0.0),
elementEdge_m(0.0),
startField_m(0.0),
endField_m(0.0),
reinitialize_m(false),
recalcRefTraj_m(false),
length_m(0.0),
gap_m(0.0),
refTrajMapSize_m(0),
refTrajMapStepSize_m(0.0),
entranceParameter1_m(0.0),
entranceParameter2_m(0.0),
entranceParameter3_m(0.0),
exitParameter1_m(0.0),
exitParameter2_m(0.0),
exitParameter3_m(0.0),
xOriginEngeEntry_m(0.0),
zOriginEngeEntry_m(0.0),
deltaBeginEntry_m(0.0),
deltaEndEntry_m(0.0),
polyOrderEntry_m(0),
xExit_m(0.0),
zExit_m(0.0),
xOriginEngeExit_m(0.0),
zOriginEngeExit_m(0.0),
deltaBeginExit_m(0.0),
deltaEndExit_m(0.0),
polyOrderExit_m(0),
cosEntranceAngle_m(1.0),
sinEntranceAngle_m(0.0),
exitEdgeAngle_m(0.0),
cosExitAngle_m(1.0),
sinExitAngle_m(0.0) {
setElType(isDipole);
}
Bend::Bend(const Bend &right):
Component(right),
messageHeader_m(" * "),
pusher_m(right.pusher_m),
fileName_m(right.fileName_m),
fieldmap_m(right.fieldmap_m),
fast_m(right.fast_m),
angle_m(right.angle_m),
aperture_m(right.aperture_m),
designEnergy_m(right.designEnergy_m),
designRadius_m(right.designRadius_m),
fieldAmplitude_m(right.fieldAmplitude_m),
bX_m(right.bX_m),
bY_m(right.bY_m),
angleGreaterThanPi_m(right.angleGreaterThanPi_m),
entranceAngle_m(right.entranceAngle_m),
exitAngle_m(right.exitAngle_m),
fieldIndex_m(right.fieldIndex_m),
elementEdge_m(right.elementEdge_m),
startField_m(right.startField_m),
endField_m(right.endField_m),
reinitialize_m(right.reinitialize_m),
recalcRefTraj_m(right.recalcRefTraj_m),
length_m(right.length_m),
gap_m(right.gap_m),
refTrajMapX_m(right.refTrajMapX_m),
refTrajMapY_m(right.refTrajMapY_m),
refTrajMapZ_m(right.refTrajMapZ_m),
refTrajMapSize_m(right.refTrajMapSize_m),
refTrajMapStepSize_m(right.refTrajMapStepSize_m),
entranceParameter1_m(right.entranceParameter1_m),
entranceParameter2_m(right.entranceParameter2_m),
entranceParameter3_m(right.entranceParameter3_m),
exitParameter1_m(right.exitParameter1_m),
exitParameter2_m(right.exitParameter2_m),
exitParameter3_m(right.exitParameter3_m),
xOriginEngeEntry_m(right.xOriginEngeEntry_m),
zOriginEngeEntry_m(right.zOriginEngeEntry_m),
deltaBeginEntry_m(right.deltaBeginEntry_m),
deltaEndEntry_m(right.deltaEndEntry_m),
polyOrderEntry_m(right.polyOrderEntry_m),
xExit_m(right.xExit_m),
zExit_m(right.zExit_m),
xOriginEngeExit_m(right.xOriginEngeExit_m),
zOriginEngeExit_m(right.zOriginEngeExit_m),
deltaBeginExit_m(right.deltaBeginExit_m),
deltaEndExit_m(right.deltaEndExit_m),
polyOrderExit_m(right.polyOrderExit_m),
cosEntranceAngle_m(right.cosEntranceAngle_m),
sinEntranceAngle_m(right.sinEntranceAngle_m),
exitEdgeAngle_m(right.exitEdgeAngle_m),
cosExitAngle_m(right.cosExitAngle_m),
sinExitAngle_m(right.sinExitAngle_m) {
setElType(isDipole);
}
Bend::Bend(const std::string &name):
Component(name),
messageHeader_m(" * "),
pusher_m(),
fileName_m(""),
fieldmap_m(NULL),
fast_m(false),
angle_m(0.0),
aperture_m(0.0),
designEnergy_m(0.0),
designRadius_m(0.0),
fieldAmplitude_m(0.0),
bX_m(0.0),
bY_m(0.0),
angleGreaterThanPi_m(false),
entranceAngle_m(0.0),
exitAngle_m(0.0),
fieldIndex_m(0.0),
elementEdge_m(0.0),
startField_m(0.0),
endField_m(0.0),
reinitialize_m(false),
recalcRefTraj_m(false),
length_m(0.0),
gap_m(0.0),
refTrajMapSize_m(0),
refTrajMapStepSize_m(0.0),
entranceParameter1_m(0.0),
entranceParameter2_m(0.0),
entranceParameter3_m(0.0),
exitParameter1_m(0.0),
exitParameter2_m(0.0),
exitParameter3_m(0.0),
xOriginEngeEntry_m(0.0),
zOriginEngeEntry_m(0.0),
deltaBeginEntry_m(0.0),
deltaEndEntry_m(0.0),
polyOrderEntry_m(0),
xExit_m(0.0),
zExit_m(0.0),
xOriginEngeExit_m(0.0),
zOriginEngeExit_m(0.0),
deltaBeginExit_m(0.0),
deltaEndExit_m(0.0),
polyOrderExit_m(0),
cosEntranceAngle_m(1.0),
sinEntranceAngle_m(0.0),
exitEdgeAngle_m(0.0),
cosExitAngle_m(1.0),
sinExitAngle_m(0.0) {
setElType(isDipole);
}
Bend::~Bend() {
}
/*
* OPAL-T Methods.
* ===============
*/
/*
* This function merely repackages the field arrays as type Vector_t and calls
* the equivalent method but with the Vector_t data types.
*/
bool Bend::apply(const size_t &i, const double &t, double E[], double B[]) {
Vector_t Ev(0.0, 0.0, 0.0);
Vector_t Bv(0.0, 0.0, 0.0);
if(apply(RefPartBunch_m->R[i], RefPartBunch_m->get_rmean(), 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 Bend::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
if(designRadius_m > 0.0) {
// Shift position to magnet frame.
Vector_t X = RefPartBunch_m->X[i];
X(2) += startField_m - elementEdge_m;
/*
* Add in transverse bend displacements. (ds is already
* accounted for.)
*/
X(0) -= dx_m;
X(1) -= dy_m;
// Get field from field map.
Vector_t eField(0.0, 0.0, 0.0);
Vector_t bField(0.0, 0.0, 0.0);
CalculateMapField(X, eField, bField);
bField *= fieldAmplitude_m;
B(0) += bField(0);
B(1) += bField(1);
B(2) += bField(2);
}
return false;
}
bool Bend::apply(const Vector_t &R,
const Vector_t &centroid,
const double &t,
Vector_t &E,
Vector_t &B) {
if(designRadius_m > 0.0) {
int index = static_cast<int>
(std::floor((R(2) - startField_m) / refTrajMapStepSize_m));
if(index > 0 && index + 1 < refTrajMapSize_m) {
// Find indices for position in pre-computed central trajectory map.
double lever = (R(2) - startField_m) / refTrajMapStepSize_m - index;
double x = (1.0 - lever) * refTrajMapX_m.at(index)
+ lever * refTrajMapX_m.at(index + 1);
double y = (1.0 - lever) * refTrajMapY_m.at(index)
+ lever * refTrajMapY_m.at(index + 1);
double z = (1.0 - lever) * refTrajMapZ_m.at(index)
+ lever * refTrajMapZ_m.at(index + 1);
// Adjust position relative to pre-computed central trajectory map.
Vector_t X(0.0, 0.0, 0.0);
X(0) = R(0) + x;
X(1) = R(1) + y;
X(2) = z;
Vector_t tempE(0.0, 0.0, 0.0);
Vector_t tempB(0.0, 0.0, 0.0);
Vector_t XInBendFrame = RotateToBendFrame(X);
/*
* Add in transverse bend displacements. (ds is already
* accounted for.)
*/
XInBendFrame(0) -= dx_m;
XInBendFrame(1) -= dy_m;
CalculateMapField(XInBendFrame, tempE, tempB);
tempB = fieldAmplitude_m * RotateOutOfBendFrame(tempB);
B(0) += tempB(0);
B(1) += tempB(1);
B(2) += tempB(2);
}
}
return false;
}
bool Bend::bends() const {
return true;
}
void Bend::finalise() {
online_m = false;
}
void Bend::goOnline(const double &) {
// Check if we need to reinitialize the bend field amplitude.
if(reinitialize_m) {
reinitialize_m = Reinitialize();
recalcRefTraj_m = false;
}
/*
* Always recalculate the reference trajectory on first call even
* if we do not reinitialize the bend. The reference trajectory
* has to be calculated at the same energy as the actual beam or
* we do not get accurate values for the magnetic field in the output
* file.
*/
if(recalcRefTraj_m) {
double angleX = 0.0;
double angleY = 0.0;
CalculateRefTrajectory(angleX, angleY);
recalcRefTraj_m = false;
}
online_m = true;
}
void Bend::getDimensions(double &sBegin, double &sEnd) const {
sBegin = startField_m;
sEnd = endField_m;
}
void Bend::initialise(PartBunch *bunch,
double &startField,
double &endField,
const double &scaleFactor) {
Inform msg(messageHeader_m.c_str(), *gmsg);
if(InitializeFieldMap(msg)) {
retrieveDesignEnergy(startField);
SetupPusher(bunch);
ReadFieldMap(msg);
SetupBendGeometry(msg, startField, endField);
double bendAngleX = 0.0;
double bendAngleY = 0.0;
CalculateRefTrajectory(bendAngleX, bendAngleY);
recalcRefTraj_m = true;
Print(msg, bendAngleX, bendAngleY);
// Pass start and end of field to calling function.
startField = startField_m;
endField = endField_m;
} else {
ERRORMSG("There is something wrong with your field map \""
<< fileName_m
<< "\"");
endField = startField - 1.0e-3;
}
}
double Bend::GetBendAngle() const {
return angle_m;
}
double Bend::GetBendRadius() const {
return designRadius_m;
}
double Bend::GetEffectiveCenter() const {
return elementEdge_m + designRadius_m * angle_m / 2.0;
}
double Bend::GetEffectiveLength() const {
return designRadius_m * angle_m;
}
std::string Bend::GetFieldMapFN() const {
return fileName_m;
}
double Bend::GetStartElement() const {
return elementEdge_m;
}
void Bend::SetAngleGreaterThanPiFlag(bool angleGreaterThanPi) {
angleGreaterThanPi_m = angleGreaterThanPi;
}
void Bend::SetAperture(double aperture) {
aperture_m = std::abs(aperture);
}
void Bend::SetBendAngle(double angle) {
angle_m = angle;
}
void Bend::SetBeta(double beta) {
Orientation_m(1) = beta;
}
void Bend::SetDesignEnergy(double energy) {
designEnergy_m = std::abs(energy);
}
void Bend::SetEntranceAngle(double entranceAngle) {
entranceAngle_m = entranceAngle;
}
void Bend::setExitAngle(double exitAngle) {
exitAngle_m = exitAngle;
}
void Bend::SetFieldAmplitude(double k0, double k0s) {
bY_m = k0;
bX_m = k0s;
}
void Bend::SetFieldMapFN(std::string fileName) {
fileName_m = fileName;
}
void Bend::SetFullGap(double gap) {
gap_m = std::abs(gap);
}
void Bend::SetK1(double k1) {
fieldIndex_m = k1;
}
void Bend::SetLength(double length) {
length_m = std::abs(length);
}
void Bend::SetRotationAboutZ(double rotation) {
Orientation_m(2) = rotation;
}
void Bend::AdjustFringeFields(double ratio) {
double delta = std::abs(entranceParameter1_m - entranceParameter2_m);
entranceParameter1_m = entranceParameter2_m - delta * ratio;
delta = std::abs(entranceParameter2_m - entranceParameter3_m);
entranceParameter3_m = entranceParameter2_m + delta * ratio;
delta = std::abs(exitParameter1_m - exitParameter2_m);
exitParameter1_m = exitParameter2_m - delta * ratio;
delta = std::abs(exitParameter2_m - exitParameter3_m);
exitParameter3_m = exitParameter2_m + delta * ratio;
}
double Bend::CalculateBendAngle() {
const double mass = RefPartBunch_m->getM();
const double gamma = designEnergy_m / mass + 1.0;
const double betaGamma = sqrt(pow(gamma, 2.0) - 1.0);
const double beta = betaGamma / gamma;
const double deltaT = RefPartBunch_m->getdT();
// Integrate through field for initial angle.
Vector_t X(0.0, 0.0, startField_m - elementEdge_m);
Vector_t P(0.0, 0.0, betaGamma);
double deltaS = 0.0;
double bendLength = endField_m - startField_m;
while(deltaS < bendLength) {
X /= Vector_t(Physics::c * deltaT);
pusher_m.push(X, P, deltaT);
X *= Vector_t(Physics::c * deltaT);
Vector_t eField(0.0, 0.0, 0.0);
Vector_t bField(0.0, 0.0, 0.0);
CalculateMapField(X, eField, bField);
bField = fieldAmplitude_m * bField;
X /= Vector_t(Physics::c * deltaT);
pusher_m.kick(X, P, eField, bField, deltaT);
pusher_m.push(X, P, deltaT);
X *= Vector_t(Physics::c * deltaT);
deltaS += deltaT * beta * Physics::c;
}
double angle = -atan2(P(0), P(2));
return angle;
}
void Bend::CalcCentralField(const Vector_t &R,
double deltaX,
double angle,
Vector_t &B) {
double nOverRho = fieldIndex_m / designRadius_m;
double expFactor = exp(-nOverRho * deltaX);
double bxBzFactor = expFactor * nOverRho * R(1);
B(0) = -bxBzFactor * cos(angle);
B(1) = expFactor * (1.0 - pow(nOverRho * R(1), 2.0) / 2.0);
B(2) = -bxBzFactor * sin(angle);
}
void Bend::CalcEngeFunction(double zNormalized,
const std::vector<double> &engeCoeff,
int polyOrder,
double &engeFunc,
double &engeFuncDeriv,
double &engeFuncSecDerivNorm) {
double expSum = 0.0;
double expSumDeriv = 0.0;
double expSumSecDeriv = 0.0;
if(polyOrder >= 2) {
expSum = engeCoeff.at(0)
+ engeCoeff.at(1) * zNormalized;
expSumDeriv = engeCoeff.at(1);
for(int index = 2; index <= polyOrder; index++) {
expSum += engeCoeff.at(index) * pow(zNormalized, index);
expSumDeriv += index * engeCoeff.at(index)
* pow(zNormalized, index - 1);
expSumSecDeriv += index * (index - 1) * engeCoeff.at(index)
* pow(zNormalized, index - 2);
}
} else if(polyOrder == 1) {
expSum = engeCoeff.at(0)
+ engeCoeff.at(1) * zNormalized;
expSumDeriv = engeCoeff.at(1);
} else
expSum = engeCoeff.at(0);
double engeExp = exp(expSum);
engeFunc = 1.0 / (1.0 + engeExp);
if(!std::isnan(engeFunc)) {
expSumDeriv /= gap_m;
expSumSecDeriv /= pow(gap_m, 2.0);
double engeExpDeriv = expSumDeriv * engeExp;
double engeExpSecDeriv = (expSumSecDeriv + pow(expSumDeriv, 2.0))
* engeExp;
double engeFuncSq = pow(engeFunc, 2.0);
engeFuncDeriv = -engeExpDeriv * engeFuncSq;
if (std::isnan(engeFuncDeriv) || std::isinf(engeFuncDeriv))
engeFuncDeriv = 0.0;
engeFuncSecDerivNorm = -engeExpSecDeriv * engeFunc
+ 2.0 * pow(engeExpDeriv, 2.0)
* engeFuncSq;
if (std::isnan(engeFuncSecDerivNorm) || std::isinf(engeFuncSecDerivNorm))
engeFuncSecDerivNorm = 0.0;
} else {
engeFunc = 0.0;
engeFuncDeriv = 0.0;
engeFuncSecDerivNorm = 0.0;
}
}
void Bend::CalcEntranceFringeField(const Vector_t &REntrance,
double deltaX,
Vector_t &B) {
double zNormalized = -REntrance(2) / gap_m;
double engeFunc = 0.0;
double engeFuncDeriv = 0.0;
double engeFuncSecDerivNorm = 0.0;
CalcEngeFunction(zNormalized,
engeCoeffsEntry_m,
polyOrderEntry_m,
engeFunc,
engeFuncDeriv,
engeFuncSecDerivNorm);
double nOverRho = fieldIndex_m / designRadius_m;
double expFactor = exp(-nOverRho * deltaX);
double trigFactor = pow(nOverRho, 2.0) + engeFuncSecDerivNorm;
double bXEntrance = -engeFunc * nOverRho * expFactor* REntrance(1);
double bYEntrance = expFactor * engeFunc
* (1.0 - trigFactor * pow(REntrance(1), 2.0) / 2.0);
double bZEntrance = -expFactor * engeFuncDeriv * REntrance(1);
B(0) = bXEntrance * cosEntranceAngle_m - bZEntrance * sinEntranceAngle_m;
B(1) = bYEntrance;
B(2) = bXEntrance * sinEntranceAngle_m + bZEntrance * cosEntranceAngle_m;
}
void Bend::CalcExitFringeField(const Vector_t &RExit, double deltaX, Vector_t &B) {
double zNormalized = RExit(2) / gap_m;
double engeFunc = 0.0;
double engeFuncDeriv = 0.0;
double engeFuncSecDerivNorm = 0.0;
CalcEngeFunction(zNormalized,
engeCoeffsExit_m,
polyOrderExit_m,
engeFunc,
engeFuncDeriv,
engeFuncSecDerivNorm);
double nOverRho = fieldIndex_m / designRadius_m;
double expFactor = exp(-nOverRho * deltaX);
double trigFactor = pow(nOverRho, 2.0) + engeFuncSecDerivNorm;
double bXExit = -engeFunc * nOverRho * expFactor* RExit(1);
double bYExit = expFactor * engeFunc
* (1.0 - trigFactor * pow(RExit(1), 2.0) / 2.0);
double bZExit = expFactor * engeFuncDeriv * RExit(1);
B(0) = bXExit * cosExitAngle_m - bZExit * sinExitAngle_m;
B(1) = bYExit;
B(2) = bXExit * sinExitAngle_m + bZExit * cosExitAngle_m;
}
void Bend::CalculateMapField(const Vector_t &R, Vector_t &E, Vector_t &B) {
E = Vector_t(0.0);
B = Vector_t(0.0);
// Vector_t REntrance(0.0, 0.0, 0.0);
// Vector_t RExit(0.0, 0.0, 0.0);
// if (IsPositionInEntranceField(R, REntrance)) {
// CalcEntranceFringeField(REntrance, 0.0, B);
// } else if (IsPositionInExitField(R, RExit)) {
// CalcExitFringeField(RExit, 0.0, B);
// } else {
// CalcCentralField(R, 0.0, 0.0, B);
// }
double deltaXEntrance = 0.0;
double deltaXExit = 0.0;
bool inEntranceRegion = InMagnetEntranceRegion(R, deltaXEntrance);