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,21 @@
 //
 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)
+    , 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();
 }
 
 
@@ -74,26 +74,21 @@ ThickTracker::ThickTracker(const Beamline &beamline,
     , 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();
 }
 
 
@@ -161,24 +156,6 @@ void ThickTracker::prepareSections() {
 
 
 
-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 +225,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() )
     {
@@ -269,16 +243,7 @@ void ThickTracker::execute() {
     
     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 disp; // stat
     disp.open("dispersion.txt");
     disp << std::setprecision(10);
@@ -289,15 +254,7 @@ void ThickTracker::execute() {
     // TODO need a way to implement Initial dispersion Values
     //dispInitialVal[0][0]=  0.5069938765;
     //dispInitialVal[0][1]= -0.1681363086;
-
-
-    setTime();
-    
-    double t = itsBunch_m->getT();
-    itsBunch_m->setT(t);
-    
     OpalData::getInstance()->setInPrepState(false);
-    selectDT();
 
     track_m();
     
@@ -353,6 +310,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);
 }
 
@@ -612,24 +577,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();
-    
-    //FIXME update time
-//     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();    
 }