Commit f9be6c67 authored by kraus's avatar kraus

Merge branch 'issue515' into 'OPAL-2.2'

Change origin of the position parameter in the getFieldStrength method in OPAL-2.2

See merge request !337
parents 1133f86f 6d76e398
......@@ -63,7 +63,6 @@ RFCavity::RFCavity(const RFCavity &right):
fieldmap_m(right.fieldmap_m),
length_m(right.length_m),
startField_m(right.startField_m),
endField_m(right.endField_m),
type_m(right.type_m),
rmin_m(right.rmin_m),
rmax_m(right.rmax_m),
......@@ -99,7 +98,6 @@ RFCavity::RFCavity(const std::string &name):
fieldmap_m(nullptr),
length_m(0.0),
startField_m(0.0),
endField_m(0.0),
type_m(SW),
rmin_m(0.0),
rmax_m(0.0),
......@@ -203,12 +201,11 @@ bool RFCavity::apply(const Vector_t &R,
const double &t,
Vector_t &E,
Vector_t &B) {
Vector_t tmpR(R(0), R(1), R(2) - startField_m);
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
if (R(2) >= startField_m &&
R(2) < startField_m + length_m) {
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
if (tmpR(2) >= 0.0 &&
tmpR(2) < length_m) {
bool outOfBounds = fieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
bool outOfBounds = fieldmap_m->getFieldstrength(R, tmpE, tmpB);
if (outOfBounds) return true;
E += (scale_m + scaleError_m) * cos(frequency_m * t + phase_m + phaseError_m) * tmpE;
......@@ -224,12 +221,11 @@ bool RFCavity::applyToReferenceParticle(const Vector_t &R,
Vector_t &E,
Vector_t &B) {
Vector_t tmpR(R(0), R(1), R(2) - startField_m);
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
if (R(2) >= startField_m &&
R(2) < startField_m + length_m) {
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
if (tmpR(2) >= 0.0 &&
tmpR(2) < length_m) {
bool outOfBounds = fieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
bool outOfBounds = fieldmap_m->getFieldstrength(R, tmpE, tmpB);
if (outOfBounds) return true;
E += scale_m * cos(frequency_m * t + phase_m) * tmpE;
......@@ -242,10 +238,10 @@ bool RFCavity::applyToReferenceParticle(const Vector_t &R,
void RFCavity::initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField) {
using Physics::two_pi;
startField_m = endField_m = 0.0;
startField_m = 0.0;
if (bunch == NULL) {
startField = startField_m;
endField = endField_m;
endField = startField_m;
return;
}
......@@ -256,8 +252,8 @@ void RFCavity::initialise(PartBunchBase<double, 3> *bunch, double &startField, d
RefPartBunch_m = bunch;
fieldmap_m = Fieldmap::getFieldmap(filename_m, fast_m);
fieldmap_m->getFieldDimensions(startField_m, endField_m, rBegin, rEnd);
if (endField_m <= startField_m) {
fieldmap_m->getFieldDimensions(startField_m, endField, rBegin, rEnd);
if (endField <= startField_m) {
throw GeneralClassicException("RFCavity::initialise",
"The length of the field map '" + filename_m + "' is zero or negativ");
}
......@@ -277,9 +273,7 @@ void RFCavity::initialise(PartBunchBase<double, 3> *bunch, double &startField, d
}
frequency_m = fieldmap_m->getFrequency();
}
length_m = endField_m - startField_m;
endField = startField + length_m;
length_m = endField - startField_m;
}
// In current version ,this function reads in the cavity voltage profile data from file.
......@@ -570,7 +564,7 @@ double RFCavity::spline(double z, double *za) {
void RFCavity::getDimensions(double &zBegin, double &zEnd) const {
zBegin = startField_m;
zEnd = endField_m;
zEnd = startField_m + length_m;
}
......@@ -822,4 +816,4 @@ void RFCavity::getElementDimensions(double &begin,
double &end) const {
double tmp;
fieldmap_m->getFieldDimensions(begin, end, tmp, tmp);
}
\ No newline at end of file
}
......@@ -234,8 +234,6 @@ protected:
double startField_m; /**< starting point of field(m)*/
private:
double endField_m;
CavityType type_m;
double rmin_m;
......@@ -519,7 +517,7 @@ CoordinateSystemTrafo RFCavity::getEdgeToBegin() const
inline
CoordinateSystemTrafo RFCavity::getEdgeToEnd() const
{
CoordinateSystemTrafo ret(Vector_t(0, 0, endField_m),
CoordinateSystemTrafo ret(Vector_t(0, 0, startField_m + length_m),
Quaternion(1, 0, 0, 0));
return ret;
......
......@@ -162,13 +162,13 @@ bool Solenoid::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);
}
bool Solenoid::apply(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B) {
const Vector_t tmpR(R(0), R(1), R(2) - startField_m);
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
bool Solenoid::apply(const Vector_t &R, const Vector_t &/*P*/, const double &/*t*/, Vector_t &/*E*/, Vector_t &B) {
if (R(2) >= startField_m
&& R(2) < startField_m + length_m) {
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
if (tmpR(2) >= 0.0 && tmpR(2) < length_m) {
const bool out_of_bounds = myFieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
if (out_of_bounds) return true;
const bool outOfBounds = myFieldmap_m->getFieldstrength(R, tmpE, tmpB);
if (outOfBounds) return true;
B += (scale_m + scaleError_m) * tmpB;
}
......@@ -176,13 +176,13 @@ bool Solenoid::apply(const Vector_t &R, const Vector_t &P, const double &t, Vec
return false;
}
bool Solenoid::applyToReferenceParticle(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B) {
const Vector_t tmpR(R(0), R(1), R(2) - startField_m);
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
bool Solenoid::applyToReferenceParticle(const Vector_t &R, const Vector_t &/*P*/, const double &/*t*/, Vector_t &/*E*/, Vector_t &B) {
if (R(2) >= startField_m
&& R(2) < startField_m + length_m) {
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
if (tmpR(2) >= 0.0 && tmpR(2) < length_m) {
const bool out_of_bounds = myFieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
if (out_of_bounds) return true;
const bool outOfBounds = myFieldmap_m->getFieldstrength(R, tmpE, tmpB);
if (outOfBounds) return true;
B += scale_m * tmpB;
}
......@@ -249,7 +249,8 @@ ElementBase::ElementType Solenoid::getType() const {
}
bool Solenoid::isInside(const Vector_t &r) const {
return (r(2) >= startField_m && r(2) < startField_m + length_m && isInsideTransverse(r) && myFieldmap_m->isInside(r));
return isInsideTransverse(r)
&& myFieldmap_m->isInside(r);
}
......@@ -261,4 +262,4 @@ void Solenoid::getElementDimensions(double &begin,
double &end) const {
begin = startField_m;
end = begin + length_m;
}
\ No newline at end of file
}
......@@ -326,8 +326,12 @@ void TravelingWave::initialise(PartBunchBase<double, 3> *bunch, double &startFie
RefPartBunch_m = bunch;
double zBegin = 0.0, zEnd = 0.0;
RFCavity::initialise(bunch, zBegin, zEnd);
if (std::abs(startField_m) > 0.0) {
throw GeneralClassicException("TravelingWave::initialise",
"The field map of a traveling wave structure has to begin at 0.0");
}
PeriodLength_m = (zEnd - zBegin) / 2.0;
PeriodLength_m = length_m / 2.0;
CellLength_m = PeriodLength_m * Mode_m;
startField_m = -0.5 * PeriodLength_m;
......@@ -520,5 +524,7 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
}
bool TravelingWave::isInside(const Vector_t &r) const {
return (isInsideTransverse(r) && r(2) >= -0.5 * PeriodLength_m && r(2) < startExitField_m);
return (isInsideTransverse(r)
&& r(2) >= -0.5 * PeriodLength_m
&& r(2) < startExitField_m);
}
......@@ -191,7 +191,7 @@ bool Astra1DDynamic::getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &
// do fourier interpolation in z-direction
const double RR2 = R(0) * R(0) + R(1) * R(1);
const double kz = two_pi * R(2) / length_m + Physics::pi;
const double kz = two_pi * (R(2) - zbegin_m) / length_m + Physics::pi;
double ez = FourCoefs_m[0];
double ezp = 0.0;
......@@ -231,8 +231,8 @@ bool Astra1DDynamic::getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &
return false;
}
bool Astra1DDynamic::getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const {
const double kz = two_pi * R(2) / length_m + Physics::pi;
bool Astra1DDynamic::getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &/*B*/, const DiffDirection &/*dir*/) const {
const double kz = two_pi * (R(2) - zbegin_m) / length_m + Physics::pi;
double ezp = 0.0;
int n = 1;
......@@ -288,4 +288,4 @@ void Astra1DDynamic::getOnaxisEz(vector<pair<double, double> > & F) {
for (int i = 0; i < num_gridpz_m; ++ i) {
F[i].second /= Ez_max;
}
}
\ No newline at end of file
}
......@@ -82,10 +82,10 @@ bool Astra1DDynamic_fast::getFieldstrength(const Vector_t &R, Vector_t &E, Vecto
// do fourier interpolation in z-direction
const double RR2 = R(0) * R(0) + R(1) * R(1);
double ez = gsl_spline_eval(onAxisInterpolants_m[0], R(2), onAxisAccel_m[0]);
double ezp = gsl_spline_eval(onAxisInterpolants_m[1], R(2), onAxisAccel_m[1]);
double ezpp = gsl_spline_eval(onAxisInterpolants_m[2], R(2), onAxisAccel_m[2]);
double ezppp = gsl_spline_eval(onAxisInterpolants_m[3], R(2), onAxisAccel_m[3]);
double ez = gsl_spline_eval(onAxisInterpolants_m[0], R(2) - zbegin_m, onAxisAccel_m[0]);
double ezp = gsl_spline_eval(onAxisInterpolants_m[1], R(2) - zbegin_m, onAxisAccel_m[1]);
double ezpp = gsl_spline_eval(onAxisInterpolants_m[2], R(2) - zbegin_m, onAxisAccel_m[2]);
double ezppp = gsl_spline_eval(onAxisInterpolants_m[3], R(2) - zbegin_m, onAxisAccel_m[3]);
// expand the field off-axis
const double f = -(ezpp + ez * xlrep_m * xlrep_m) / 16.;
......@@ -103,8 +103,8 @@ bool Astra1DDynamic_fast::getFieldstrength(const Vector_t &R, Vector_t &E, Vecto
return false;
}
bool Astra1DDynamic_fast::getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const {
double ezp = gsl_spline_eval(onAxisInterpolants_m[1], R(2), onAxisAccel_m[1]);
bool Astra1DDynamic_fast::getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &/*B*/, const DiffDirection &/*dir*/) const {
double ezp = gsl_spline_eval(onAxisInterpolants_m[1], R(2) - zbegin_m, onAxisAccel_m[1]);
E(2) += ezp;
......@@ -149,8 +149,8 @@ void Astra1DDynamic_fast::getOnaxisEz(std::vector<std::pair<double, double> > &
for(int i = 0; i < num_gridpz_m; ++ i) {
F[i].first = hz_m * i;
interpretLine<double>(in, F[i].second);
if(fabs(F[i].second) > Ez_max) {
Ez_max = fabs(F[i].second);
if(std::abs(F[i].second) > Ez_max) {
Ez_max = std::abs(F[i].second);
}
}
in.close();
......@@ -205,4 +205,4 @@ int Astra1DDynamic_fast::stripFileHeader(std::ifstream &file) {
interpretLine<double>(file, tmpDouble);
return accuracy;
}
\ No newline at end of file
}
......@@ -184,7 +184,7 @@ bool Astra1DElectroStatic::getFieldstrength(const Vector_t &R, Vector_t &E, Vect
// do fourier interpolation in z-direction
const double RR2 = R(0) * R(0) + R(1) * R(1);
const double kz = two_pi * (R(2) - length_m) / length_m;
const double kz = two_pi * (R(2) - zbegin_m) / length_m + Physics::pi;
double ez = FourCoefs_m[0];
double ezp = 0.0;
......
......@@ -76,10 +76,10 @@ bool Astra1DElectroStatic_fast::getFieldstrength(const Vector_t &R, Vector_t &E,
// do fourier interpolation in z-direction
const double RR2 = R(0) * R(0) + R(1) * R(1);
double ez = gsl_spline_eval(onAxisInterpolants_m[0], R(2), onAxisAccel_m[0]);
double ezp = gsl_spline_eval(onAxisInterpolants_m[1], R(2), onAxisAccel_m[1]);
double ezpp = gsl_spline_eval(onAxisInterpolants_m[2], R(2), onAxisAccel_m[2]);
double ezppp = gsl_spline_eval(onAxisInterpolants_m[3], R(2), onAxisAccel_m[3]);
double ez = gsl_spline_eval(onAxisInterpolants_m[0], R(2) - zbegin_m, onAxisAccel_m[0]);
double ezp = gsl_spline_eval(onAxisInterpolants_m[1], R(2) - zbegin_m, onAxisAccel_m[1]);
double ezpp = gsl_spline_eval(onAxisInterpolants_m[2], R(2) - zbegin_m, onAxisAccel_m[2]);
double ezppp = gsl_spline_eval(onAxisInterpolants_m[3], R(2) - zbegin_m, onAxisAccel_m[3]);
// expand to off-axis
const double EfieldR = -ezp / 2. + ezppp / 16. * RR2;
......
......@@ -113,8 +113,8 @@ void Astra1DMagnetoStatic::readMap() {
// the sequence of z-position should be strictly increasing
// drop sampling points that don't comply to this
if(zvals[i] - tmpDouble > 1e-10) {
if(fabs(RealValues[i]) > Bz_max) {
Bz_max = fabs(RealValues[i]);
if(std::abs(RealValues[i]) > Bz_max) {
Bz_max = std::abs(RealValues[i]);
}
tmpDouble = zvals[i];
++ i; // increment i only if sampling point is accepted
......@@ -176,7 +176,7 @@ bool Astra1DMagnetoStatic::getFieldstrength(const Vector_t &R, Vector_t &E, Vect
// do fourier interpolation in z-direction
const double RR2 = R(0) * R(0) + R(1) * R(1);
const double kz = two_pi * R(2) / length_m + Physics::pi;
const double kz = two_pi * (R(2) - zbegin_m) / length_m + Physics::pi;
double ez = FourCoefs_m[0];
double ezp = 0.0;
......
......@@ -76,10 +76,10 @@ bool Astra1DMagnetoStatic_fast::getFieldstrength(const Vector_t &R, Vector_t &E,
// do fourier interpolation in z-direction
const double RR2 = R(0) * R(0) + R(1) * R(1);
double bz = gsl_spline_eval(onAxisInterpolants_m[0], R(2), onAxisAccel_m[0]);
double bzp = gsl_spline_eval(onAxisInterpolants_m[1], R(2), onAxisAccel_m[1]);
double bzpp = gsl_spline_eval(onAxisInterpolants_m[2], R(2), onAxisAccel_m[2]);
double bzppp = gsl_spline_eval(onAxisInterpolants_m[3], R(2), onAxisAccel_m[3]);
double bz = gsl_spline_eval(onAxisInterpolants_m[0], R(2) - zbegin_m, onAxisAccel_m[0]);
double bzp = gsl_spline_eval(onAxisInterpolants_m[1], R(2) - zbegin_m, onAxisAccel_m[1]);
double bzpp = gsl_spline_eval(onAxisInterpolants_m[2], R(2) - zbegin_m, onAxisAccel_m[2]);
double bzppp = gsl_spline_eval(onAxisInterpolants_m[3], R(2) - zbegin_m, onAxisAccel_m[3]);
// expand to off-axis
const double BfieldR = -bzp / 2. + bzppp / 16. * RR2;
......
......@@ -72,7 +72,7 @@ bool FM1DDynamic::getFieldstrength(const Vector_t &R, Vector_t &E,
Vector_t &B) const {
std::vector<double> fieldComponents;
computeFieldOnAxis(R(2), fieldComponents);
computeFieldOnAxis(R(2) - zBegin_m, fieldComponents);
computeFieldOffAxis(R, E, B, fieldComponents);
return false;
......@@ -83,7 +83,7 @@ bool FM1DDynamic::getFieldDerivative(const Vector_t &R,
Vector_t &B,
const DiffDirection &dir) const {
double kz = Physics::two_pi * R(2) / length_m + Physics::pi;
double kz = Physics::two_pi * (R(2) - zBegin_m) / length_m + Physics::pi;
double eZPrime = 0.0;
int coefIndex = 1;
......
......@@ -100,7 +100,7 @@ bool FM1DDynamic_fast::getFieldstrength(const Vector_t &R, Vector_t &E,
Vector_t &B) const {
std::vector<double> fieldComponents;
computeFieldOnAxis(R(2), fieldComponents);
computeFieldOnAxis(R(2) - zBegin_m, fieldComponents);
computeFieldOffAxis(R, E, B, fieldComponents);
return false;
......@@ -111,7 +111,7 @@ bool FM1DDynamic_fast::getFieldDerivative(const Vector_t &R,
Vector_t &B,
const DiffDirection &dir) const {
E(2) += gsl_spline_eval(onAxisFieldPInterpolants_m, R(2),
E(2) += gsl_spline_eval(onAxisFieldPInterpolants_m, R(2) - zBegin_m,
onAxisFieldPAccel_m);
return false;
......
......@@ -72,7 +72,7 @@ bool FM1DElectroStatic::getFieldstrength(const Vector_t &R, Vector_t &E,
Vector_t &B) const {
std::vector<double> fieldComponents;
computeFieldOnAxis(R(2), fieldComponents);
computeFieldOnAxis(R(2) - zBegin_m, fieldComponents);
computeFieldOffAxis(R, E, B, fieldComponents);
return false;
......@@ -83,7 +83,7 @@ bool FM1DElectroStatic::getFieldDerivative(const Vector_t &R,
Vector_t &B,
const DiffDirection &dir) const {
double kz = Physics::two_pi * R(2) / length_m + Physics::pi;
double kz = Physics::two_pi * (R(2) - zBegin_m) / length_m + Physics::pi;
double eZPrime = 0.0;
int coefIndex = 1;
......
......@@ -100,7 +100,7 @@ bool FM1DElectroStatic_fast::getFieldstrength(const Vector_t &R, Vector_t &E,
Vector_t &B) const {
std::vector<double> fieldComponents;
computeFieldOnAxis(R(2), fieldComponents);
computeFieldOnAxis(R(2) - zBegin_m, fieldComponents);
computeFieldOffAxis(R, E, B, fieldComponents);
return false;
......@@ -112,7 +112,7 @@ bool FM1DElectroStatic_fast::getFieldDerivative(const Vector_t &R,
Vector_t &B,
const DiffDirection &dir) const {
E(2) += gsl_spline_eval(onAxisFieldPInterpolants_m, R(2),
E(2) += gsl_spline_eval(onAxisFieldPInterpolants_m, R(2) - zBegin_m,
onAxisFieldPAccel_m);
return false;
......
......@@ -72,7 +72,7 @@ bool FM1DMagnetoStatic::getFieldstrength(const Vector_t &R, Vector_t &E,
Vector_t &B) const {
std::vector<double> fieldComponents;
computeFieldOnAxis(R(2), fieldComponents);
computeFieldOnAxis(R(2) - zBegin_m, fieldComponents);
computeFieldOffAxis(R, E, B, fieldComponents);
return false;
......@@ -83,7 +83,7 @@ bool FM1DMagnetoStatic::getFieldDerivative(const Vector_t &R,
Vector_t &B,
const DiffDirection &dir) const {
double kz = Physics::two_pi * R(2) / length_m + Physics::pi;
double kz = Physics::two_pi * (R(2) - zBegin_m) / length_m + Physics::pi;
double bZPrime = 0.0;
int coefIndex = 1;
......
......@@ -99,7 +99,7 @@ bool FM1DMagnetoStatic_fast::getFieldstrength(const Vector_t &R, Vector_t &E,
Vector_t &B) const {
std::vector<double> fieldComponents;
computeFieldOnAxis(R(2), fieldComponents);
computeFieldOnAxis(R(2) - zBegin_m, fieldComponents);
computeFieldOffAxis(R, E, B, fieldComponents);
return false;
......@@ -111,7 +111,7 @@ bool FM1DMagnetoStatic_fast::getFieldDerivative(const Vector_t &R,
Vector_t &B,
const DiffDirection &dir) const {
B(2) += gsl_spline_eval(onAxisFieldPInterpolants_m, R(2),
B(2) += gsl_spline_eval(onAxisFieldPInterpolants_m, R(2) - zBegin_m,
onAxisFieldPAccel_m);
return false;
......
......@@ -193,11 +193,11 @@ bool FM2DDynamic::getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B)
// do bi-linear interpolation
const double RR = sqrt(R(0) * R(0) + R(1) * R(1));
const int indexr = (int)floor(RR / hr_m);
const int indexr = (int)std::floor(RR / hr_m);
const double leverr = RR / hr_m - indexr;
const int indexz = (int)floor(R(2) / hz_m);
const double leverz = R(2) / hz_m - indexz;
const int indexz = (int)std::floor((R(2) - zbegin_m) / hz_m);
const double leverz = (R(2) - zbegin_m) / hz_m - indexz;
if((indexz < 0) || (indexz + 2 > num_gridpz_m))
return false;
......
......@@ -121,7 +121,7 @@ void FM2DElectroStatic::readMap() {
FieldstrengthEr_m[i + j * num_gridpz_m],
FieldstrengthEz_m[i + j * num_gridpz_m]);
}
if (fabs(FieldstrengthEz_m[i]) > Ezmax) Ezmax = fabs(FieldstrengthEz_m[i]);
if (std::abs(FieldstrengthEz_m[i]) > Ezmax) Ezmax = std::abs(FieldstrengthEz_m[i]);
}
} else {
for (int j = 0; j < num_gridpr_m; ++ j) {
......@@ -169,11 +169,11 @@ bool FM2DElectroStatic::getFieldstrength(const Vector_t &R, Vector_t &E, Vector_
// do bi-linear interpolation
const double RR = sqrt(R(0) * R(0) + R(1) * R(1));
const int indexr = (int)floor(RR / hr_m);
const int indexr = (int)std::floor(RR / hr_m);
const double leverr = (RR / hr_m) - indexr;
const int indexz = (int)floor((R(2)) / hz_m);
const double leverz = (R(2) / hz_m) - indexz;
const int indexz = (int)std::floor((R(2) - zbegin_m) / hz_m);
const double leverz = (R(2) - zbegin_m) / hz_m - indexz;
if ((indexz < 0) || (indexz + 2 > num_gridpz_m))
return false;
......
......@@ -169,11 +169,11 @@ bool FM2DMagnetoStatic::getFieldstrength(const Vector_t &R, Vector_t &E, Vector_
// do bi-linear interpolation
const double RR = sqrt(R(0) * R(0) + R(1) * R(1));
const int indexr = abs((int)floor(RR / hr_m));
const int indexr = std::abs((int)std::floor(RR / hr_m));
const double leverr = (RR / hr_m) - indexr;
const int indexz = abs((int)floor((R(2)) / hz_m));
const double leverz = (R(2) / hz_m) - indexz;
const int indexz = std::abs((int)std::floor((R(2) - zbegin_m) / hz_m));
const double leverz = (R(2) - zbegin_m) / hz_m - indexz;
if((indexz < 0) || (indexz + 2 > num_gridpz_m))
return false;
......@@ -207,10 +207,10 @@ bool FM2DMagnetoStatic::getFieldDerivative(const Vector_t &R, Vector_t &E, Vecto
const double RR = sqrt(R(0) * R(0) + R(1) * R(1));
const int indexr = (int)floor(RR / hr_m);
const int indexr = (int)std::floor(RR / hr_m);
const double leverr = (RR / hr_m) - indexr;
const int indexz = (int)floor((R(2)) / hz_m);
const int indexz = (int)std::floor((R(2)) / hz_m);
const double leverz = (R(2) / hz_m) - indexz;
if((indexz < 0) || (indexz + 2 > num_gridpz_m))
......
......@@ -242,14 +242,14 @@ void FM3DDynamic::freeMap() {
}
bool FM3DDynamic::getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const {
const unsigned int index_x = static_cast<int>(floor((R(0) - xbegin_m) / hx_m));
const unsigned int index_x = static_cast<int>(std::floor((R(0) - xbegin_m) / hx_m));
const double lever_x = (R(0) - xbegin_m) / hx_m - index_x;
const unsigned int index_y = static_cast<int>(floor((R(1) - ybegin_m) / hy_m));
const unsigned int index_y = static_cast<int>(std::floor((R(1) - ybegin_m) / hy_m));
const double lever_y = (R(1) - ybegin_m) / hy_m - index_y;
const unsigned int index_z = (int)floor(R(2) / hz_m);
const double lever_z = R(2) / hz_m - index_z;
const unsigned int index_z = (int)std::floor((R(2) - zbegin_m)/ hz_m);
const double lever_z = (R(2) - zbegin_m) / hz_m - index_z;
if(index_z >= num_gridpz_m - 2) {
return false;
......
This diff is collapsed.
This diff is collapsed.
......@@ -40,6 +40,7 @@ FM3DMagnetoStatic::FM3DMagnetoStatic(std::string aFilename):
normalize_m = (tmpString == "TRUE");
}
parsing_passed = parsing_passed &&
interpretLine<double, double, unsigned int>(file, xbegin_m, xend_m, num_gridpx_m);
parsing_passed = parsing_passed &&
......@@ -61,8 +62,6 @@ FM3DMagnetoStatic::FM3DMagnetoStatic(std::string aFilename):
file.close();
if(!parsing_passed) {
disableFieldmapWarning();
zend_m = zbegin_m - 1e-3;
throw GeneralClassicException("FM3DMagnetoStatic::FM3DMagnetoStatic",
"An error occured when reading the fieldmap '" + Filename_m + "'");
} else {
......@@ -83,9 +82,8 @@ FM3DMagnetoStatic::FM3DMagnetoStatic(std::string aFilename):
}
} else {
noFieldmapWarning();
zbegin_m = 0.0;
zend_m = -1e-3;
throw GeneralClassicException("FM3DMagnetoStatic::FM3DMagnetoStatic",
"An error occured when reading the fieldmap '" + Filename_m + "'");
}
}
......@@ -164,9 +162,12 @@ void FM3DMagnetoStatic::freeMap() {
}
}
bool FM3DMagnetoStatic::getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const {
if (isInside(R))
B += interpolateTrilinearly(R);
bool FM3DMagnetoStatic::getFieldstrength(const Vector_t &R, Vector_t &/*E*/, Vector_t &B) const {
if (!isInside(R)) {
return true;
}
B += interpolateTrilinearly(R);
return false;
}
......@@ -235,4 +236,4 @@ void FM3DMagnetoStatic::getInfo(Inform *msg) {
<< " xini= " << xbegin_m << " xfinal= " << xend_m
<< " yini= " << ybegin_m << " yfinal= " << yend_m
<< " zini= " << zbegin_m << " zfinal= " << zend_m << " (m) " << endl;
}
\ No newline at end of file
}
......@@ -60,7 +60,6 @@ private:
double zbegin_m;
double zend_m;
double length_m;
double hx_m; /**< length between points in grid, x-direction */
double hy_m; /**< length between points in grid, y-direction */
......@@ -75,10 +74,9 @@ private:
inline
bool FM3DMagnetoStatic::isInside(const Vector_t &r) const
{
return ((r(0) >= xbegin_m && r(0) < xend_m) &&
(r(1) >= ybegin_m && r(1) < yend_m) &&
(r(2) >= 0.0 && r(2) < length_m));
// (r(2) >= zbegin_m && r(2) < zend_m));
return ((r(0) >= xbegin_m && r(0) < xend_m)
&& (r(1) >= ybegin_m && r(1) < yend_m)
&& (r(2) >= zbegin_m && r(2) < zend_m));
}
inline
......@@ -103,8 +101,10 @@ FM3DMagnetoStatic::IndexTriplet FM3DMagnetoStatic::getIndex(const Vector_t &X) c
IndexTriplet idx;
idx.i = std::floor((X(0) - xbegin_m) / hx_m);
idx.j = std::floor((X(1) - ybegin_m) / hy_m);
idx.k = std::floor(X(2) / hz_m);
// idx.k = std::floor((X(2) - zbegin_m) / hz_m);
idx.k = std::floor((X(2) - zbegin_m)/ hz_m);
PAssert_LT(idx.i, num_gridpx_m - 1);
PAssert_LT(idx.j, num_gridpy_m - 1);
PAssert_LT(idx.k, num_gridpz_m - 1);
idx.weight(0) = (X(0) - xbegin_m) / hx_m - idx.i;
idx.weight(1) = (X(1) - ybegin_m) / hy_m - idx.j;
......
......@@ -81,7 +81,6 @@ inline
bool FM3DMagnetoStaticExtended::isInside(const Vector_t &r) const
{
return r(2) >= 0.0 && r(2) < length_m && r(0) >= xbegin_m && r(0) < xend_m && std::abs(r(1)) < yend_m;
// return r(2) >= zbegin_m && r(2) < zend_m && r(0) >= xbegin_m && r(0) < xend_m && std::abs(r(1)) < yend_m;
}
inline
......@@ -98,14 +97,16 @@ FM3DMagnetoStaticExtended::IndexTriplet FM3DMagnetoStaticExtended::getIndex(cons
IndexTriplet idx;
idx.i = std::floor((X(0) - xbegin_m) / hx_m);
idx.j = std::floor(std::abs(X(1)) / hy_m);
idx.k = std::floor((X(2)) / hz_m);
// idx.k = std::floor((X(2) - zbegin_m) / hz_m);
idx.k = std::floor((X(2) - zbegin_m) / hz_m);
PAssert_LT(idx.i, num_gridpx_m - 1);
PAssert_LT(idx.j, num_gridpy_m - 1);
PAssert_LT(idx.k, num_gridpz_m - 1);
idx.weight(0) = (X(0) - xbegin_m) / hx_m - idx.i;
idx.weight(1) = std::abs(X(1)) / hy_m - idx.j;
idx.weight(2)