Commit 058f85d9 authored by Steve Russell's avatar Steve Russell
Browse files

Fixed two bugs in the 1D field classes, FM1DDynamic, FM1DDynamic_fast,

FM1DElectroStatic, FM1DElectroStatic_fast, FM1DMagnetoStatic and FM1DMagnetoStatic_fast.

The first bug had to do with field normalization. Although the 1D field profile
derivatives were normalized properly, the actual field was not. Instead of
looping through all of the field values, the loop was over the Fourier 
accuracy parameter. This typically had the effect that the on axis field was not
normalized at all. For input fields that were already normalized to 1, this
had no effect. However, if the input field was not then the field expansion
was incorrect.

The second had to do with the Fourier expansion used to calculate the 
derivatives of the on axis field. The code reflects the input
field in order to ensure that it is periodic. For best accuracy, this 
reflection should be about the first point of the input data (z = 0). However,
what was actually implemented was a reflection about -delta / 2, where delta
is the 1D grid spacing. This caused an error that was small, and became
smaller with decreasing delta, but significant. Using the new implementation,
the reconstruction of the field file for the Envelope-Tracker-Phase-1
(FINSS-RGUN.dat) is about two orders of magnitude more accurate than when
using the old Fourier expansion implementation. Even though this error was
small, there was a noticeable difference in the beam envelope when running
the Envelope-Tracker-Phase-1 and ExternalFieldTest regression tests. The
reference files for both of the tests were updated using the new field
classes.

A quick glance shows that the Fourier expansion bug may also be present in the
1D Astra field classes, but I need to look more carefully to be sure.

Finally, all six classes were cleaned up to conform to the new coding style.
parent a4084ea6
......@@ -179,7 +179,7 @@ void RBend::addKR(int i, double t, Vector_t &K) {
DiffDirection zdir(DZ);
myFieldmap->getFieldstrength(tmpA, tmpE, tmpB);
myFieldmap->getFieldstrength_fdiff(tmpA, tmpE_diff, tmpB_diff, zdir);
myFieldmap->getFieldDerivative(tmpA, tmpE_diff, tmpB_diff, zdir);
double g = RefPartBunch_m->getGamma(i);
......@@ -580,7 +580,7 @@ double RBend::calculateRefTrajectory(const double zBegin) {
bool EntryFringe_passed = false;
double PathLengthEntryFringe = 0.0; // in S coordinates. This value is different from zBegin due to the curvature!
if(map_m != NULL) delete map_m;
if(map_m != NULL) delete [] map_m;
map_step_size_m = betagamma / gamma * Physics::c * dt;
map_size_m = static_cast<int>(floor(length_m / 2. * Physics::pi / map_step_size_m));
......
......@@ -206,6 +206,7 @@ private:
double design_energy_m;
double angle_m; // Bend angle of central trajectory particle (radians).
double *map_m;
int map_size_m;
double map_step_size_m;
......
This diff is collapsed.
......@@ -71,7 +71,7 @@ Solenoid::Solenoid(const string &name):
Solenoid::~Solenoid() {
// Fieldmap::deleteFieldmap(filename_m);
// Fieldmap::deleteFieldmap(filename_m);
}
......@@ -133,7 +133,7 @@ void Solenoid::addKT(int i, double t, Vector_t &K) {
myFieldmap_m->getFieldstrength(tmpA, tmpE, tmpB);
// get derivation of B in z-direction
myFieldmap_m->getFieldstrength_fdiff(tmpA, tmpE_diff, tmpB_diff, zdir);
myFieldmap_m->getFieldDerivative(tmpA, tmpE_diff, tmpB_diff, zdir);
double bz = scale_m * tmpB(2);
double g = RefPartBunch_m->getGamma(i);
......
......@@ -172,7 +172,7 @@ void TravelingWave::addKR(int i, double t, Vector_t &K) {
wtf = frequency_m * t + phase_m;
CoreFieldmap_m->getFieldstrength(tmpA0, tmpE, tmpB);
CoreFieldmap_m->getFieldstrength_fdiff(tmpA0, tmpE_diff, tmpB_diff, zdir);
CoreFieldmap_m->getFieldDerivative(tmpA0, tmpE_diff, tmpB_diff, zdir);
k = scale_m * (tmpE_diff(2) * cos(wtf) - b * frequency_m * tmpE(2) * sin(wtf) / Physics::c);
} else if(tmpA0(2) < startExitField_m) {
......@@ -188,7 +188,7 @@ void TravelingWave::addKR(int i, double t, Vector_t &K) {
tmpB_diff = 0.0;
CoreFieldmap_m->getFieldstrength(tmpA0, tmpE, tmpB);
CoreFieldmap_m->getFieldstrength_fdiff(tmpA0, tmpE_diff, tmpB_diff, zdir);
CoreFieldmap_m->getFieldDerivative(tmpA0, tmpE_diff, tmpB_diff, zdir);
k = scaleCore_m * (tmpE_diff(2) * cos(wtf) - b * frequency_m * tmpE(2) * sin(wtf) / Physics::c);
wtf = frequency_m * t + phaseCore2_m;
......@@ -201,7 +201,7 @@ void TravelingWave::addKR(int i, double t, Vector_t &K) {
tmpB_diff = 0.0;
CoreFieldmap_m->getFieldstrength(tmpA0, tmpE, tmpB);
CoreFieldmap_m->getFieldstrength_fdiff(tmpA0, tmpE_diff, tmpB_diff, zdir);
CoreFieldmap_m->getFieldDerivative(tmpA0, tmpE_diff, tmpB_diff, zdir);
k += scaleCore_m * (tmpE_diff(2) * cos(wtf) - b * frequency_m * tmpE(2) * sin(wtf) / Physics::c);
} else {
......@@ -214,7 +214,7 @@ void TravelingWave::addKR(int i, double t, Vector_t &K) {
tmpB_diff = 0.0;
CoreFieldmap_m->getFieldstrength(tmpA0, tmpE, tmpB);
CoreFieldmap_m->getFieldstrength_fdiff(tmpA0, tmpE_diff, tmpB_diff, zdir);
CoreFieldmap_m->getFieldDerivative(tmpA0, tmpE_diff, tmpB_diff, zdir);
k = scale_m * (tmpE_diff(2) * cos(wtf) - b * frequency_m * tmpE(2) * sin(wtf) / Physics::c);
}
......@@ -523,18 +523,18 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
E.resize(N3, E0);
E2.resize(N3, E0);
for(int i = 1; i < N1; ++ i) {
E[i] = E0 + (F[i].first + F[i-1].first) / 2. * scale_m / mass;
E[i] = E0 + (F[i].first + F[i - 1].first) / 2. * scale_m / mass;
E2[i] = E[i];
}
for(int i = N1; i < N3 - N1 + 1; ++ i) {
int I = (i - N1) % N2 + N1;
double Z = (F[I].first + F[I-1].first) / 2. + floor((i - N1) / N2) * Dz;
double Z = (F[I].first + F[I - 1].first) / 2. + floor((i - N1) / N2) * Dz;
E[i] = E0 + Z * scaleCore_m / mass;
E2[i] = E[i];
}
for(int i = N3 - N1 + 1; i < N3; ++ i) {
int I = i - N3 - 1 + 2 * N1 + N2;
double Z = (F[I].first + F[I-1].first) / 2. + ((NumCells_m - 1) * Mode_m - 1) * Dz;
double Z = (F[I].first + F[I - 1].first) / 2. + ((NumCells_m - 1) * Mode_m - 1) * Dz;
E[i] = E0 + Z * scale_m / mass;
E2[i] = E[i];
}
......@@ -542,23 +542,23 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
for(int iter = 0; iter < 10; ++ iter) {
A = B = 0.0;
for(int i = 1; i < N1; ++ i) {
t[i] = t[i-1] + getdT(i, i, E, F, mass);
t2[i] = t2[i-1] + getdT(i, i, E2, F, mass);
t[i] = t[i - 1] + getdT(i, i, E, F, mass);
t2[i] = t2[i - 1] + getdT(i, i, E2, F, mass);
A += scale_m * (1. + frequency_m * (t2[i] - t[i]) / dphi) * getdA(i, i, t, 0.0, F);
B += scale_m * (1. + frequency_m * (t2[i] - t[i]) / dphi) * getdB(i, i, t, 0.0, F);
}
for(int i = N1; i < N3 - N1 + 1; ++ i) {
int I = (i - N1) % N2 + N1;
int J = (i - N1 + N4) % N2 + N1;
t[i] = t[i-1] + getdT(i, I, E, F, mass);
t2[i] = t2[i-1] + getdT(i, I, E2, F, mass);
t[i] = t[i - 1] + getdT(i, I, E, F, mass);
t2[i] = t2[i - 1] + getdT(i, I, E2, F, mass);
A += scaleCore_m * (1. + frequency_m * (t2[i] - t[i]) / dphi) * (getdA(i, I, t, phaseC1, F) + getdA(i, J, t, phaseC2, F));
B += scaleCore_m * (1. + frequency_m * (t2[i] - t[i]) / dphi) * (getdB(i, I, t, phaseC1, F) + getdB(i, J, t, phaseC2, F));
}
for(int i = N3 - N1 + 1; i < N3; ++ i) {
int I = i - N3 - 1 + 2 * N1 + N2;
t[i] = t[i-1] + getdT(i, I, E, F, mass);
t2[i] = t2[i-1] + getdT(i, I, E2, F, mass);
t[i] = t[i - 1] + getdT(i, I, E, F, mass);
t2[i] = t2[i - 1] + getdT(i, I, E2, F, mass);
A += scale_m * (1. + frequency_m * (t2[i] - t[i]) / dphi) * getdA(i, I, t, phaseE, F);
B += scale_m * (1. + frequency_m * (t2[i] - t[i]) / dphi) * getdB(i, I, t, phaseE, F);
}
......@@ -572,21 +572,21 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
tmp_phi += Physics::pi;
}
if(fabs(phi - tmp_phi) < frequency_m * (t[N3-1] - t[0]) / N3) {
if(fabs(phi - tmp_phi) < frequency_m * (t[N3 - 1] - t[0]) / N3) {
for(int i = 1; i < N1; ++ i) {
E[i] = E[i-1] + q * scale_m * getdE(i, i, t, phi, F);
E[i] = E[i - 1] + q * scale_m * getdE(i, i, t, phi, F);
}
for(int i = N1; i < N3 - N1 + 1; ++ i) {
int I = (i - N1) % N2 + N1;
int J = (i - N1 + N4) % N2 + N1;
E[i] = E[i-1] + q * scaleCore_m * (getdE(i, I, t, phi + phaseC1, F) + getdE(i, J, t, phi + phaseC2, F));
E[i] = E[i - 1] + q * scaleCore_m * (getdE(i, I, t, phi + phaseC1, F) + getdE(i, J, t, phi + phaseC2, F));
}
for(int i = N3 - N1 + 1; i < N3; ++ i) {
int I = i - N3 - 1 + 2 * N1 + N2;
E[i] = E[i-1] + q * scale_m * getdE(i, I, t, phi + phaseE, F);
E[i] = E[i - 1] + q * scale_m * getdE(i, I, t, phi + phaseE, F);
}
msg << "estimated phi= " << tmp_phi << " rad, "
<< "Ekin= " << E[N3-1] << " MeV" << endl;
<< "Ekin= " << E[N3 - 1] << " MeV" << endl;
// stringstream fn;
// fn << getName() << ".dat";
// ofstream ckrtest(fn.str().c_str());
......@@ -628,31 +628,31 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
for(int i = 1; i < N1; ++ i) {
E[i] = E[i-1] + q * scale_m * getdE(i, i, t, phi, F);
E2[i] = E2[i-1] + q * scale_m * getdE(i, i, t, phi + dphi, F); // should I use here t or t2?
t[i] = t[i-1] + getdT(i, i, E, F, mass);
t2[i] = t2[i-1] + getdT(i, i, E2, F, mass);
E[i] = E[i-1] + q * scale_m * getdE(i, i, t, phi, F);
E2[i] = E2[i-1] + q * scale_m * getdE(i, i, t2, phi + dphi, F);
E[i] = E[i - 1] + q * scale_m * getdE(i, i, t, phi, F);
E2[i] = E2[i - 1] + q * scale_m * getdE(i, i, t, phi + dphi, F); // should I use here t or t2?
t[i] = t[i - 1] + getdT(i, i, E, F, mass);
t2[i] = t2[i - 1] + getdT(i, i, E2, F, mass);
E[i] = E[i - 1] + q * scale_m * getdE(i, i, t, phi, F);
E2[i] = E2[i - 1] + q * scale_m * getdE(i, i, t2, phi + dphi, F);
}
for(int i = N1; i < N3 - N1 + 1; ++ i) {
int I = (i - N1) % N2 + N1;
int J = (i - N1 + N4) % N2 + N1;
E[i] = E[i-1] + q * scaleCore_m * (getdE(i, I, t, phi + phaseC1, F) + getdE(i, J, t, phi + phaseC2, F));
E2[i] = E2[i-1] + q * scaleCore_m * (getdE(i, I, t, phi + phaseC1 + dphi, F) + getdE(i, J, t, phi + phaseC2 + dphi, F)); //concerning t: see above
t[i] = t[i-1] + getdT(i, I, E, F, mass);
t2[i] = t2[i-1] + getdT(i, I, E2, F, mass);
E[i] = E[i-1] + q * scaleCore_m * (getdE(i, I, t, phi + phaseC1, F) + getdE(i, J, t, phi + phaseC2, F));
E2[i] = E2[i-1] + q * scaleCore_m * (getdE(i, I, t2, phi + phaseC1 + dphi, F) + getdE(i, J, t2, phi + phaseC2 + dphi, F));
E[i] = E[i - 1] + q * scaleCore_m * (getdE(i, I, t, phi + phaseC1, F) + getdE(i, J, t, phi + phaseC2, F));
E2[i] = E2[i - 1] + q * scaleCore_m * (getdE(i, I, t, phi + phaseC1 + dphi, F) + getdE(i, J, t, phi + phaseC2 + dphi, F)); //concerning t: see above
t[i] = t[i - 1] + getdT(i, I, E, F, mass);
t2[i] = t2[i - 1] + getdT(i, I, E2, F, mass);
E[i] = E[i - 1] + q * scaleCore_m * (getdE(i, I, t, phi + phaseC1, F) + getdE(i, J, t, phi + phaseC2, F));
E2[i] = E2[i - 1] + q * scaleCore_m * (getdE(i, I, t2, phi + phaseC1 + dphi, F) + getdE(i, J, t2, phi + phaseC2 + dphi, F));
}
for(int i = N3 - N1 + 1; i < N3; ++ i) {
int I = i - N3 - 1 + 2 * N1 + N2;
E[i] = E[i-1] + q * scale_m * getdE(i, I, t, phi + phaseE, F);
E2[i] = E2[i-1] + q * scale_m * getdE(i, I, t, phi + phaseE + dphi, F); //concerning t: see above
t[i] = t[i-1] + getdT(i, I, E, F, mass);
t2[i] = t2[i-1] + getdT(i, I, E2, F, mass);
E[i] = E[i-1] + q * scale_m * getdE(i, I, t, phi + phaseE, F);
E2[i] = E2[i-1] + q * scale_m * getdE(i, I, t2, phi + phaseE + dphi, F);
E[i] = E[i - 1] + q * scale_m * getdE(i, I, t, phi + phaseE, F);
E2[i] = E2[i - 1] + q * scale_m * getdE(i, I, t, phi + phaseE + dphi, F); //concerning t: see above
t[i] = t[i - 1] + getdT(i, I, E, F, mass);
t2[i] = t2[i - 1] + getdT(i, I, E2, F, mass);
E[i] = E[i - 1] + q * scale_m * getdE(i, I, t, phi + phaseE, F);
E2[i] = E2[i - 1] + q * scale_m * getdE(i, I, t2, phi + phaseE + dphi, F);
}
// msg << ", Ekin= " << E[N3-1] << " MeV" << endl;
}
......
......@@ -697,13 +697,13 @@ void PartBunch::computeSelfFields(int binNumber) {
if(fs_m->hasValidSolver()) {
/// Scatter charge onto space charge grid.
this->Q *= this->dt;
if (!interpolationCacheSet_m) {
if (interpolationCache_m.size() < getLocalNum()) {
if(!interpolationCacheSet_m) {
if(interpolationCache_m.size() < getLocalNum()) {
interpolationCache_m.create(getLocalNum() - interpolationCache_m.size());
} else {
interpolationCache_m.destroy(interpolationCache_m.size() - getLocalNum(),
getLocalNum(),
true);
getLocalNum(),
true);
}
interpolationCacheSet_m = true;
......@@ -2432,20 +2432,20 @@ void PartBunch::boundp_destroy() {
const int checkfactor = Options::remotePartDel;
// check the bunch if its full size is larger than checkfactor times of its rms size
if (checkfactor != -1) {
if(len[0] > checkfactor * rrms_m[0] || len[1] > checkfactor * rrms_m[1] || len[2] > checkfactor * rrms_m[2]) {
for(unsigned int ii = 0; ii < this->getLocalNum(); ii++) {
// delete the particle if the ditance to the beam center is larger than 8 times of beam's rms size
if(abs(R[ii](0) - rmean_m(0)) > checkfactor * rrms_m[0] || abs(R[ii](1) - rmean_m(1)) > checkfactor * rrms_m[1] || abs(R[ii](2) - rmean_m(2)) > checkfactor * rrms_m[2]) {
// put particle onto deletion list
destroy(1, ii);
//update bin parameter
if(weHaveBins()) countLost[Bin[ii]] += 1 ;
gmsgAll << "REMOTE PARTICLE DELETION: ID = " << ID[ii] << ", R = " << R[ii] << ", beam rms = " << rrms_m << endl;
if(checkfactor != -1) {
if(len[0] > checkfactor * rrms_m[0] || len[1] > checkfactor * rrms_m[1] || len[2] > checkfactor * rrms_m[2]) {
for(unsigned int ii = 0; ii < this->getLocalNum(); ii++) {
// delete the particle if the ditance to the beam center is larger than 8 times of beam's rms size
if(abs(R[ii](0) - rmean_m(0)) > checkfactor * rrms_m[0] || abs(R[ii](1) - rmean_m(1)) > checkfactor * rrms_m[1] || abs(R[ii](2) - rmean_m(2)) > checkfactor * rrms_m[2]) {
// put particle onto deletion list
destroy(1, ii);
//update bin parameter
if(weHaveBins()) countLost[Bin[ii]] += 1 ;
gmsgAll << "REMOTE PARTICLE DELETION: ID = " << ID[ii] << ", R = " << R[ii] << ", beam rms = " << rrms_m << endl;
}
}
}
}
}
for(int i = 0; i < dimIdx; i++) {
rmax_m[i] += dh_m * abs(rmax_m[i] - rmin_m[i]);
......@@ -2714,7 +2714,7 @@ void PartBunch::pop() {
if(!bunchStashed_m) return;
size_t Nloc = getLocalNum();
if (getTotalNum() > 0) {
if(getTotalNum() > 0) {
destroy(Nloc, 0);
}
update();
......@@ -2748,3 +2748,34 @@ void PartBunch::pop() {
update();
}
double PartBunch::getZPos() {
if(sum(PType != 0)) {
size_t numberOfPrimaryParticles = 0;
double zAverage = 0.0;
if(getLocalNum() != 0) {
for(size_t partIndex = 0; partIndex < getLocalNum(); partIndex++) {
if(PType[partIndex] == 0) {
zAverage += X[partIndex](2);
numberOfPrimaryParticles++;
}
}
if(numberOfPrimaryParticles != 0)
zAverage /= numberOfPrimaryParticles;
}
reduce(zAverage, zAverage, OpAddAssign());
zAverage /= Ippl::getNodes();
return zAverage;
} else {
if(getTotalNum() > 0)
return sum(X(2)) / getTotalNum();
else
return 0.0;
}
}
void PartBunch::getXBounds(Vector_t &xMin, Vector_t &xMax) {
bounds(X, xMin, xMax);
}
......@@ -107,8 +107,8 @@ public:
bool addDistributions(std::vector<Distribution *> distributions, size_t numberOfParticles);
void setGridIsFixed();
bool isGridFixed();
void setGridIsFixed();
bool isGridFixed();
/*
......@@ -161,7 +161,7 @@ public:
void calcGammas_cycl();
/** \brief Get gamma of one bin */
double getBinGamma(int bin);
double getBinGamma(int bin);
/** \brief Set the charge of one bin to the value of q and all other to zero */
void setBinCharge(int bin, double q);
......@@ -173,7 +173,7 @@ public:
double getMaxdEBins();
/** \brief calculates back the max/min of the efield on the grid */
std::pair<Vector_t, Vector_t> getEExtrema();
std::pair<Vector_t, Vector_t> getEExtrema();
/*
......@@ -181,11 +181,11 @@ public:
*/
const Mesh_t &getMesh() const;
const Mesh_t &getMesh() const;
Mesh_t &getMesh();
Mesh_t &getMesh();
FieldLayout_t &getFieldLayout();
FieldLayout_t &getFieldLayout();
void setBCAllOpen();
......@@ -294,6 +294,12 @@ public:
*/
double get_sPos();
/// Get average z position from local "lab frame" coordinates, X.
double getZPos();
/// Get bounds of local "lab frame" coordinates, X.
void getXBounds(Vector_t &xMin, Vector_t &xMax);
double get_phase() const;
double get_gamma() const;
......@@ -374,7 +380,7 @@ public:
/// step in a TRACK command
inline void setLocalTrackStep(long long n) {localTrackStep_m = n;}
inline void incTrackSteps() {globalTrackStep_m++;localTrackStep_m++;}
inline void incTrackSteps() {globalTrackStep_m++; localTrackStep_m++;}
inline long long getLocalTrackStep() const {return localTrackStep_m;}
inline void setNumBunch(int n);
......@@ -566,7 +572,7 @@ private:
bool interpolationCacheSet_m;
ParticleAttrib<CacheDataCIC<double,3U> > interpolationCache_m;
ParticleAttrib<CacheDataCIC<double, 3U> > interpolationCache_m;
/// counter to store the distributin dump
int distDump_m;
......@@ -588,7 +594,7 @@ private:
ParticleAttrib<short> stash_ptype_m;
bool bunchStashed_m;
PartBunch& operator=(const PartBunch&) = delete;
PartBunch &operator=(const PartBunch &) = delete;
///
int fieldDBGStep_m;
......@@ -830,17 +836,17 @@ std::pair<Vector_t, Vector_t> PartBunch::getEExtrema() {
}
inline
const Mesh_t & PartBunch::getMesh() const {
const Mesh_t &PartBunch::getMesh() const {
return getLayout().getLayout().getMesh();
}
inline
Mesh_t & PartBunch::getMesh() {
Mesh_t &PartBunch::getMesh() {
return getLayout().getLayout().getMesh();
}
inline
FieldLayout_t & PartBunch::getFieldLayout() {
FieldLayout_t &PartBunch::getFieldLayout() {
return dynamic_cast<FieldLayout_t &>(getLayout().getLayout().getFieldLayout());
}
......@@ -948,7 +954,7 @@ double PartBunch::getT() const {
inline
void PartBunch::resetInterpolationCache(bool clearCache) {
interpolationCacheSet_m = false;
if (clearCache) {
if(clearCache) {
interpolationCache_m.destroy(interpolationCache_m.size(), 0, true);
}
}
......
......@@ -86,7 +86,7 @@ void Astra1DDynamic::readMap() {
double Ez_max = 0.0;
double dz = (zend_m - zbegin_m) / (num_gridpz_m - 1);
double *RealValues = new double[2*num_gridpz_m];
double *RealValues = new double[2 * num_gridpz_m];
double *zvals = new double[num_gridpz_m];
gsl_spline *Ez_interpolant = gsl_spline_alloc(gsl_interp_cspline, num_gridpz_m);
......@@ -95,7 +95,7 @@ void Astra1DDynamic::readMap() {
gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(2 * num_gridpz_m);
gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(2 * num_gridpz_m);
FourCoefs_m = new double[2*accuracy_m - 1];
FourCoefs_m = new double[2 * accuracy_m - 1];
// read in and parse field map
in.open(Filename_m.c_str());
......@@ -150,7 +150,7 @@ void Astra1DDynamic::readMap() {
delete[] zvals;
delete[] RealValues;
INFOMSG( typeset_msg("read in fieldmap '" + Filename_m + "'", "info") << endl);
INFOMSG(typeset_msg("read in fieldmap '" + Filename_m + "'", "info") << endl);
}
}
......@@ -184,13 +184,13 @@ bool Astra1DDynamic::getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &
somefactor = 1.0;
coskzl = cos(kz * l);
sinkzl = sin(kz * l);
ez += (FourCoefs_m[n] * coskzl - FourCoefs_m[n+1] * sinkzl);
ez += (FourCoefs_m[n] * coskzl - FourCoefs_m[n + 1] * sinkzl);
somefactor *= somefactor_base;
ezp += somefactor * (-FourCoefs_m[n] * sinkzl - FourCoefs_m[n+1] * coskzl);
ezp += somefactor * (-FourCoefs_m[n] * sinkzl - FourCoefs_m[n + 1] * coskzl);
somefactor *= somefactor_base;
ezpp += somefactor * (-FourCoefs_m[n] * coskzl + FourCoefs_m[n+1] * sinkzl);
ezpp += somefactor * (-FourCoefs_m[n] * coskzl + FourCoefs_m[n + 1] * sinkzl);
somefactor *= somefactor_base;
ezppp += somefactor * (FourCoefs_m[n] * sinkzl + FourCoefs_m[n+1] * coskzl);
ezppp += somefactor * (FourCoefs_m[n] * sinkzl + FourCoefs_m[n + 1] * coskzl);
}
// expand the field off-axis
const double f = -(ezpp + ez * xlrep_m * xlrep_m) / 16.;
......@@ -208,13 +208,13 @@ bool Astra1DDynamic::getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &
return false;
}
bool Astra1DDynamic::getFieldstrength_fdiff(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const {
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;
double ezp = 0.0;
int n = 1;
for(int l = 1; l < accuracy_m; l++, n += 2)
ezp += two_pi / length_m * l * (-FourCoefs_m[n] * sin(kz * l) - FourCoefs_m[n+1] * cos(kz * l));
ezp += two_pi / length_m * l * (-FourCoefs_m[n] * sin(kz * l) - FourCoefs_m[n + 1] * cos(kz * l));
E(2) += ezp;
......@@ -253,7 +253,7 @@ void Astra1DDynamic::getOnaxisEz(vector<pair<double, double> > & F) {
ifstream in(Filename_m.c_str());
interpreteLine<string, int>(in, tmpString, tmpInt);
interpreteLine<double>(in, tmpDouble);
for(int i = 0; i < num_gridpz_m; ++ i) {
interpreteLine<double, double>(in, F[i].first, F[i].second);
if(fabs(F[i].second) > Ez_max) {
......
......@@ -7,7 +7,7 @@ class Astra1DDynamic: public Fieldmap {
public:
virtual bool getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const;
virtual bool getFieldstrength_fdiff(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const;
virtual bool getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const;
virtual void getFieldDimensions(double &zBegin, double &zEnd, double &rBegin, double &rEnd) const;
virtual void getFieldDimensions(double &xIni, double &xFinal, double &yIni, double &yFinal, double &zIni, double &zFinal) const;
virtual void swap();
......
......@@ -12,7 +12,7 @@ using Physics::c;
using Physics::two_pi;
Astra1DDynamic_fast::Astra1DDynamic_fast(string aFilename):
Fieldmap(aFilename){
Fieldmap(aFilename) {
Inform msg("*1DD ");
ifstream file;
int tmpInt;
......@@ -97,7 +97,7 @@ void Astra1DDynamic_fast::readMap() {
double coskzl, sinkzl;
double *higherDerivatives[3];
double *zvals;
double *RealValues = new double[2*num_gridpz_m];
double *RealValues = new double[2 * num_gridpz_m];
onAxisField_m = new double[num_gridpz_m];
zvals_m = new double[num_gridpz_m];
......@@ -157,7 +157,7 @@ void Astra1DDynamic_fast::readMap() {
RealValues[ii ++] = onAxisField_m[num_gridpz_m - 1];
// prepend mirror sampling points such that field values are periodic for sure
// ii == 2*num_gridpz_m at the moment
-- ii;
-- ii;
for(int i = 0; i < num_gridpz_m; ++ i, -- ii) {
RealValues[i] = RealValues[ii];
}
......@@ -165,7 +165,7 @@ void Astra1DDynamic_fast::readMap() {
gsl_fft_real_transform(RealValues, 1, 2 * num_gridpz_m, real, work);
// normalize to Ez_max = 1 MV/m
RealValues[0] /= (2 * num_gridpz_m);
RealValues[0] /= (2 * num_gridpz_m);
for(int i = 1; i < 2 * accuracy - 1; i++) {
RealValues[i] /= num_gridpz_m;
}
......@@ -181,12 +181,12 @@ void Astra1DDynamic_fast::readMap() {
interior_derivative = base;
coskzl = cos(kz * l);
sinkzl = sin(kz * l);
higherDerivatives[0][i] += interior_derivative * (-RealValues[n] * sinkzl - RealValues[n+1] * coskzl);
higherDerivatives[0][i] += interior_derivative * (-RealValues[n] * sinkzl - RealValues[n + 1] * coskzl);
interior_derivative *= base;
higherDerivatives[1][i] += interior_derivative * (-RealValues[n] * coskzl + RealValues[n+1] * sinkzl);
higherDerivatives[1][i] += interior_derivative * (-RealValues[n] * coskzl + RealValues[n + 1] * sinkzl);
interior_derivative *= base;
higherDerivatives[2][i] += interior_derivative * (RealValues[n] * sinkzl + RealValues[n+1] * coskzl);
higherDerivatives[2][i] += interior_derivative * (RealValues[n] * sinkzl + RealValues[n + 1] * coskzl);
}
}
......@@ -194,7 +194,7 @@ void Astra1DDynamic_fast::readMap() {
gsl_fft_real_wavetable_free(real);
for(int i = 1; i < 4; ++i) {
gsl_spline_init(onAxisInterpolants_m[i], zvals, higherDerivatives[i-1], num_gridpz_m);
gsl_spline_init(onAxisInterpolants_m[i], zvals, higherDerivatives[i - 1], num_gridpz_m);
}
delete[] RealValues;
......@@ -203,7 +203,7 @@ void Astra1DDynamic_fast::readMap() {
delete[] higherDerivatives[2];