Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Open sidebar
OPAL
src
Commits
0530b5ea
Commit
0530b5ea
authored
Apr 29, 2020
by
ext-calvo_p
Browse files
Resolve "cleanup src: remove using namespace std"
parent
e7812d0b
Changes
106
Hide whitespace changes
Inline
Side-by-side
Showing
20 changed files
with
293 additions
and
313 deletions
+293
-313
src/AbstractObjects/OpalData.cpp
src/AbstractObjects/OpalData.cpp
+8
-8
src/AbstractObjects/RangeRep.cpp
src/AbstractObjects/RangeRep.cpp
+0
-1
src/Algorithms/LieMapper.cpp
src/Algorithms/LieMapper.cpp
+5
-6
src/Algorithms/MPSplitIntegrator.cpp
src/Algorithms/MPSplitIntegrator.cpp
+3
-5
src/Algorithms/MapAnalyser.cpp
src/Algorithms/MapAnalyser.cpp
+1
-1
src/Algorithms/ParallelCyclotronTracker.cpp
src/Algorithms/ParallelCyclotronTracker.cpp
+56
-58
src/Algorithms/ThickMapper.cpp
src/Algorithms/ThickMapper.cpp
+6
-7
src/Algorithms/TransportMapper.cpp
src/Algorithms/TransportMapper.cpp
+7
-8
src/Aperture/Aperture.cpp
src/Aperture/Aperture.cpp
+52
-50
src/Aperture/Split.cpp
src/Aperture/Split.cpp
+57
-57
src/BasicActions/DumpEMFields.cpp
src/BasicActions/DumpEMFields.cpp
+11
-12
src/BasicActions/DumpFields.cpp
src/BasicActions/DumpFields.cpp
+3
-3
src/Classic/AbsBeamline/BeamStripping.cpp
src/Classic/AbsBeamline/BeamStripping.cpp
+2
-2
src/Classic/AbsBeamline/Bend2D.cpp
src/Classic/AbsBeamline/Bend2D.cpp
+54
-56
src/Classic/AbsBeamline/CCollimator.cpp
src/Classic/AbsBeamline/CCollimator.cpp
+1
-2
src/Classic/AbsBeamline/Corrector.cpp
src/Classic/AbsBeamline/Corrector.cpp
+2
-3
src/Classic/AbsBeamline/Cyclotron.cpp
src/Classic/AbsBeamline/Cyclotron.cpp
+18
-21
src/Classic/AbsBeamline/Degrader.cpp
src/Classic/AbsBeamline/Degrader.cpp
+4
-7
src/Classic/AbsBeamline/ElementBase.cpp
src/Classic/AbsBeamline/ElementBase.cpp
+1
-2
src/Classic/AbsBeamline/FlexibleCollimator.cpp
src/Classic/AbsBeamline/FlexibleCollimator.cpp
+2
-4
No files found.
src/AbstractObjects/OpalData.cpp
View file @
0530b5ea
...
...
@@ -45,7 +45,7 @@
// /DTA
#define MAX_NUM_INSTANCES 10
using
namespace
std
;
// Class OpalData::ClearReference
// ------------------------------------------------------------------------
...
...
@@ -524,8 +524,8 @@ void OpalData::define(Object *newObject) {
if
(
table
->
isDependent
(
name
))
{
if
(
Options
::
info
)
{
cerr
<<
endl
<<
"Erasing dependent table
\"
"
<<
tableName
<<
"
\"
."
<<
endl
<<
endl
;
std
::
cerr
<<
std
::
endl
<<
"Erasing dependent table
\"
"
<<
tableName
<<
"
\"
."
<<
std
::
endl
;
}
// Remove table from directory.
...
...
@@ -599,8 +599,8 @@ void OpalData::makeDirty(Object *obj) {
void
OpalData
::
printNames
(
std
::
ostream
&
os
,
const
std
::
string
&
pattern
)
{
int
column
=
0
;
RegularExpression
regex
(
pattern
);
os
<<
endl
<<
"Object names matching the pattern
\"
"
<<
pattern
<<
"
\"
:"
<<
endl
;
os
<<
std
::
endl
<<
"Object names matching the pattern
\"
"
<<
pattern
<<
"
\"
:"
<<
std
::
endl
;
for
(
ObjectDir
::
const_iterator
index
=
p
->
mainDirectory
.
begin
();
index
!=
p
->
mainDirectory
.
end
();
++
index
)
{
...
...
@@ -617,14 +617,14 @@ void OpalData::printNames(std::ostream &os, const std::string &pattern) {
column
++
;
}
while
((
column
%
20
)
!=
0
);
}
else
{
os
<<
endl
;
os
<<
std
::
endl
;
column
=
0
;
}
}
}
if
(
column
)
os
<<
endl
;
os
<<
endl
;
if
(
column
)
os
<<
std
::
endl
;
os
<<
std
::
endl
;
}
...
...
src/AbstractObjects/RangeRep.cpp
View file @
0530b5ea
...
...
@@ -19,7 +19,6 @@
#include "AbstractObjects/RangeRep.h"
#include "AbstractObjects/Element.h"
#include <iostream>
using
namespace
std
;
// Class RangeRep
...
...
src/Algorithms/LieMapper.cpp
View file @
0530b5ea
...
...
@@ -54,7 +54,6 @@
class
Beamline
;
class
PartData
;
using
Physics
::
c
;
typedef
FTps
<
double
,
6
>
Series
;
...
...
@@ -183,7 +182,7 @@ void LieMapper::visitMonitor(const Monitor &corr) {
void
LieMapper
::
visitMultipole
(
const
Multipole
&
mult
)
{
double
length
=
mult
.
getElementLength
()
*
flip_s
;
const
BMultipoleField
&
field
=
mult
.
getField
();
double
scale
=
(
flip_B
*
itsReference
.
getQ
()
*
c
)
/
itsReference
.
getP
();
double
scale
=
(
flip_B
*
itsReference
.
getQ
()
*
Physics
::
c
)
/
itsReference
.
getP
();
if
(
length
)
{
// Normal case: Finite-length multipole, field coefficients are B.
...
...
@@ -231,7 +230,7 @@ void LieMapper::visitRBend(const RBend &bend) {
// Geometry.
const
RBendGeometry
&
geometry
=
bend
.
getGeometry
();
double
length
=
flip_s
*
geometry
.
getElementLength
();
double
scale
=
(
flip_B
*
itsReference
.
getQ
()
*
c
)
/
itsReference
.
getP
();
double
scale
=
(
flip_B
*
itsReference
.
getQ
()
*
Physics
::
c
)
/
itsReference
.
getP
();
const
BMultipoleField
&
field
=
bend
.
getField
();
if
(
length
)
{
...
...
@@ -313,7 +312,7 @@ void LieMapper::visitRFCavity(const RFCavity &as) {
// Compute Hamiltonian.
static
const
Series
t
=
Series
::
makeVariable
(
AbstractMapper
::
T
);
Series
H
=
peak
*
cos
(
as
.
getPhase
()
+
(
freq
/
c
)
*
t
);
Series
H
=
peak
*
cos
(
as
.
getPhase
()
+
(
freq
/
Physics
::
c
)
*
t
);
// Build map.
DragtFinnMap
<
3
>
theMap
=
DragtFinnMap
<
3
>::
factorSimple
(
H
);
...
...
@@ -333,7 +332,7 @@ void LieMapper::visitRFQuadrupole(const RFQuadrupole &rfq) {
void
LieMapper
::
visitSBend
(
const
SBend
&
bend
)
{
const
PlanarArcGeometry
&
geometry
=
bend
.
getGeometry
();
double
length
=
flip_s
*
geometry
.
getElementLength
();
double
scale
=
(
flip_B
*
itsReference
.
getQ
()
*
c
)
/
itsReference
.
getP
();
double
scale
=
(
flip_B
*
itsReference
.
getQ
()
*
Physics
::
c
)
/
itsReference
.
getP
();
const
BMultipoleField
&
field
=
bend
.
getField
();
if
(
length
)
{
...
...
@@ -417,7 +416,7 @@ void LieMapper::visitSolenoid(const Solenoid &solenoid) {
double
length
=
flip_s
*
solenoid
.
getElementLength
();
if
(
length
)
{
double
ks
=
(
flip_B
*
itsReference
.
getQ
()
*
solenoid
.
getBz
()
*
c
)
/
double
ks
=
(
flip_B
*
itsReference
.
getQ
()
*
solenoid
.
getBz
()
*
Physics
::
c
)
/
(
2.0
*
itsReference
.
getP
());
if
(
ks
)
{
...
...
src/Algorithms/MPSplitIntegrator.cpp
View file @
0530b5ea
...
...
@@ -29,8 +29,6 @@
#include "FixedAlgebra/FVps.h"
#include "Physics/Physics.h"
using
Physics
::
c
;
// Class MPSplitIntegrator
// ------------------------------------------------------------------------
...
...
@@ -82,7 +80,7 @@ void MPSplitIntegrator::trackMap(FVps<double, 6> &map,
double
length
=
itsMultipole
->
getElementLength
();
if
(
revTrack
)
length
=
-
length
;
BMultipoleField
field
=
itsMultipole
->
getField
();
double
scale
=
(
ref
.
getQ
()
*
c
)
/
(
ref
.
getP
());
double
scale
=
(
ref
.
getQ
()
*
Physics
::
c
)
/
(
ref
.
getP
());
if
(
revBeam
)
scale
=
-
scale
;
if
(
length
)
{
...
...
@@ -112,7 +110,7 @@ void MPSplitIntegrator::trackParticle(OpalParticle &part, const PartData &ref,
double
length
=
itsMultipole
->
getElementLength
();
if
(
revTrack
)
length
=
-
length
;
BMultipoleField
field
=
itsMultipole
->
getField
();
double
scale
=
(
ref
.
getQ
()
*
c
)
/
(
ref
.
getP
());
double
scale
=
(
ref
.
getQ
()
*
Physics
::
c
)
/
(
ref
.
getP
());
if
(
revBeam
)
scale
=
-
scale
;
if
(
length
)
{
...
...
@@ -143,7 +141,7 @@ void MPSplitIntegrator::trackBunch(PartBunchBase<double, 3> *bunch,
double
length
=
itsMultipole
->
getElementLength
();
if
(
revTrack
)
length
=
-
length
;
BMultipoleField
field
=
itsMultipole
->
getField
();
double
scale
=
(
ref
.
getQ
()
*
c
)
/
(
ref
.
getP
());
double
scale
=
(
ref
.
getQ
()
*
Physics
::
c
)
/
(
ref
.
getP
());
if
(
revBeam
)
scale
=
-
scale
;
if
(
length
)
{
...
...
src/Algorithms/MapAnalyser.cpp
View file @
0530b5ea
...
...
@@ -306,7 +306,7 @@ void MapAnalyser::normalizeEigen_m(cfMatrix_t& eigenVecM, cfMatrix_t& invEigenVe
for
(
int
j
=
0
;
j
<
6
;
j
+=
2
){
temp
+=
2
*
(
eigenVecM
[
j
][
i
]
*
std
::
conj
(
eigenVecM
[
j
+
1
][
i
])).
imag
();
}
temp
=
std
::
f
abs
(
temp
);
temp
=
std
::
abs
(
temp
);
if
(
temp
>
1e-10
){
for
(
int
j
=
0
;
j
<
6
;
j
++
){
...
...
src/Algorithms/ParallelCyclotronTracker.cpp
View file @
0530b5ea
...
...
@@ -23,6 +23,7 @@
#include <limits>
#include <vector>
#include <numeric>
#include <cmath>
#include "AbstractObjects/Element.h"
#include "AbstractObjects/OpalData.h"
...
...
@@ -82,9 +83,6 @@
#include "Structure/DataSink.h"
#include "Structure/LossDataSink.h"
using
Physics
::
pi
;
using
Physics
::
q_e
;
constexpr
double
c_mmtns
=
Physics
::
c
*
1.0e-6
;
// m/s --> mm/ns
...
...
@@ -348,13 +346,13 @@ void ParallelCyclotronTracker::visitRing(const Ring &ring) {
referencePz
=
0.0
;
referencePtot
=
itsReference
.
getGamma
()
*
itsReference
.
getBeta
();
referencePt
=
sqrt
(
referencePtot
*
referencePtot
-
referencePr
*
referencePr
);
referencePt
=
std
::
sqrt
(
referencePtot
*
referencePtot
-
referencePr
*
referencePr
);
if
(
referencePtot
<
0.0
)
referencePt
*=
-
1.0
;
sinRefTheta_m
=
sin
(
referenceTheta
*
Physics
::
deg2rad
);
cosRefTheta_m
=
cos
(
referenceTheta
*
Physics
::
deg2rad
);
sinRefTheta_m
=
std
::
sin
(
referenceTheta
*
Physics
::
deg2rad
);
cosRefTheta_m
=
std
::
cos
(
referenceTheta
*
Physics
::
deg2rad
);
double
BcParameter
[
8
]
=
{};
// zero initialise array
...
...
@@ -447,7 +445,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
}
else
{
referencePt
=
sqrt
(
insqrt
);
referencePt
=
std
::
sqrt
(
insqrt
);
}
if
(
referencePtot
<
0.0
)
...
...
@@ -483,8 +481,8 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
}
}
sinRefTheta_m
=
sin
(
referenceTheta
*
Physics
::
deg2rad
);
cosRefTheta_m
=
cos
(
referenceTheta
*
Physics
::
deg2rad
);
sinRefTheta_m
=
std
::
sin
(
referenceTheta
*
Physics
::
deg2rad
);
cosRefTheta_m
=
std
::
cos
(
referenceTheta
*
Physics
::
deg2rad
);
*
gmsg
<<
endl
;
*
gmsg
<<
"* Bunch global starting position:"
<<
endl
;
...
...
@@ -1366,7 +1364,7 @@ void ParallelCyclotronTracker::MtsTracker() {
// INFOMSG("No space charge Effects are included!"<<endl;);
if
((
step_m
%
Options
::
repartFreq
*
100
)
==
0
)
{
//TODO: why * 100?
Vector_t
const
meanP
=
calcMeanP
();
double
const
phi
=
calculateAngle
(
meanP
(
0
),
meanP
(
1
))
-
0.5
*
pi
;
double
const
phi
=
calculateAngle
(
meanP
(
0
),
meanP
(
1
))
-
0.5
*
Physics
::
pi
;
Vector_t
const
meanR
=
calcMeanR
();
globalToLocal
(
itsBunch_m
->
R
,
phi
,
meanR
);
itsBunch_m
->
updateNumTotal
();
...
...
@@ -1571,8 +1569,8 @@ bool ParallelCyclotronTracker::checkGapCross(Vector_t Rold, Vector_t Rnew,
bool
ParallelCyclotronTracker
::
RFkick
(
RFCavity
*
rfcavity
,
const
double
t
,
const
double
dt
,
const
int
Pindex
){
// For OPAL 2.0: As long as the RFCavity is in mm, we have to change R to mm here -DW
double
radius
=
s
qrt
(
pow
(
1000.0
*
itsBunch_m
->
R
[
Pindex
](
0
),
2.0
)
+
pow
(
1000.0
*
itsBunch_m
->
R
[
Pindex
](
1
),
2.0
)
-
pow
(
rfcavity
->
getPerpenDistance
()
,
2.0
));
double
radius
=
s
td
::
sqrt
(
std
::
pow
(
1000.0
*
itsBunch_m
->
R
[
Pindex
](
0
),
2.0
)
+
std
::
pow
(
1000.0
*
itsBunch_m
->
R
[
Pindex
](
1
),
2.0
)
-
std
::
pow
(
rfcavity
->
getPerpenDistance
()
,
2.0
));
double
rmin
=
rfcavity
->
getRmin
();
double
rmax
=
rfcavity
->
getRmax
();
double
nomalRadius
=
(
radius
-
rmin
)
/
(
rmax
-
rmin
);
...
...
@@ -1740,8 +1738,8 @@ void ParallelCyclotronTracker::globalToLocal(ParticleAttrib<Vector_t> & particle
IpplTimings
::
startTimer
(
TransformTimer_m
);
particleVectors
-=
translationToGlobal
;
Tenzor
<
double
,
3
>
const
rotation
(
cos
(
phi
),
sin
(
phi
),
0
,
-
sin
(
phi
),
cos
(
phi
),
0
,
Tenzor
<
double
,
3
>
const
rotation
(
std
::
cos
(
phi
),
std
::
sin
(
phi
),
0
,
-
std
::
sin
(
phi
),
std
::
cos
(
phi
),
0
,
0
,
0
,
1
);
// clockwise rotation
for
(
unsigned
int
i
=
0
;
i
<
itsBunch_m
->
getLocalNum
();
++
i
)
{
...
...
@@ -1754,8 +1752,8 @@ void ParallelCyclotronTracker::globalToLocal(ParticleAttrib<Vector_t> & particle
void
ParallelCyclotronTracker
::
localToGlobal
(
ParticleAttrib
<
Vector_t
>
&
particleVectors
,
double
phi
,
Vector_t
const
translationToGlobal
)
{
IpplTimings
::
startTimer
(
TransformTimer_m
);
Tenzor
<
double
,
3
>
const
rotation
(
cos
(
phi
),
-
sin
(
phi
),
0
,
sin
(
phi
),
cos
(
phi
),
0
,
Tenzor
<
double
,
3
>
const
rotation
(
std
::
cos
(
phi
),
-
std
::
sin
(
phi
),
0
,
std
::
sin
(
phi
),
std
::
cos
(
phi
),
0
,
0
,
0
,
1
);
// counter-clockwise rotation
for
(
unsigned
int
i
=
0
;
i
<
itsBunch_m
->
getLocalNum
();
++
i
)
{
...
...
@@ -1900,9 +1898,9 @@ inline void ParallelCyclotronTracker::normalizeQuaternion(Quaternion_t & quatern
double
tolerance
=
1.0e-10
;
double
length2
=
dot
(
quaternion
,
quaternion
);
if
(
f
abs
(
length2
)
>
tolerance
&&
f
abs
(
length2
-
1.0
f
)
>
tolerance
)
{
if
(
std
::
abs
(
length2
)
>
tolerance
&&
std
::
abs
(
length2
-
1.0
f
)
>
tolerance
)
{
double
length
=
sqrt
(
length2
);
double
length
=
std
::
sqrt
(
length2
);
quaternion
/=
length
;
}
}
...
...
@@ -1912,9 +1910,9 @@ inline void ParallelCyclotronTracker::normalizeVector(Vector_t & vector) {
double
tolerance
=
1.0e-10
;
double
length2
=
dot
(
vector
,
vector
);
if
(
f
abs
(
length2
)
>
tolerance
&&
f
abs
(
length2
-
1.0
f
)
>
tolerance
)
{
if
(
std
::
abs
(
length2
)
>
tolerance
&&
std
::
abs
(
length2
-
1.0
f
)
>
tolerance
)
{
double
length
=
sqrt
(
length2
);
double
length
=
std
::
sqrt
(
length2
);
vector
/=
length
;
}
}
...
...
@@ -1922,8 +1920,8 @@ inline void ParallelCyclotronTracker::normalizeVector(Vector_t & vector) {
inline
void
ParallelCyclotronTracker
::
rotateAroundZ
(
ParticleAttrib
<
Vector_t
>
&
particleVectors
,
double
const
phi
)
{
// Clockwise rotation of particles 'particleVectors' by 'phi' around Z axis
Tenzor
<
double
,
3
>
const
rotation
(
cos
(
phi
),
sin
(
phi
),
0
,
-
sin
(
phi
),
cos
(
phi
),
0
,
Tenzor
<
double
,
3
>
const
rotation
(
std
::
cos
(
phi
),
std
::
sin
(
phi
),
0
,
-
std
::
sin
(
phi
),
std
::
cos
(
phi
),
0
,
0
,
0
,
1
);
for
(
unsigned
int
i
=
0
;
i
<
itsBunch_m
->
getLocalNum
();
++
i
)
{
...
...
@@ -1935,8 +1933,8 @@ inline void ParallelCyclotronTracker::rotateAroundZ(ParticleAttrib<Vector_t> & p
inline
void
ParallelCyclotronTracker
::
rotateAroundZ
(
Vector_t
&
myVector
,
double
const
phi
)
{
// Clockwise rotation of single vector 'myVector' by 'phi' around Z axis
Tenzor
<
double
,
3
>
const
rotation
(
cos
(
phi
),
sin
(
phi
),
0
,
-
sin
(
phi
),
cos
(
phi
),
0
,
Tenzor
<
double
,
3
>
const
rotation
(
std
::
cos
(
phi
),
std
::
sin
(
phi
),
0
,
-
std
::
sin
(
phi
),
std
::
cos
(
phi
),
0
,
0
,
0
,
1
);
myVector
=
dot
(
rotation
,
myVector
);
...
...
@@ -1946,8 +1944,8 @@ inline void ParallelCyclotronTracker::rotateAroundX(ParticleAttrib<Vector_t> & p
// Clockwise rotation of particles 'particleVectors' by 'psi' around X axis
Tenzor
<
double
,
3
>
const
rotation
(
1
,
0
,
0
,
0
,
cos
(
psi
),
sin
(
psi
),
0
,
-
sin
(
psi
),
cos
(
psi
));
0
,
std
::
cos
(
psi
),
std
::
sin
(
psi
),
0
,
-
std
::
sin
(
psi
),
std
::
cos
(
psi
));
for
(
unsigned
int
i
=
0
;
i
<
itsBunch_m
->
getLocalNum
();
++
i
)
{
...
...
@@ -1959,8 +1957,8 @@ inline void ParallelCyclotronTracker::rotateAroundX(Vector_t & myVector, double
// Clockwise rotation of single vector 'myVector' by 'psi' around X axis
Tenzor
<
double
,
3
>
const
rotation
(
1
,
0
,
0
,
0
,
cos
(
psi
),
sin
(
psi
),
0
,
-
sin
(
psi
),
cos
(
psi
));
0
,
std
::
cos
(
psi
),
std
::
sin
(
psi
),
0
,
-
std
::
sin
(
psi
),
std
::
cos
(
psi
));
myVector
=
dot
(
rotation
,
myVector
);
}
...
...
@@ -1972,12 +1970,12 @@ inline void ParallelCyclotronTracker::getQuaternionTwoVectors(Vector_t u, Vector
normalizeVector
(
v
);
double
k_cos_theta
=
dot
(
u
,
v
);
double
k
=
sqrt
(
dot
(
u
,
u
)
*
dot
(
v
,
v
));
double
k
=
std
::
sqrt
(
dot
(
u
,
u
)
*
dot
(
v
,
v
));
double
tolerance1
=
1.0e-5
;
double
tolerance2
=
1.0e-8
;
Vector_t
resultVectorComponent
;
if
(
f
abs
(
k_cos_theta
/
k
+
1.0
)
<
tolerance1
)
{
if
(
std
::
abs
(
k_cos_theta
/
k
+
1.0
)
<
tolerance1
)
{
// u and v are almost exactly antiparallel so we need to do
// 180 degree rotation around any vector orthogonal to u
...
...
@@ -1989,13 +1987,13 @@ inline void ParallelCyclotronTracker::getQuaternionTwoVectors(Vector_t u, Vector
resultVectorComponent
=
cross
(
u
,
zaxis
);
}
double
halfAngle
=
0.5
*
pi
;
double
sinHalfAngle
=
sin
(
halfAngle
);
double
halfAngle
=
0.5
*
Physics
::
pi
;
double
sinHalfAngle
=
std
::
sin
(
halfAngle
);
resultVectorComponent
*=
sinHalfAngle
;
k
=
0.0
;
k_cos_theta
=
cos
(
halfAngle
);
k_cos_theta
=
std
::
cos
(
halfAngle
);
}
else
{
...
...
@@ -2025,7 +2023,7 @@ bool ParallelCyclotronTracker::push(double h) {
bool
flagNeedUpdate
=
false
;
for
(
unsigned
int
i
=
0
;
i
<
itsBunch_m
->
getLocalNum
();
++
i
)
{
Vector_t
const
oldR
=
itsBunch_m
->
R
[
i
];
double
const
gamma
=
sqrt
(
1.0
+
dot
(
itsBunch_m
->
P
[
i
],
itsBunch_m
->
P
[
i
]));
double
const
gamma
=
std
::
sqrt
(
1.0
+
dot
(
itsBunch_m
->
P
[
i
],
itsBunch_m
->
P
[
i
]));
double
const
c_gamma
=
Physics
::
c
/
gamma
;
Vector_t
const
v
=
itsBunch_m
->
P
[
i
]
*
c_gamma
;
itsBunch_m
->
R
[
i
]
+=
h
*
v
;
...
...
@@ -2038,7 +2036,7 @@ bool ParallelCyclotronTracker::push(double h) {
if
(
distOld
>
0.0
)
tagCrossing
=
true
;
}
if
(
tagCrossing
)
{
double
const
dt1
=
distOld
/
sqrt
(
dot
(
v
,
v
));
double
const
dt1
=
distOld
/
std
::
sqrt
(
dot
(
v
,
v
));
double
const
dt2
=
h
-
dt1
;
// Retrack particle from the old postion to cavity gap point
...
...
@@ -2065,7 +2063,7 @@ bool ParallelCyclotronTracker::kick(double h) {
bool
flagNeedUpdate
=
false
;
BorisPusher
pusher
;
double
const
q
=
itsBunch_m
->
Q
[
0
]
/
q_e
;
// For now all particles have the same charge
double
const
q
=
itsBunch_m
->
Q
[
0
]
/
Physics
::
q_e
;
// For now all particles have the same charge
double
const
M
=
itsBunch_m
->
M
[
0
]
*
1.0e9
;
// For now all particles have the same rest energy
for
(
unsigned
int
i
=
0
;
i
<
itsBunch_m
->
getLocalNum
();
++
i
)
{
...
...
@@ -2246,10 +2244,10 @@ bool ParallelCyclotronTracker::deleteParticle(bool flagNeedUpdate){
Vector_t
const
meanP
=
calcMeanP
();
// Bunch (local) azimuth at meanR w.r.t. y-axis
double
const
phi
=
calculateAngle
(
meanP
(
0
),
meanP
(
1
))
-
0.5
*
pi
;
double
const
phi
=
calculateAngle
(
meanP
(
0
),
meanP
(
1
))
-
0.5
*
Physics
::
pi
;
// Bunch (local) elevation at meanR w.r.t. xy plane
double
const
psi
=
0.5
*
pi
-
acos
(
meanP
(
2
)
/
sqrt
(
dot
(
meanP
,
meanP
)));
double
const
psi
=
0.5
*
Physics
::
pi
-
std
::
acos
(
meanP
(
2
)
/
std
::
sqrt
(
dot
(
meanP
,
meanP
)));
// For statistics data, the bunch is transformed into a local coordinate system
// with meanP in direction of y-axis -DW
...
...
@@ -2396,10 +2394,10 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
Vector_t
const
meanP
=
calcMeanP
();
// Bunch (local) azimuth at meanR w.r.t. y-axis
double
const
phi
=
calculateAngle
(
meanP
(
0
),
meanP
(
1
))
-
0.5
*
pi
;
double
const
phi
=
calculateAngle
(
meanP
(
0
),
meanP
(
1
))
-
0.5
*
Physics
::
pi
;
// Bunch (local) elevation at meanR w.r.t. xy plane
double
const
psi
=
0.5
*
pi
-
acos
(
meanP
(
2
)
/
sqrt
(
dot
(
meanP
,
meanP
)));
double
const
psi
=
0.5
*
Physics
::
pi
-
std
::
acos
(
meanP
(
2
)
/
std
::
sqrt
(
dot
(
meanP
,
meanP
)));
double
radius
=
std
::
sqrt
(
meanR
[
0
]
*
meanR
[
0
]
+
meanR
[
1
]
*
meanR
[
1
]);
// [m]
...
...
@@ -2606,10 +2604,10 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
Vector_t
meanP
=
calcMeanP
();
// Bunch (local) azimuth at meanR w.r.t. y-axis
phi
=
calculateAngle
(
meanP
(
0
),
meanP
(
1
))
-
0.5
*
pi
;
phi
=
calculateAngle
(
meanP
(
0
),
meanP
(
1
))
-
0.5
*
Physics
::
pi
;
// Bunch (local) elevation at meanR w.r.t. xy plane
psi
=
0.5
*
pi
-
acos
(
meanP
(
2
)
/
sqrt
(
dot
(
meanP
,
meanP
)));
psi
=
0.5
*
Physics
::
pi
-
std
::
acos
(
meanP
(
2
)
/
std
::
sqrt
(
dot
(
meanP
,
meanP
)));
// Rotate so Pmean is in positive y direction. No shift, so that normalized emittance and
// unnormalized emittance as well as centroids are calculated correctly in
...
...
@@ -2662,10 +2660,10 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
if
(
Options
::
psDumpFrame
!=
Options
::
GLOBAL
)
{
// -------------------- ----------- Do Transformations ---------------------------------- //
// Bunch (local) azimuth at meanR w.r.t. y-axis
phi
=
calculateAngle
(
meanP
(
0
),
meanP
(
1
))
-
0.5
*
pi
;
phi
=
calculateAngle
(
meanP
(
0
),
meanP
(
1
))
-
0.5
*
Physics
::
pi
;
// Bunch (local) elevation at meanR w.r.t. xy plane
psi
=
0.5
*
pi
-
acos
(
meanP
(
2
)
/
sqrt
(
dot
(
meanP
,
meanP
)));
psi
=
0.5
*
Physics
::
pi
-
std
::
acos
(
meanP
(
2
)
/
std
::
sqrt
(
dot
(
meanP
,
meanP
)));
// Rotate so Pmean is in positive y direction. No shift, so that normalized emittance and
// unnormalized emittance as well as centroids are calculated correctly in
...
...
@@ -2720,18 +2718,18 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
meanP
=
itsBunch_m
->
P
[
0
];
}
double
const
betagamma_temp
=
sqrt
(
dot
(
meanP
,
meanP
));
double
const
betagamma_temp
=
std
::
sqrt
(
dot
(
meanP
,
meanP
));
double
const
E
=
itsBunch_m
->
get_meanKineticEnergy
();
// Bunch (global) angle w.r.t. x-axis (cylinder coordinates)
double
const
theta
=
atan2
(
meanR
(
1
),
meanR
(
0
));
double
const
theta
=
std
::
atan2
(
meanR
(
1
),
meanR
(
0
));
// Bunch (local) azimuth at meanR w.r.t. y-axis
double
const
phi
=
calculateAngle
(
meanP
(
0
),
meanP
(
1
))
-
0.5
*
pi
;
double
const
phi
=
calculateAngle
(
meanP
(
0
),
meanP
(
1
))
-
0.5
*
Physics
::
pi
;
// Bunch (local) elevation at meanR w.r.t. xy plane
double
const
psi
=
0.5
*
pi
-
acos
(
meanP
(
2
)
/
sqrt
(
dot
(
meanP
,
meanP
)));
double
const
psi
=
0.5
*
Physics
::
pi
-
std
::
acos
(
meanP
(
2
)
/
std
::
sqrt
(
dot
(
meanP
,
meanP
)));
// ---------------- Re-calculate reference values in format of input values ----------------- //
// Position:
...
...
@@ -2743,8 +2741,8 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
// Momentum in Theta-hat, R-hat, Z-hat coordinates at position meanR:
referencePtot
=
betagamma_temp
;
referencePz
=
meanP
(
2
);
referencePr
=
meanP
(
0
)
*
cos
(
theta
)
+
meanP
(
1
)
*
sin
(
theta
);
referencePt
=
sqrt
(
referencePtot
*
referencePtot
-
\
referencePr
=
meanP
(
0
)
*
std
::
cos
(
theta
)
+
meanP
(
1
)
*
std
::
sin
(
theta
);
referencePt
=
std
::
sqrt
(
referencePtot
*
referencePtot
-
\
referencePz
*
referencePz
-
referencePr
*
referencePr
);
// ----- Calculate the external fields at the center of the bunch (Cave: Global Frame) ----- //
...
...
@@ -2890,7 +2888,7 @@ std::tuple<double, double, double> ParallelCyclotronTracker::initializeTracking_
double
t
=
itsBunch_m
->
getT
()
*
1.0e9
;
// current time (s --> ns)
double
oldReferenceTheta
=
referenceTheta
*
Physics
::
deg2rad
;
// init here, reset each step
setup_m
.
deltaTheta
=
pi
/
(
setup_m
.
stepsPerTurn
);
// half of the average angle per step
setup_m
.
deltaTheta
=
Physics
::
pi
/
(
setup_m
.
stepsPerTurn
);
// half of the average angle per step
//int stepToLastInj = itsBunch_m->getSteptoLastInj(); // TODO: Do we need this? -DW
...
...
@@ -2988,8 +2986,8 @@ void ParallelCyclotronTracker::finalizeTracking_m(dvector_t& Ttime,
{
for
(
size_t
ii
=
0
;
ii
<
(
itsBunch_m
->
getLocalNum
());
ii
++
)
{
if
(
itsBunch_m
->
ID
[
ii
]
==
0
)
{
double
FinalMomentum2
=
pow
(
itsBunch_m
->
P
[
ii
](
0
),
2.0
)
+
pow
(
itsBunch_m
->
P
[
ii
](
1
),
2.0
)
+
pow
(
itsBunch_m
->
P
[
ii
](
2
),
2.0
);
double
FinalEnergy
=
(
sqrt
(
1.0
+
FinalMomentum2
)
-
1.0
)
*
itsBunch_m
->
getM
()
*
1.0e-6
;
double
FinalMomentum2
=
std
::
pow
(
itsBunch_m
->
P
[
ii
](
0
),
2.0
)
+
std
::
pow
(
itsBunch_m
->
P
[
ii
](
1
),
2.0
)
+
std
::
pow
(
itsBunch_m
->
P
[
ii
](
2
),
2.0
);
double
FinalEnergy
=
(
std
::
sqrt
(
1.0
+
FinalMomentum2
)
-
1.0
)
*
itsBunch_m
->
getM
()
*
1.0e-6
;
*
gmsg
<<
"* Final energy of reference particle = "
<<
FinalEnergy
<<
" [MeV]"
<<
endl
;
*
gmsg
<<
"* Total phase space dump number(includes the initial distribution) = "
<<
lastDumpedStep_m
+
1
<<
endl
;
*
gmsg
<<
"* One can restart simulation from the last dump step (--restart "
<<
lastDumpedStep_m
<<
")"
<<
endl
;
...
...
@@ -3004,7 +3002,7 @@ void ParallelCyclotronTracker::finalizeTracking_m(dvector_t& Ttime,
{
// Calculate tunes after tracking.
*
gmsg
<<
endl
;
*
gmsg
<<
"* **************** The result for tune calulation (NO space charge) ******************* *"
<<
endl
*
gmsg
<<
"* **************** The result for tune cal
c
ulation (NO space charge) ******************* *"
<<
endl
<<
"* Number of tracked turns: "
<<
TturnNumber
.
back
()
<<
endl
;
double
nur
,
nuz
;
getTunes
(
Ttime
,
Tdeltr
,
Tdeltz
,
TturnNumber
.
back
(),
nur
,
nuz
);
...
...
@@ -3072,8 +3070,8 @@ void ParallelCyclotronTracker::seoMode_m(double& t, const double dt, bool& /*fin
}
double
OldTheta
=
calculateAngle
(
itsBunch_m
->
R
[
i
](
0
),
itsBunch_m
->
R
[
i
](
1
));
r_tuning
[
i
]
=
itsBunch_m
->
R
[
i
](
0
)
*
cos
(
OldTheta
)
+
itsBunch_m
->
R
[
i
](
1
)
*
sin
(
OldTheta
);
r_tuning
[
i
]
=
itsBunch_m
->
R
[
i
](
0
)
*
std
::
cos
(
OldTheta
)
+
itsBunch_m
->
R
[
i
](
1
)
*
std
::
sin
(
OldTheta
);
z_tuning
[
i
]
=
itsBunch_m
->
R
[
i
](
2
);
...
...
@@ -3281,8 +3279,8 @@ void ParallelCyclotronTracker::gapCrossKick_m(size_t i, double t,
itsBunch_m
->
cavityGapCrossed
[
i
]
=
true
;
double
oldMomentum2
=
dot
(
Pold
,
Pold
);
double
oldBetgam
=
sqrt
(
oldMomentum2
);
double
oldGamma
=
sqrt
(
1.0
+
oldMomentum2
);
double
oldBetgam
=
std
::
sqrt
(
oldMomentum2
);
double
oldGamma
=
std
::
sqrt
(
1.0
+
oldMomentum2
);
double
oldBeta
=
oldBetgam
/
oldGamma
;
double
dt1
=
DistOld
/
(
Physics
::
c
*
oldBeta
*
1.0e-6
);
// ns
double
dt2
=
dt
-
dt1
;
...
...
src/Algorithms/ThickMapper.cpp
View file @
0530b5ea
...
...
@@ -55,7 +55,6 @@
class
Beamline
;
class
PartData
;
using
Physics
::
c
;
#define PSdim 6
typedef
FVector
<
double
,
PSdim
>
Vector
;
...
...
@@ -106,7 +105,7 @@ void ThickMapper::visitCorrector(const Corrector &corr) {
// Apply kick.
double
scale
=
(
flip_s
*
flip_B
*
itsReference
.
getQ
()
*
c
)
/
itsReference
.
getP
();
(
flip_s
*
flip_B
*
itsReference
.
getQ
()
*
Physics
::
c
)
/
itsReference
.
getP
();
const
BDipoleField
&
field
=
corr
.
getField
();
itsMap
[
PX
]
-=
field
.
getBy
()
*
scale
;
itsMap
[
PY
]
+=
field
.
getBx
()
*
scale
;
...
...
@@ -153,7 +152,7 @@ void ThickMapper::visitMultipole(const Multipole &mult) {
//std::cerr << "==> In ThickMapper::visitMultipole(const Multipole &mult)" << std::endl;
// Geometry and Field
double
length
=
flip_s
*
mult
.
getElementLength
();
double
scale
=
(
flip_B
*
itsReference
.
getQ
()
*
c
)
/
itsReference
.
getP
();
double
scale
=
(
flip_B
*
itsReference
.
getQ
()
*
Physics
::
c
)
/
itsReference
.
getP
();
const
BMultipoleField
&
field
=
mult
.
getField
();
int
order
=
field
.
order
();
...
...
@@ -220,7 +219,7 @@ void ThickMapper::visitRBend(const RBend &bend) {
// Geometry and Field.
const
RBendGeometry
&
geometry
=
bend
.
getGeometry
();
double
length
=
flip_s
*
geometry
.
getElementLength
();
double
scale
=
(
flip_B
*
itsReference
.
getQ
()
*
c
)
/
itsReference
.
getP
();
double
scale
=
(
flip_B
*
itsReference
.
getQ
()
*
Physics
::
c
)
/
itsReference
.
getP
();
const
BMultipoleField
&
field
=
bend
.
getField
();
int
order
=
field
.
order
();
double
beta_inv
=
1.0
/
itsReference
.
getBeta
();
...
...
@@ -336,7 +335,7 @@ void ThickMapper::visitRFCavity(const RFCavity &as) {
double
peak
=
flip_s
*
as
.
getAmplitude
()
/
itsReference
.
getP
();
Series
pt
=
itsMap
[
PT
]
+
1.0
;
Series
speed
=
(
c
*
pt
)
/
sqrt
(
pt
*
pt
+
kin
*
kin
);
Series
speed
=
(
Physics
::
c
*
pt
)
/
sqrt
(
pt
*
pt
+
kin
*
kin
);
Series
phase
=
as
.
getPhase
()
+
freq
*
itsMap
[
T
]
/
speed
;
itsMap
[
PT
]
+=
peak
*
sin
(
phase
)
/
pt
;
...
...
@@ -356,7 +355,7 @@ void ThickMapper::visitSBend(const SBend &bend) {
// Geometry and Field.
const
PlanarArcGeometry
&
geometry
=
bend
.
getGeometry
();
double
length
=
flip_s
*
geometry
.
getElementLength
();
double
scale
=
(
flip_B
*
itsReference
.
getQ
()
*
c
)
/
itsReference
.
getP
();
double
scale
=
(
flip_B
*
itsReference
.
getQ
()
*
Physics
::
c
)
/
itsReference
.
getP
();
const
BMultipoleField
&
field
=
bend
.
getField
();
int
order
=
field
.
order
();
double
beta_inv
=
1.0
/
itsReference
.
getBeta
();
...
...
@@ -480,7 +479,7 @@ void ThickMapper::visitSolenoid(const Solenoid &solenoid) {
double
length
=
flip_s
*
solenoid
.
getElementLength
();
if
(
length
!=
0.0
)
{
double
ks
=
(
flip_B
*
itsReference
.
getQ
()
*
solenoid
.
getBz
()
*
c
)
/
double
ks
=
(
flip_B
*
itsReference
.
getQ
()
*
solenoid
.
getBz
()
*
Physics
::
c
)
/
(
2.0
*
itsReference
.
getP
());
if
(
ks
)
{
...
...
src/Algorithms/TransportMapper.cpp
View file @
0530b5ea
...
...
@@ -51,9 +51,8 @@
#include "FixedAlgebra/TransportFun.h"
#include "FixedAlgebra/TransportMath.h"
#include "Physics/Physics.h"
#include <cmath>
using
Physics
::
c
;
#include <cmath>
typedef
FTps
<
double
,
2
>
Series2
;
typedef
TransportFun
<
double
,
6
>
TptFun
;
...
...
@@ -141,7 +140,7 @@ void TransportMapper::visitCorrector(const Corrector &corr) {
if
(
length
)
applyDrift
(
length
/
2.0
);
// Apply kick.
double
scale
=
(
flip_B
*
itsReference
.
getQ
()
*
c
)
/
itsReference
.
getP
();
double
scale
=
(
flip_B
*
itsReference
.
getQ
()
*
Physics
::
c
)
/
itsReference
.
getP
();
const
BDipoleField
&
field
=
corr
.
getField
();
itsMap
[
PX
]
-=
field
.
getBy
()
*
scale
;