Code indexing in gitaly is broken and leads to code not being visible to the user. We work on the issue with highest priority.

Skip to content
Snippets Groups Projects

Opal maps

Merged Opal maps
OPAL-maps into master
Merged snuverink_j requested to merge OPAL-maps into master
Compare
and
Show latest version
2 files
+ 127
− 108
Compare changes
  • Side-by-side
  • Inline
Files
2
src/Algorithms/ThickTracker.cpp
+ 116
− 106
@@ -31,9 +31,6 @@
#include "Physics/Physics.h"
// debug info
#define PHIL_WRITE 1
//
// Class ThickTracker
// ------------------------------------------------------------------------
@@ -43,6 +40,8 @@ ThickTracker::ThickTracker(const Beamline &beamline,
bool revBeam, bool revTrack)
: Tracker(beamline, reference, revBeam, revTrack)
, hamiltonian_m(1)
, RefPartR_m(0.0)
, RefPartP_m(0.0)
, itsOpalBeamline_m(beamline.getOrigin3D(), beamline.getCoordTransformationTo())
, zstart_m(0.0)
, zstop_m(0.0)
@@ -71,6 +70,8 @@ ThickTracker::ThickTracker(const Beamline &beamline,
const int& truncOrder)
: Tracker(beamline, bunch, reference, revBeam, revTrack)
, hamiltonian_m(truncOrder)
, RefPartR_m(0.0)
, RefPartP_m(0.0)
, itsDataSink_m(&ds)
, itsOpalBeamline_m(beamline.getOrigin3D(), beamline.getCoordTransformationTo())
, zstart_m(zstart)
@@ -123,39 +124,6 @@ void ThickTracker::prepareSections() {
// void ThickTracker::updateReferenceParticle() {
// const double dt = std::min(itsBunch_m->getT(), itsBunch_m->getdT());
// const double scaleFactor = Physics::c * dt;
// Vector_t Ef(0.0), Bf(0.0);
//
// IndexMap::value_t elements = itsOpalBeamline_m.getElements(referenceToLabCSTrafo_m.transformTo(RefPartR_m));
// IndexMap::value_t::const_iterator it = elements.begin();
// const IndexMap::value_t::const_iterator end = elements.end();
//
// for (; it != end; ++ it) {
// CoordinateSystemTrafo refToLocalCSTrafo = itsOpalBeamline_m.getCSTrafoLab2Local((*it)) * referenceToLabCSTrafo_m;
//
// Vector_t localR = refToLocalCSTrafo.transformTo(RefPartR_m);
// Vector_t localP = refToLocalCSTrafo.rotateTo(RefPartP_m);
// Vector_t localE(0.0), localB(0.0);
//
// if ((*it)->applyToReferenceParticle(localR,
// localP,
// itsBunch_m->getT() - 0.5 * dt,
// localE,
// localB)) {
// *gmsg << level1 << "The reference particle hit an element" << endl;
// globalEOL_m = true;
// }
//
// Ef += refToLocalCSTrafo.rotateFrom(localE);
// Bf += refToLocalCSTrafo.rotateFrom(localB);
// }
// }
/*
//TODO complete and test fringe fields
void ThickTracker::insertFringeField(SBend* pSBend, lstruct_t& mBL,
@@ -241,14 +209,12 @@ void ThickTracker::execute() {
this->fillGaps_m();
//FIXME in local system only
RefPartR_m = Vector_t(0.0);
RefPartP_m = euclidean_norm(itsBunch_m->get_pmean_Distribution()) * Vector_t(0, 0, 1);
itsBunch_m->set_sPos(zstart_m);
#ifdef PHIL_WRITE
std::ofstream disp; // stat
disp.open("dispersion.txt");
disp << std::setprecision(10);
#endif
// IpplTimings::startTimer(mapCreation_m);
// TODO need a way to implement Initial dispersion Values
@@ -395,6 +361,30 @@ void ThickTracker::track_m()
}
ThickTracker::particle_t
ThickTracker::particleToVector_m(const Vector_t& R,
const Vector_t& P) const
{
particle_t particle;
for (int d = 0; d < 3; ++d) {
particle[2 * d] = R(d);
particle[2 *d + 1] = P(d);
}
return particle;
}
void ThickTracker::vectorToParticle_m(const particle_t& particle,
Vector_t& R,
Vector_t& P) const
{
for (int d = 0; d < 3; ++d) {
R(d) = particle[2 * d];
P(d) = particle[2 *d + 1];
}
}
void ThickTracker::write_m(const map_t& map) {
if ( Ippl::myNode() == 0 ) {
@@ -422,45 +412,59 @@ void ThickTracker::write_m(const map_t& map) {
void ThickTracker::advanceParticles_m(const map_t& map) {
double betagamma = itsBunch_m->getInitialBeta() * itsBunch_m->getInitialGamma();
particle_t particle;
for (std::size_t ip = 0; ip < itsBunch_m->getLocalNum(); ++ip) {
for (int d = 0; d < 3; ++d) {
particle[2 * d] = itsBunch_m->R[ip](d);
particle[2 *d + 1] = itsBunch_m->P[ip](d);
}
particle_t particle = this->particleToVector_m(itsBunch_m->R[ip],
itsBunch_m->P[ip]);
this->updateParticle_m(particle, map);
this->vectorToParticle_m(particle,
itsBunch_m->R[ip],
itsBunch_m->P[ip]);
}
// update reference particle
particle_t particle = this->particleToVector_m(RefPartR_m,
RefPartP_m);
this->updateParticle_m(particle, map);
this->vectorToParticle_m(particle,
RefPartR_m,
RefPartP_m);
}
//Units
double pParticle= std::sqrt( particle[1] * particle[1] +
particle[3] * particle[3] +
particle[5] * particle[5]); // [beta gamma]
particle[1] /= betagamma;
particle[3] /= betagamma;
void ThickTracker::updateParticle_m(particle_t& particle,
const map_t& map)
{
double betagamma = itsBunch_m->getInitialBeta() * itsBunch_m->getInitialGamma();
//Units
double pParticle= std::sqrt( particle[1] * particle[1] +
particle[3] * particle[3] +
particle[5] * particle[5]); // [beta gamma]
particle[1] /= betagamma;
particle[3] /= betagamma;
particle[5] = std::sqrt( 1.0 + pParticle * pParticle ) / betagamma
- 1.0 / itsBunch_m->getInitialBeta();
particle[5] = std::sqrt( 1.0 + pParticle * pParticle ) / betagamma
- 1.0 / itsBunch_m->getInitialBeta();
//Apply element map
particle = map * particle;
//Apply element map
particle = map * particle;
//Units back
particle[1] *= betagamma;
particle[3] *= betagamma;
//Units back
particle[1] *= betagamma;
particle[3] *= betagamma;
pParticle = std::sqrt( ( (particle[5] + 1.0 / itsBunch_m->getInitialBeta()) * betagamma)
* ( (particle[5] + 1./itsBunch_m->getInitialBeta() ) * betagamma) - 1.0);
pParticle = std::sqrt( ( (particle[5] + 1.0 / itsBunch_m->getInitialBeta()) * betagamma)
* ( (particle[5] + 1./itsBunch_m->getInitialBeta() ) * betagamma) - 1.0);
particle[5] = std::sqrt(pParticle * pParticle -
particle[1] * particle[1] -
particle[3] * particle[3] );
itsBunch_m->set_part(particle, ip);
}
particle[5] = std::sqrt(pParticle * pParticle -
particle[1] * particle[1] -
particle[3] * particle[3] );
}
@@ -485,42 +489,48 @@ void ThickTracker::advanceDispersion_m(fMatrix_t tempMatrix,
FMatrix<double, 1, 4> initialVal,
double pos)
{
#ifdef PHIL_WRITE
std::ofstream disp;
disp.open ("dispersion.txt",std::ios::app);
disp << std::setprecision(16);
#else
Inform msg("ThickTracker", *gmsg);
msg << "Activate PHIL_WRITE to use this Function! " << endl;
#endif
FMatrix<double, 2, 2> subx, suby;
FMatrix<double, 2, 1> subdx, subdy;
FMatrix<double, 2, 1> dxi, dyi, dx, dy;
for (int i=0; i<2; i++){
dxi[i][0]= initialVal[0][i];
dyi[i][0]= initialVal[0][i+2];
subdx[i][0]=tempMatrix[i][5]*itsBunch_m->getInitialBeta();
subdy[i][0]=tempMatrix[i+2][5]*itsBunch_m->getInitialBeta();
for (int j=0; j<2; j++){
subx[i][j]=tempMatrix[i][j];
suby[i][j]=tempMatrix[i+2][j+2];
if ( Ippl::myNode() == 0 ) {
static bool first = true;
std::ofstream out;
std::string fn = OpalData::getInstance()->getInputBasename() + ".dispersion";
if ( first ) {
first = false;
out.open(fn, std::ios::out);
} else {
out.open(fn, std::ios::app);
}
out << std::setprecision(16);
FMatrix<double, 2, 2> subx, suby;
FMatrix<double, 2, 1> subdx, subdy;
FMatrix<double, 2, 1> dxi, dyi, dx, dy;
for (int i=0; i<2; i++){
dxi[i][0]= initialVal[0][i];
dyi[i][0]= initialVal[0][i+2];
subdx[i][0]=tempMatrix[i][5]*itsBunch_m->getInitialBeta();
subdy[i][0]=tempMatrix[i+2][5]*itsBunch_m->getInitialBeta();
for (int j=0; j<2; j++){
subx[i][j]=tempMatrix[i][j];
suby[i][j]=tempMatrix[i+2][j+2];
}
}
dx= subx*dxi + subdx;
dy= suby*dyi + subdy;
out <<pos<< "\t" << dx[0][0] << "\t" << dx[0][1] << "\t" << dy[0][0] << "\t" << dy[0][1] << std::endl;
}
dx= subx*dxi + subdx;
dy= suby*dyi + subdy;
disp <<pos<< "\t" << dx[0][0] << "\t" << dx[0][1] << "\t" << dy[0][0] << "\t" << dy[0][1] << std::endl;
}