Opal maps

This MR will implement Map tracking for:

• Quadrupole
• Dipole
• Drifts
• Bends

TODO list:

• Code implementation (@ganz_p): ThickTracker, Hamiltonian, MapAnalyser
• Doxygen documentation (@ganz_p)
• Documentation (@ganz_p)
• Code cleaning (@frey_m)
• Regression tests (@snuverink_j)
• Examples (@ganz_p)
• Jupyter notebook (@ganz_p)
• Use beam line order

Remaining issues:

• Output dumps in local coordinates only
• Dipole until 1st order
• Tune calculation not stable
• Space charge calculations needs to be verified
• Fringe fields for Bends

src/Algorithms/ThickTracker.cpp

@@ -40,21 +40,23 @@
 //
 ThickTracker::ThickTracker(const Beamline &beamline,
                            const PartData &reference,
-                           bool revBeam, bool revTrack):
-  Tracker(beamline, reference, revBeam, revTrack),
-  hamiltonian_m(1),
-  itsOpalBeamline_m(beamline.getOrigin3D(), beamline.getCoordTransformationTo()),
-  zstart_m(0.0),
-  threshold_m(100.0 * std::numeric_limits<double>::epsilon()),
-  dtCurrentTrack_m(0.0),
-  truncOrder_m(1), // linear
-  mapCreation_m(IpplTimings::getTimer("mapCreation")),
-  mapCombination_m(IpplTimings::getTimer("mapCombination")),
-  mapTracking_m(IpplTimings::getTimer("mapTracking"))
+                           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)
+    , threshold_m(100.0 * std::numeric_limits<double>::epsilon())
+    , truncOrder_m(1) // linear
+    , mapCreation_m(IpplTimings::getTimer("mapCreation"))
+    , mapCombination_m(IpplTimings::getTimer("mapCombination"))
+    , mapTracking_m(IpplTimings::getTimer("mapTracking"))
 {
-  CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
-                                 beamline.getCoordTransformationTo());
-  referenceToLabCSTrafo_m = labToRef.inverted();
+    CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
+                                   beamline.getCoordTransformationTo());
+    referenceToLabCSTrafo_m = labToRef.inverted();
 }
 
 
@@ -71,29 +73,26 @@ 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)
+    , zstop_m(zstop[0])
     , threshold_m(1.0e-6)
-    , dtCurrentTrack_m(0.0)
     , truncOrder_m(truncOrder)
     , mapCreation_m(IpplTimings::getTimer("mapCreation"))
     , mapCombination_m(IpplTimings::getTimer("mapCombination"))
     , mapTracking_m(IpplTimings::getTimer("mapTracking"))
 {
-  CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
-                                 beamline.getCoordTransformationTo());
-  referenceToLabCSTrafo_m = labToRef.inverted();
-
-  for (std::vector<unsigned long long>::const_iterator it = maxSteps.begin(); it != maxSteps.end(); ++ it) {
-    localTrackSteps_m.push(*it);
-  }
-  for (std::vector<double>::const_iterator it = dt.begin(); it != dt.end(); ++ it) {
-    dtAllTracks_m.push(*it);
-  }
-  for (std::vector<double>::const_iterator it = zstop.begin(); it != zstop.end(); ++ it) {
-    zStop_m.push(*it);
-  }
+    if ( zstop.size() > 1 )
+        throw OpalException("ThickTracker::ThickTracker()",
+                            "Multiple tracks not yet supported.");
+    
+    
+    CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
+                                   beamline.getCoordTransformationTo());
+    referenceToLabCSTrafo_m = labToRef.inverted();
 }
 
 
@@ -128,57 +127,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);
-//     }
-// }
-
-
-
-
-void ThickTracker::selectDT() {
-    if (itsBunch_m->getIfBeamEmitting()) {
-        double dt = itsBunch_m->getEmissionDeltaT();
-        itsBunch_m->setdT(dt);
-    } else {
-        double dt = dtCurrentTrack_m;
-        itsBunch_m->setdT(dt);
-    }
-}
-
-void ThickTracker::changeDT() {
-    selectDT();
-    const unsigned int localNum = itsBunch_m->getLocalNum();
-    for (unsigned int i = 0; i < localNum; ++ i) {
-        itsBunch_m->dt[i] = itsBunch_m->getdT();
-    }
-}
-
 /*
 //TODO complete and test fringe fields
 void ThickTracker::insertFringeField(SBend* pSBend, lstruct_t& mBL,
@@ -248,10 +196,7 @@ void ThickTracker::execute() {
     OpalData::getInstance()->setInPrepState(true);
 
     OpalData::getInstance()->setGlobalPhaseShift(0.0);
-    //Beam *beam = Beam::find(Attributes::getString(itsAttr[BEAM]));
-    dtCurrentTrack_m = itsBunch_m->getdT();
-
-
+    
     if ( OpalData::getInstance()->hasPriorTrack() ||
          OpalData::getInstance()->inRestartRun() )
     {
@@ -267,27 +212,13 @@ 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);
 
-
-    //files for analysis
-
 #ifdef PHIL_WRITE
-
-    msg<<"P: " <<itsBunch_m->getP()<< endl;
-    msg<<"beta: " << itsBunch_m->getInitialBeta()<< endl;
-    msg<<"gamma: " << itsBunch_m->getInitialGamma()<< endl;
-    msg<<"E0:  " << itsBunch_m->getM() <<endl;
-
-    std::ofstream tmap;
-    tmap.open ("TransferMap.txt");
-    tmap << std::setprecision(16);
-
-
-    std::ofstream tplot; // stat
-    tplot.open ("tunePlot.txt");
-    tmap << std::setprecision(16);
-
     std::ofstream disp; // stat
     disp.open("dispersion.txt");
     disp << std::setprecision(10);
@@ -298,19 +229,7 @@ void ThickTracker::execute() {
     // TODO need a way to implement Initial dispersion Values
     //dispInitialVal[0][0]=  0.5069938765;
     //dispInitialVal[0][1]= -0.1681363086;
-
-
-    tplot.close();
-    tplot.open("tunePlot.txt", std::ios::app);
-
-
-    setTime();
-    
-    double t = itsBunch_m->getT();
-    itsBunch_m->setT(t);
-    
     OpalData::getInstance()->setInPrepState(false);
-    selectDT();
 
     track_m();
     
@@ -366,6 +285,14 @@ void ThickTracker::fillGaps_m() {
         currentEnd = pos + std::get<2>(*el);
         ++el;
     }
+    
+    double length = std::abs(zstop_m - currentEnd);
+    if ( length > threshold_m ) {
+        //FIXME if gamma changes this is an error!!!
+        double gamma = itsReference.getGamma();
+        tmp.push_back(std::make_tuple(hamiltonian_m.drift(gamma), 1, length));
+    }
+    
     elements_m.swap(tmp);
 }
 
@@ -373,11 +300,7 @@ void ThickTracker::fillGaps_m() {
 void ThickTracker::track_m()
 {
     IpplTimings::startTimer(mapTracking_m);
-#ifdef PHIL_WRITE
-    std::ofstream tplot;
-    tplot.open("tunePlot.txt", std::ios::app);
-#endif
-
+    
     fMatrix_t tFMatrix;
     for (int i=0; i<6; i++){
         tFMatrix[i][i]=1.;
@@ -447,6 +370,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 ) {
@@ -474,45 +421,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] );
 }
 
 
@@ -587,10 +548,9 @@ void ThickTracker::dump_m() {
     
     const std::size_t step = itsBunch_m->getGlobalTrackStep();
     
-    bool psDump   = (step + 1) % Options::psDumpFreq;
-    bool statDump = (step + 1) % Options::statDumpFreq;
+    bool psDump   = ((step + 1) % Options::psDumpFreq == 0);
+    bool statDump = ((step + 1) % Options::statDumpFreq == 0);
     
-
     // Sample fields at (xmin, ymin, zmin), (xmax, ymax, zmax) and the centroid location. We
     // are sampling the electric and magnetic fields at the back, front and
     // center of the beam.
@@ -630,23 +590,22 @@ void ThickTracker::dump_m() {
 void ThickTracker::update_m(const double& spos,
                             const std::size_t& step)
 {
-    itsBunch_m->set_sPos(spos);
-    itsBunch_m->setGlobalTrackStep(step);
-    itsBunch_m->calcBeamParameters();
-    itsBunch_m->calcEMean();
-    
-//     itsBunch_m->setT();
-    
+    // update dt and t
+    double ds = spos - itsBunch_m->get_sPos();    
+    double gamma = itsBunch_m->get_gamma();
+    double beta  = std::sqrt(1.0 - 1.0 / (gamma * gamma) );
+    double dt = ds / Physics::c / beta;
+    itsBunch_m->setdT(dt);
+    itsBunch_m->setT(itsBunch_m->getT() + dt);
     
-//     changeDT();
-
-    
-}
-
-
-void ThickTracker::setTime() {
     const unsigned int localNum = itsBunch_m->getLocalNum();
     for (unsigned int i = 0; i < localNum; ++i) {
-        itsBunch_m->dt[i] = itsBunch_m->getdT();
+        itsBunch_m->dt[i] = dt;
     }
+    
+    
+    itsBunch_m->set_sPos(spos);
+    itsBunch_m->setGlobalTrackStep(step);
+    itsBunch_m->calcBeamParameters();
+    itsBunch_m->calcEMean();    
 }