Commit cfba4cff authored by ext-calvo_p's avatar ext-calvo_p

Part of "Static code analysis of OPAL"

parent be8cf664
......@@ -21,8 +21,9 @@
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "AbsBeamline/BeamlineVisitor.h"
#include "AbsBeamline/BeamStripping.h"
#include "AbsBeamline/BeamlineVisitor.h"
#include "AbsBeamline/Cyclotron.h"
#include "Algorithms/PartBunchBase.h"
#include "Fields/Fieldmap.h"
......@@ -30,22 +31,19 @@
#include "Solvers/BeamStrippingPhysics.hh"
#include "Solvers/ParticleMatterInteractionHandler.hh"
#include "Utilities/LogicalError.h"
#include "Utilities/Options.h"
#include "Utilities/Util.h"
#include "Utilities/GeneralClassicException.h"
#include <memory>
#include <fstream>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <cstdio>
#include "Ippl.h"
#include "gsl/gsl_spline.h"
#include "gsl/gsl_interp.h"
#define CHECK_BSTP_FSCANF_EOF(arg) if(arg == EOF)\
#define CHECK_BSTP_FSCANF_EOF(arg) if (arg == EOF)\
throw GeneralClassicException("BeamStripping::getPressureFromFile",\
"fscanf returned EOF at " #arg);
......@@ -104,7 +102,7 @@ void BeamStripping::setPressure(double pressure) {
}
double BeamStripping::getPressure() const {
if(pressure_m > 0.0)
if (pressure_m > 0.0)
return pressure_m;
else {
throw LogicalError("BeamStripping::getPressure",
......@@ -117,7 +115,7 @@ void BeamStripping::setTemperature(double temperature) {
}
double BeamStripping::getTemperature() const {
if(temperature_m > 0.0)
if (temperature_m > 0.0)
return temperature_m;
else {
throw LogicalError("BeamStripping::getTemperature",
......@@ -138,7 +136,7 @@ void BeamStripping::setPScale(double ps) {
}
double BeamStripping::getPScale() const {
if(pscale_m > 0.0)
if (pscale_m > 0.0)
return pscale_m;
else {
throw LogicalError("BeamStripping::getPScale",
......@@ -151,7 +149,7 @@ void BeamStripping::setResidualGas(std::string gas) {
}
std::string BeamStripping::getResidualGas() const {
if(gas_m == "H2" || gas_m == "AIR")
if (gas_m == "H2" || gas_m == "AIR")
return gas_m;
else {
throw GeneralClassicException("BeamStripping::getResidualGas",
......@@ -184,10 +182,9 @@ bool BeamStripping::checkBeamStripping(PartBunchBase<double, 3> *bunch, Cyclotro
maxz_m = 1000 * cycl->getMaxZ();
minz_m = 1000 * cycl->getMinZ();
int pflag = 0;
size_t tempnum = bunch->getLocalNum();
for (unsigned int i = 0; i < tempnum; ++i) {
pflag = checkPoint(bunch->R[i](0), bunch->R[i](1), bunch->R[i](2));
int pflag = checkPoint(bunch->R[i](0), bunch->R[i](1), bunch->R[i](2));
if ( (pflag != 0) && (bunch->Bin[i] != -1) )
flagNeedUpdate = true;
}
......@@ -224,7 +221,7 @@ void BeamStripping::initialise(PartBunchBase<double, 3> *bunch, const double &sc
for (size_t i = 0; i < bunch->getLocalNum(); ++i) {
bunch->M[i] = bunch->getM()*1E-9;
bunch->Q[i] = bunch->getQ() * Physics::q_e;
if(bunch->weHaveBins())
if (bunch->weHaveBins())
bunch->Bin[bunch->getLocalNum()-1] = bunch->Bin[i];
}
}
......@@ -271,7 +268,7 @@ double BeamStripping::checkPressure(const double &x, const double &y) {
double pressure = 0.0;
if(pmapfn_m != "") {
if (pmapfn_m != "") {
const double rad = std::sqrt(x * x + y * y);
const double xir = (rad - PP.rmin) / PP.delr;
......@@ -282,13 +279,13 @@ double BeamStripping::checkPressure(const double &x, const double &y) {
// wr2 : the relative distance to the outer path radius
const double wr2 = 1.0 - wr1;
const double tempv = atan(y / x);
const double tempv = std::atan(y / x);
double tet = tempv;
if((x < 0) && (y >= 0)) tet = Physics::pi + tempv;
else if((x < 0) && (y <= 0)) tet = Physics::pi + tempv;
else if((x > 0) && (y <= 0)) tet = Physics::two_pi + tempv;
else if((x == 0) && (y > 0)) tet = Physics::pi / 2.0;
else if((x == 0) && (y < 0)) tet = 1.5 * Physics::pi;
if ((x < 0) && (y >= 0)) tet = Physics::pi + tempv;
else if ((x < 0) && (y <= 0)) tet = Physics::pi + tempv;
else if ((x > 0) && (y <= 0)) tet = Physics::two_pi + tempv;
else if ((x == 0) && (y > 0)) tet = Physics::pi / 2.0;
else if ((x == 0) && (y < 0)) tet = 1.5 * Physics::pi;
// the actual angle of particle
tet = tet * Physics::rad2deg;
......@@ -314,12 +311,12 @@ double BeamStripping::checkPressure(const double &x, const double &y) {
r1t2 = idx(ir, it + 1);
r2t2 = idx(ir + 1, it + 1);
if((it >= 0) && (ir >= 0) && (it < PField.ntetS) && (ir < PField.nrad)) {
if ((it >= 0) && (ir >= 0) && (it < PField.ntetS) && (ir < PField.nrad)) {
pressure = (PField.pfld[r1t1] * wr2 * wt2 +
PField.pfld[r2t1] * wr1 * wt2 +
PField.pfld[r1t2] * wr2 * wt1 +
PField.pfld[r2t2] * wr1 * wt1);
if(pressure <= 0.0) {
if (pressure <= 0.0) {
*gmsg << level4 << getName() << ": Pressure data from file zero." << endl;
*gmsg << level4 << getName() << ": Take constant value through BeamStripping::getPressure" << endl;
pressure = getPressure();
......@@ -345,7 +342,7 @@ double BeamStripping::checkPressure(const double &x, const double &y) {
// Calculates radius of initial grid (dimensions in [m]!)
void BeamStripping::initR(double rmin, double dr, int nrad) {
PP.rarr.resize(nrad);
for(int i = 0; i < nrad; i++)
for (int i = 0; i < nrad; i++)
PP.rarr[i] = rmin + i * dr;
PP.delr = dr;
......@@ -361,31 +358,31 @@ void BeamStripping::getPressureFromFile(const double &scaleFactor) {
PP.Pfact = scaleFactor;
if((f = fopen(pmapfn_m.c_str(), "r")) == NULL) {
if ((f = std::fopen(pmapfn_m.c_str(), "r")) == NULL) {
throw GeneralClassicException("BeamStripping::getPressureFromFile",
"failed to open file '" + pmapfn_m + "', please check if it exists");
}
CHECK_BSTP_FSCANF_EOF(fscanf(f, "%lf", &PP.rmin));
CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%lf", &PP.rmin));
*gmsg << "* --- Minimal radius of measured pressure map: " << PP.rmin << " [mm]" << endl;
CHECK_BSTP_FSCANF_EOF(fscanf(f, "%lf", &PP.delr));
CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%lf", &PP.delr));
//if the value is negative, the actual value is its reciprocal.
if(PP.delr < 0.0) PP.delr = 1.0 / (-PP.delr);
if (PP.delr < 0.0) PP.delr = 1.0 / (-PP.delr);
*gmsg << "* --- Stepsize in radial direction: " << PP.delr << " [mm]" << endl;
CHECK_BSTP_FSCANF_EOF(fscanf(f, "%lf", &PP.tetmin));
CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%lf", &PP.tetmin));
*gmsg << "* --- Minimal angle of measured pressure map: " << PP.tetmin << " [deg]" << endl;
CHECK_BSTP_FSCANF_EOF(fscanf(f, "%lf", &PP.dtet));
CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%lf", &PP.dtet));
//if the value is negative, the actual value is its reciprocal.
if(PP.dtet < 0.0) PP.dtet = 1.0 / (-PP.dtet);
if (PP.dtet < 0.0) PP.dtet = 1.0 / (-PP.dtet);
*gmsg << "* --- Stepsize in azimuthal direction: " << PP.dtet << " [deg]" << endl;
CHECK_BSTP_FSCANF_EOF(fscanf(f, "%d", &PField.ntet));
CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%d", &PField.ntet));
*gmsg << "* --- Grid points along azimuth (ntet): " << PField.ntet << endl;
CHECK_BSTP_FSCANF_EOF(fscanf(f, "%d", &PField.nrad));
CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%d", &PField.nrad));
*gmsg << "* --- Grid points along radius (nrad): " << PField.nrad << endl;
*gmsg << "* --- Maximum radial position: " << PP.rmin + (PField.nrad-1)*PP.delr << " [mm]" << endl;
PP.rmin *= 0.001; // mm --> m
......@@ -399,15 +396,15 @@ void BeamStripping::getPressureFromFile(const double &scaleFactor) {
*gmsg << "* --- Total stored grid point number ((ntet+1) * nrad) : " << PField.ntot << endl;
*gmsg << "* --- Escale factor: " << PP.Pfact << endl;
for(int i = 0; i < PField.nrad; i++) {
for(int k = 0; k < PField.ntet; k++) {
CHECK_BSTP_FSCANF_EOF(fscanf(f, "%16lE", &(PField.pfld[idx(i, k)])));
for (int i = 0; i < PField.nrad; i++) {
for (int k = 0; k < PField.ntet; k++) {
CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%16lE", &(PField.pfld[idx(i, k)])));
PField.pfld[idx(i, k)] *= PP.Pfact;
}
}
*gmsg << "*" << endl;
fclose(f);
std::fclose(f);
}
#undef CHECK_BSTP_FSCANF_EOF
\ No newline at end of file
......@@ -135,7 +135,7 @@ private:
std::string pmapfn_m; /// stores the filename of the pressure map
double pscale_m; /// a scale factor for the P-field
double temperature_m; /// K
double stop_m;
bool stop_m; /// Flag if particles should be stripped or stopped
///@}
///@{ size limits took from cyclotron
......
This diff is collapsed.
......@@ -67,7 +67,7 @@ Degrader::Degrader(const std::string &name):
Degrader::~Degrader() {
if(online_m)
if (online_m)
goOffline();
}
......@@ -91,7 +91,7 @@ bool Degrader::apply(const size_t &i, const double &t, Vector_t &/*E*/, Vector_t
const Vector_t &P = RefPartBunch_m->P[i];
const double recpgamma = Physics::c * RefPartBunch_m->dt[i] / std::sqrt(1.0 + dot(P, P));
if(isInMaterial(R(2))) {
if (isInMaterial(R(2))) {
//if particle was allready marked as -1 (it means it should have gone into degrader but didn't)
//set the label to -2 (will not go into degrader and will be deleted when particles per core > 2)
if (RefPartBunch_m->Bin[i] < 0)
......
......@@ -369,7 +369,6 @@ double RFCavity::getCycFrequency()const {
void RFCavity::getMomentaKick(const double normalRadius, double momentum[], const double t, const double dtCorrt, const int PID, const double restMass, const int chargenumber) {
double derivate;
double Voltage;
double momentum2 = momentum[0] * momentum[0] + momentum[1] * momentum[1] + momentum[2] * momentum[2];
double betgam = std::sqrt(momentum2);
......@@ -377,15 +376,14 @@ void RFCavity::getMomentaKick(const double normalRadius, double momentum[], cons
double gamma = std::sqrt(1.0 + momentum2);
double beta = betgam / gamma;
Voltage = spline(normalRadius, &derivate) * scale_m * 1.0e6; // V
double Voltage = spline(normalRadius, &derivate) * scale_m * 1.0e6; // V
double transit_factor = 0.0;
double Ufactor = 1.0;
double frequency = frequency_m * frequency_td_m->getValue(t);
if (gapwidth_m > 0.0) {
transit_factor = 0.5 * frequency * gapwidth_m * 1.0e-3 / (Physics::c * beta);
double transit_factor = 0.5 * frequency * gapwidth_m * 1.0e-3 / (Physics::c * beta);
Ufactor = std::sin(transit_factor) / transit_factor;
}
......@@ -543,8 +541,6 @@ double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const
gsl_spline *onAxisInterpolants;
gsl_interp_accel *onAxisAccel;
unsigned int N;
double A, B;
double phi = 0.0, tmp_phi, dphi = 0.5 * Physics::pi / 180.;
double dz = 1.0, length = 0.0;
fieldmap_m->getOnaxisEz(G);
......@@ -566,7 +562,7 @@ double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const
G.clear();
N = (int)floor(length / dz + 1);
unsigned int N = (int)floor(length / dz + 1);
dz = length / N;
F.resize(N);
......@@ -589,7 +585,8 @@ double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const
}
for (int iter = 0; iter < 10; ++ iter) {
A = B = 0.0;
double A = 0.0;
double B = 0.0;
for (unsigned int i = 1; i < N; ++ i) {
t[i] = t[i - 1] + getdT(i, E, dz, mass);
t2[i] = t2[i - 1] + getdT(i, E2, dz, mass);
......
......@@ -101,7 +101,7 @@ bool Stripper::doPreCheck(PartBunchBase<double, 3> *bunch) {
double ymax = std::max(std::abs(rmin(1)), std::abs(rmax(1)));
double rbunch_max = std::hypot(xmax, ymax);
if(rbunch_max > rmin_m - 10.0) {
if (rbunch_max > rmin_m - 10.0) {
return true;
}
return false;
......@@ -116,13 +116,13 @@ bool Stripper::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, co
size_t tempnum = bunch->getLocalNum();
Inform gmsgALL("OPAL", INFORM_ALL_NODES);
for(unsigned int i = 0; i < tempnum; ++i) {
if(bunch->PType[i] != ParticleType::REGULAR) continue;
for (unsigned int i = 0; i < tempnum; ++i) {
if (bunch->PType[i] != ParticleType::REGULAR) continue;
double tangle = calculateIncidentAngle(bunch->P[i](0), bunch->P[i](1));
changeWidth(bunch, i, tstep, tangle);
int pflag = checkPoint(bunch->R[i](0), bunch->R[i](1));
if(pflag == 0) continue;
if (pflag == 0) continue;
// dist1 > 0, right hand, dt > 0; dist1 < 0, left hand, dt < 0
double dist1 = (A_m * bunch->R[i](0) + B_m * bunch->R[i](1) + C_m) / R_m * 1.0e-3; // [m]
......@@ -142,7 +142,7 @@ bool Stripper::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, co
} else {
gmsgALL << level4 << getName() << ": Particle " << bunch->ID[i] << " collide in stripper " << getName() << endl;
// change charge and mass of PartData when the reference particle hits the stripper.
if(bunch->ID[i] == 0)
if (bunch->ID[i] == 0)
bunch->setPType(ParticleType::STRIPPED);
// change the mass and charge
......@@ -161,7 +161,7 @@ bool Stripper::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, co
bunch->M[index] = bunch->M[i];
// once the particle is stripped, change PType from 0 to 1 as a flag so as to avoid repetitive stripping.
bunch->PType[index] = ParticleType::STRIPPED;
if(bunch->weHaveBins())
if (bunch->weHaveBins())
bunch->Bin[index] = bunch->Bin[i];
count++;
......@@ -175,7 +175,7 @@ bool Stripper::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, co
bool Stripper::doFinaliseCheck(PartBunchBase<double, 3> *bunch, bool flagNeedUpdate) {
reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign());
if(!stop_m){
if (!stop_m){
// change charge and mass of PartData when the reference particle hits the stripper.
if (bunch->getPType() == ParticleType::STRIPPED) {
bunch->resetM(opmass_m * 1.0e9); // GeV -> eV
......
......@@ -268,53 +268,51 @@ ElementBase::ElementType TravelingWave::getType() const {
double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, const double &q, const double &mass) {
std::vector<double> t, E, t2, E2;
std::vector<std::pair<double, double> > F;
double Dz;
int N1, N2, N3, N4;
double A, B;
double phi = 0.0, tmp_phi, dphi = 0.5 * Physics::pi / 180.;
double phaseC1 = phaseCore1_m - phase_m;
double phaseC2 = phaseCore2_m - phase_m;
double phaseE = phaseExit_m - phase_m;
fieldmap_m->getOnaxisEz(F);
if(F.size() == 0) return 0.0;
if (F.size() == 0) return 0.0;
N1 = static_cast<int>(std::floor(F.size() / 4.)) + 1;
N2 = F.size() - 2 * N1 + 1;
N3 = 2 * N1 + static_cast<int>(std::floor((NumCells_m - 1) * N2 * Mode_m)) - 1;
N4 = static_cast<int>(std::round(N2 * Mode_m));
Dz = F[N1 + N2].first - F[N1].first;
int N1 = static_cast<int>(std::floor(F.size() / 4.)) + 1;
int N2 = F.size() - 2 * N1 + 1;
int N3 = 2 * N1 + static_cast<int>(std::floor((NumCells_m - 1) * N2 * Mode_m)) - 1;
int N4 = static_cast<int>(std::round(N2 * Mode_m));
double Dz = F[N1 + N2].first - F[N1].first;
t.resize(N3, t0);
t2.resize(N3, t0);
E.resize(N3, E0);
E2.resize(N3, E0);
for(int i = 1; i < N1; ++ i) {
for (int i = 1; i < N1; ++ i) {
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) {
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. + std::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) {
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;
E[i] = E0 + Z * scale_m / mass;
E2[i] = E[i];
}
for(int iter = 0; iter < 10; ++ iter) {
A = B = 0.0;
for(int i = 1; i < N1; ++ i) {
for (int iter = 0; iter < 10; ++ iter) {
double A = 0.0;
double 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);
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) {
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);
......@@ -322,7 +320,7 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
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) {
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);
......@@ -330,25 +328,25 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
B += scale_m * (1. + frequency_m * (t2[i] - t[i]) / dphi) * getdB(i, I, t, phaseE, F);
}
if(std::abs(B) > 0.0000001) {
if (std::abs(B) > 0.0000001) {
tmp_phi = std::atan(A / B);
} else {
tmp_phi = Physics::pi / 2;
}
if(q * (A * std::sin(tmp_phi) + B * std::cos(tmp_phi)) < 0) {
if (q * (A * std::sin(tmp_phi) + B * std::cos(tmp_phi)) < 0) {
tmp_phi += Physics::pi;
}
if(std::abs(phi - tmp_phi) < frequency_m * (t[N3 - 1] - t[0]) / N3) {
for(int i = 1; i < N1; ++ i) {
if (std::abs(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);
}
for(int i = N1; i < N3 - N1 + 1; ++ i) {
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));
}
for(int i = N3 - N1 + 1; i < N3; ++ i) {
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);
}
......@@ -362,7 +360,7 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
phi = tmp_phi - std::round(tmp_phi / Physics::two_pi) * Physics::two_pi;
for(int i = 1; i < N1; ++ i) {
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);
......@@ -370,7 +368,7 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
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) {
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));
......@@ -380,7 +378,7 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
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) {
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
......@@ -405,4 +403,4 @@ bool TravelingWave::isInside(const Vector_t &r) const {
return (isInsideTransverse(r)
&& r(2) >= -0.5 * PeriodLength_m
&& r(2) < startExitField_m);
}
}
\ No newline at end of file
......@@ -20,7 +20,6 @@ Astra1DDynamic::Astra1DDynamic(std::string aFilename):
int skippedValues = 0;
std::string tmpString;
double tmpDouble;
double tmpDouble2;
Type = TAstraDynamic;
......@@ -53,7 +52,7 @@ Astra1DDynamic::Astra1DDynamic(std::string aFilename):
parsing_passed = parsing_passed &&
interpretLine<double, double>(file, zbegin_m, tmpDouble);
tmpDouble2 = zbegin_m;
double tmpDouble2 = zbegin_m;
while(!file.eof() && parsing_passed) {
parsing_passed = interpretLine<double, double>(file, zend_m, tmpDouble, false);
if (zend_m - tmpDouble2 > 1e-10) {
......@@ -193,16 +192,13 @@ bool Astra1DDynamic::getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &
double ezp = 0.0;
double ezpp = 0.0;
double ezppp = 0.0;
double somefactor_base, somefactor;
double coskzl;
double sinkzl;
int n = 1;
for (int l = 1; l < accuracy_m ; ++ l, n += 2) {
somefactor_base = Physics::two_pi / length_m * l; // = \frac{d(kz*l)}{dz}
somefactor = 1.0;
coskzl = cos(kz * l);
sinkzl = sin(kz * l);
double somefactor_base = Physics::two_pi / length_m * l; // = \frac{d(kz*l)}{dz}
double somefactor = 1.0;
double coskzl = cos(kz * l);
double sinkzl = sin(kz * l);
ez += (FourCoefs_m[n] * coskzl - FourCoefs_m[n + 1] * sinkzl);
somefactor *= somefactor_base;
ezp += somefactor * (-FourCoefs_m[n] * sinkzl - FourCoefs_m[n + 1] * coskzl);
......
......@@ -19,13 +19,12 @@ Astra1DElectroStatic::Astra1DElectroStatic(std::string aFilename)
int skippedValues = 0;
std::string tmpString;
double tmpDouble;
double tmpDouble2;
Type = TAstraElectroStatic;
// open field map, parse it and disable element on error
file.open(Filename_m.c_str());
if(file.good()) {
if (file.good()) {
bool parsing_passed = true;
try {
parsing_passed = interpretLine<std::string, int>(file, tmpString, accuracy_m);
......@@ -48,12 +47,12 @@ Astra1DElectroStatic::Astra1DElectroStatic(std::string aFilename)
parsing_passed = parsing_passed &&
interpretLine<double, double>(file, zbegin_m, tmpDouble);
tmpDouble2 = zbegin_m;
double tmpDouble2 = zbegin_m;
while(!file.eof() && parsing_passed) {
parsing_passed = interpretLine<double, double>(file, zend_m, tmpDouble, false);
if(zend_m - tmpDouble2 > 1e-10) {
if (zend_m - tmpDouble2 > 1e-10) {
tmpDouble2 = zend_m;
} else if(parsing_passed) {
} else if (parsing_passed) {
++ skippedValues;
}
}
......@@ -61,7 +60,7 @@ Astra1DElectroStatic::Astra1DElectroStatic(std::string aFilename)
num_gridpz_m = lines_read_m - 2 - skippedValues;
lines_read_m = 0;
if(!parsing_passed && !file.eof()) {
if (!parsing_passed && !file.eof()) {
disableFieldmapWarning();
zend_m = zbegin_m - 1e-3;
throw GeneralClassicException("Astra1DElectroStatic::Astra1DElectroStatic",
......@@ -81,7 +80,7 @@ Astra1DElectroStatic::~Astra1DElectroStatic() {
}
void Astra1DElectroStatic::readMap() {
if(FourCoefs_m == NULL) {
if (FourCoefs_m == NULL) {
// declare variables and allocate memory
std::ifstream in;
......@@ -109,12 +108,12 @@ void Astra1DElectroStatic::readMap() {
in.open(Filename_m.c_str());
getLine(in, tmpString);
for(int i = 0; i < num_gridpz_m && parsing_passed; /* skip increment of i here */) {
for (int i = 0; i < num_gridpz_m && parsing_passed; /* skip increment of i here */) {
parsing_passed = interpretLine<double, double>(in, zvals[i], RealValues[i]);
// 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(std::abs(RealValues[i]) > Ez_max) {
if (zvals[i] - tmpDouble > 1e-10) {
if (std::abs(RealValues[i]) > Ez_max) {
Ez_max = std::abs(RealValues[i]);
}
tmpDouble = zvals[i];
......@@ -128,14 +127,14 @@ void Astra1DElectroStatic::readMap() {
// get equidistant sampling from the, possibly, non-equidistant sampling
// using cubic spline for this
int ii = num_gridpz_m;
for(int i = 0; i < num_gridpz_m - 1; ++ i, ++ ii) {
for (int i = 0; i < num_gridpz_m - 1; ++ i, ++ ii) {
double z = zbegin_m + dz * i;
RealValues[ii] = gsl_spline_eval(Ez_interpolant, z, Ez_accel);
}
RealValues[ii ++] = RealValues[num_gridpz_m - 1];
// prepend mirror sampling points such that field values are periodic for sure
-- ii; // ii == 2*num_gridpz_m at the moment
for(int i = 0; i < num_gridpz_m; ++ i, -- ii) {
for (int i = 0; i < num_gridpz_m; ++ i, -- ii) {
RealValues[i] = RealValues[ii];
}
......@@ -149,7 +148,7 @@ void Astra1DElectroStatic::readMap() {
// normalize to Ez_max = 1 MV/m
FourCoefs_m[0] = 1e6 * RealValues[0] / (Ez_max * num_gridpz_m); // factor 1e6 due to conversion MV/m to V/m
for(int i = 1; i < 2 * accuracy_m - 1; i++) {
for (int i = 1; i < 2 * accuracy_m - 1; i++) {
FourCoefs_m[i] = 1e6 * 2.* RealValues[i] / (Ez_max * num_gridpz_m);
}
......@@ -168,7 +167,7 @@ void Astra1DElectroStatic::readMap() {
}
void Astra1DElectroStatic::freeMap() {
if(FourCoefs_m != NULL) {
if (FourCoefs_m != NULL) {
delete[] FourCoefs_m;
FourCoefs_m = NULL;
......@@ -187,16 +186,13 @@ bool Astra1DElectroStatic::getFieldstrength(const Vector_t &R, Vector_t &E, Vect
double ezp = 0.0;
double ezpp = 0.0;
double ezppp = 0.0;
double somefactor_base, somefactor;
double coskzl;
double sinkzl;
int n = 1;
for(int l = 1; l < accuracy_m ; l++, n += 2) {
somefactor_base = Physics::two_pi / length_m * l; // = \frac{d(kz*l)}{dz}
somefactor = 1.0;
coskzl = cos(kz * l);
sinkzl = sin(kz * l);
for (int l = 1; l < accuracy_m ; l++, n += 2) {
double somefactor_base = Physics::two_pi / length_m * l; // = \frac{d(kz*l)}{dz}
double somefactor = 1.0;
double coskzl = cos(kz * l);