39 #include "gsl/gsl_spline.h"
40 #include "gsl/gsl_interp.h"
45 #define CHECK_BSTP_FSCANF_EOF(arg) if(arg == EOF)\
46 throw GeneralClassicException("BeamStripping::getPressureFromFile",\
47 "fscanf returned EOF at " #arg);
64 pressure_m(right.pressure_m),
65 pmapfn_m(right.pmapfn_m),
66 pscale_m(right.pscale_m),
67 temperature_m(right.temperature_m),
110 "Pressure must not be zero");
123 "Temperature must not be zero");
152 "Residual gas " +
gas_m +
" not found");
175 bool isDead = (pflag != 0);
180 const int turnnumber,
const double t,
const double tstep) {
182 bool flagNeedUpdate =
false;
186 std::pair<Vector_t, double> boundingSphere;
187 boundingSphere.first = 0.5 * (rmax + rmin);
188 boundingSphere.second =
euclidean_norm(rmax - boundingSphere.first);
197 for (
unsigned int i = 0; i < tempnum; ++i) {
198 pflag =
checkPoint(bunch->
R[i](0), bunch->
R[i](1), bunch->
R[i](2));
199 if ( (pflag != 0) && (bunch->
Bin[i] != -1) )
200 flagNeedUpdate =
true;
207 return flagNeedUpdate;
232 for (
size_t i = 0; i < bunch->
getLocalNum(); ++i) {
233 bunch->
M[i] = bunch->
getM()*1E-9;
243 *gmsg <<
"* Finalize Beam Stripping" <<
endl;
264 return "BeamStripping";
270 double rpos =
sqrt(x * x + y * y);
281 double pressure = 0.0;
284 const double rad =
sqrt(x * x + y * y);
288 const int ir = (
int)xir;
290 const double wr1 = xir - (double)ir;
292 const double wr2 = 1.0 - wr1;
294 const double tempv =
atan(y / x);
295 double tet = tempv, tet_map, xit;
297 if((x < 0) && (y >= 0)) tet =
pi + tempv;
298 else if((x < 0) && (y <= 0)) tet =
pi + tempv;
299 else if((x > 0) && (y <= 0)) tet = 2.0 *
pi + tempv;
300 else if((x == 0) && (y > 0)) tet =
pi / 2.0;
301 else if((x == 0) && (y < 0)) tet = 1.5 *
pi;
304 tet = tet /
pi * 180.0;
308 double symmetry = 1.0;
309 tet_map =
fmod(tet, 360.0 / symmetry);
312 const double wt1 = xit - (double)it;
313 const double wt2 = 1.0 - wt1;
319 int r1t1, r2t1, r1t2, r2t2;
325 r2t1 =
idx(ir + 1, it);
326 r1t2 =
idx(ir, it + 1);
327 r2t2 =
idx(ir + 1, it + 1);
344 "Pressure must not be zero");
351 for(
int i = 0; i < nrad; i++)
352 PP.
rarr[i] = rmin + i * dr;
362 *gmsg <<
"* " <<
endl;
363 *gmsg <<
"* Reading pressure field map " <<
pmapfn_m <<
endl;
367 if((f = fopen(
pmapfn_m.c_str(),
"r")) == NULL) {
369 "failed to open file '" +
pmapfn_m +
"', please check if it exists");
373 *gmsg <<
"* --- Minimal radius of measured pressure map: " <<
PP.
rmin <<
" [mm]" <<
endl;
379 *gmsg <<
"* --- Stepsize in radial direction: " <<
PP.
delr <<
" [mm]" <<
endl;
383 *gmsg <<
"* --- Minimal angle of measured pressure map: " <<
PP.
tetmin <<
" [deg]" <<
endl;
388 *gmsg <<
"* --- Stepsize in azimuthal direction: " <<
PP.
dtet <<
" [deg]" <<
endl;
391 *gmsg <<
"* --- Grid points along azimuth (ntet): " <<
PField.
ntet <<
endl;
394 *gmsg <<
"* --- Grid points along radius (nrad): " <<
PField.
nrad <<
endl;
399 *gmsg <<
"* --- Total stored grid point number: " <<
PField.
ntot <<
endl;
402 *gmsg <<
"* --- Escale factor: " <<
PP.
Pfact <<
endl;
410 *gmsg <<
"*" <<
endl;
415 #undef CHECK_BSTP_FSCANF_EOF
virtual void apply(PartBunchBase< double, 3 > *bunch, const std::pair< Vector_t, double > &boundingSphere, size_t numParticlesInSimulation=0)=0
virtual double getMaxR() const
virtual ParticleMatterInteractionHandler * getParticleMatterInteraction() const
virtual std::string getPressureMapFN() const
void initR(double rmin, double dr, int nrad)
int checkPoint(const double &x, const double &y, const double &z)
double pscale_m
stores the filename of the pressure map
double minr_m
size limits took from cyclotron
void getPressureFromFile(const double &scaleFactor)
virtual bool bends() const
virtual void visitBeamStripping(const BeamStripping &)=0
Apply the algorithm to a beam stripping interaction.
Interface for a Cyclotron.
double getPressure() const
virtual bool checkBeamStripping(PartBunchBase< double, 3 > *bunch, Cyclotron *cycl, const int turnnumber, const double t, const double tstep)
ParticleAttrib< double > Q
virtual void goOnline(const double &kineticEnergy)
#define CHECK_BSTP_FSCANF_EOF(arg)
ParticleMatterInteractionHandler * parmatint_m
virtual double getMinZ() const
std::string gas_m
parameters for BeamStripping
bool reduce(Communicate &, InputIterator, InputIterator, OutputIterator, const ReduceOp &, bool *IncludeVal=0)
virtual double getMinR() const
virtual void accept(BeamlineVisitor &) const
Apply visitor to BeamStripping.
int idx(int irad, int ktet)
constexpr double pi
The value of .
virtual double getElementLength() const
Get design length.
PETE_TUTree< FnArcTan, typename T::PETE_Expr_t > atan(const PETE_Expr< T > &l)
void setPScale(double ps)
void setStop(bool stopflag)
virtual bool getStop() const
void setTemperature(double temperature)
PartBunchBase< double, 3 > * RefPartBunch_m
void setPressureMapFN(std::string pmapfn)
double getTemperature() const
std::string getBeamStrippingShape()
double checkPressure(const double &x, const double &y)
PETE_TBTree< FnFmod, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > fmod(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
size_t getLocalNum() const
void setPressure(double pressure)
double temperature_m
a scale factor for the P-field
virtual double getPScale() const
Tps< T > sqrt(const Tps< T > &x)
Square root.
virtual ElementBase::ElementType getType() const
Get element type std::string.
T rad(T x)
Convert degrees to radians.
void setResidualGas(std::string gas)
T euclidean_norm(const Vector< T > &)
Euclidean norm.
constexpr double q_e
The elementary charge in As.
double getQ() const
Access to reference data.
virtual std::string getResidualGas() const
ParticleAttrib< int > Bin
ParticleAttrib< double > M
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B)
virtual double getMaxZ() const
Interface for a single beam element.
std::vector< double > pfld
void get_bounds(Vector_t &rmin, Vector_t &rmax)
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField)
Inform & endl(Inform &inf)
virtual void getDimensions(double &zBegin, double &zEnd) const
std::vector< double > rarr