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

@@ -196,7 +196,7 @@ ThickTracker::map_t ThickTracker::getSpaceChargeMap(double& length){
 /**Hamiltonian for a linear Thin Lens fringe field approximation
  * \f[ H_{ThinLens} = \frac{1}{2} (x^2 - y^2) k_0 \tan \left( \Psi \right) \f]
  */
-void ThickTracker::insertThinLensApproxFringe(std::list<structMapTracking>& mBL, double phi, double& K0, double& pos){
+void ThickTracker::insertThinLensApproxFringe(lstruct_t& mBL, double phi, double& K0, double& pos){
 
     structMapTracking ThinLensFringeApprox;
     series_t H;
@@ -220,13 +220,13 @@ void ThickTracker::insertThinLensApproxFringe(std::list<structMapTracking>& mBL,
 
 
 //TODO complete and test fringe fields
-void ThickTracker::insertFringeField(SBend* pSBend, std::list<structMapTracking>& mBL,
+void ThickTracker::insertFringeField(SBend* pSBend, lstruct_t& mBL,
         double& beta0, double& gamma0, double& P0, double& q, std::array<double,2>& paramFringe, std::string e){
     
 
     series_t H;
     map_t tempMap, fieldMap;
-    std::list<structMapTracking>::iterator mBLit;
+    lstruct_t::iterator mBLit;
     double lenFringe = paramFringe[0] + paramFringe[1];
 
     //double intHeight = 0.1; //:TODO: change or argument for that value
@@ -277,7 +277,7 @@ void ThickTracker::insertFringeField(SBend* pSBend, std::list<structMapTracking>
 /** \f[H_{Drift}= \frac{\delta}{\beta_0} -
 * \sqrt{\left(\frac{1}{\beta_0} + \delta \right)^2 -p_x^2 -p_y^2 - \frac{1}{\left(\beta_0 \gamma_0\right)^2 } } \f]
 */
-void ThickTracker::fillDrift( std::list<structMapTracking>& mapBeamLine,double& elementPos,
+void ThickTracker::fillDrift( lstruct_t& mapBeamLine,double& elementPos,
         double& undefSpace){
 
     Inform msg("ThickTracker", *gmsg);
@@ -341,7 +341,7 @@ void ThickTracker::fillDrift( std::list<structMapTracking>& mapBeamLine,double&
 */
 
 void ThickTracker::defMapTrackingElement(std::shared_ptr<Component> element, structMapTracking& elSrct,
-        std::list<structMapTracking>& mBL){
+        lstruct_t& mBL){
 
 #ifdef PHIL_WRITE
     std::ofstream mbl;
@@ -440,7 +440,7 @@ void ThickTracker::defMapTrackingElement(std::shared_ptr<Component> element, str
                     std::cout << "we are here 1"<< std::endl;
                     if (entrFringe[0]>0 && !mBL.empty()){
                         std::cout << "we are here 2"<< std::endl;
-                        std::list<structMapTracking>::iterator mBLit;
+                        lstruct_t::iterator mBLit;
                         //redefine previous element
                         mBLit = mBL.end();
                         mBLit--;
@@ -563,13 +563,6 @@ void ThickTracker::execute() {
     prepareSections();
     
     
-    msg << std::setprecision(10);
-    msg << itsReference.getE() << endl;
-    msg << itsReference.getGamma() * itsReference.getBeta() << endl;
-    msg << itsBunch_m->getInitialBeta() * itsBunch_m->getInitialGamma() << endl;
-    msg << *itsBunch_m << endl;
-    msg << std::setprecision(10);
-
     bool consSC = itsBunch_m->hasFieldSolver();
     std::string fieldSolver = itsBunch_m->getFieldSolverType();
 
@@ -592,14 +585,11 @@ void ThickTracker::execute() {
 
     map_t combinedMaps; //Final Transfer map_t
     series_t H;         //createHamiltonian
-    OpalParticle part;
 
     structMapTracking mapTrackingElement;
 
-    std::list<structMapTracking> mapBeamLine;
-    std::list<structMapTracking>::iterator mapBeamLineit;
-
-    std::vector<map_t> mapVec; // not necessary
+    lstruct_t mapBeamLine;
+    lstruct_t::iterator mapBeamLineit;
 
     double positionMapTracking = zstart_m;
 
@@ -726,69 +716,24 @@ void ThickTracker::execute() {
     //=================================
     IpplTimings::stopTimer(mapCreation_m);
 
-    IpplTimings::startTimer(mapCombination_m);
-    // combine Maps
-    std::size_t totalSlices=0;
-
-    fMatrix_t tFMatrix, sFMatrix;
-    double position=zstart_m;
-
-    for (int i=0; i<6; i++){
-        tFMatrix[i][i]=1.;
-    }
-
-    FMatrix<double, 1, 4> dispInitialVal;
-    for (int i=2; i<4; i++){
-        dispInitialVal[0][i]=0;
-    }
-
+    
     // TODO need a way to implement Initial dispersion Values
     //dispInitialVal[0][0]=  0.5069938765;
     //dispInitialVal[0][1]= -0.1681363086;
 
 
-    trackDisp(tFMatrix, dispInitialVal, position);
-
     tplot.close();
     tplot.open("tunePlot.txt", std::ios::app);
 
 
-    //combined map_t
-    for (mapBeamLineit=mapBeamLine.begin(); mapBeamLineit != mapBeamLine.end(); ++mapBeamLineit) {
-
-        for (std::size_t slice=0; slice < mapBeamLineit->nSlices; slice++){
-
-            combinedMaps = mapBeamLineit->elementMap * combinedMaps;
-            combinedMaps = combinedMaps.truncate(truncOrder_m);
-            tFMatrix= combinedMaps.linearTerms();
-            sFMatrix= mapBeamLineit->elementMap.linearTerms();
-
-
-            if (!considerSC_m) {
-                tplot << "sPos=  " << position + mapBeamLineit->stepSize << std::endl;
-                //linTAnalyze(tFMatrix);
-            }
-
-            position+= mapBeamLineit->stepSize;
-            totalSlices++;
-
-
-#ifdef PHIL_WRITE
-            trackDisp(tFMatrix, dispInitialVal, position);
-#endif
-        }
-    }
-    IpplTimings::stopTimer(mapCombination_m);
-
-
-#ifdef PHIL_WRITE
-
-    //eliminate higher order terms (occurring through multiplication)
-
-	tmap << "Transfermap" << std::endl;
-    tmap << "Gantry2" << std::endl;
-    tmap << combinedMaps << std::endl;
-#endif
+// #ifdef PHIL_WRITE
+// 
+//     //eliminate higher order terms (occurring through multiplication)
+// 
+// 	tmap << "Transfermap" << std::endl;
+//     tmap << "Gantry2" << std::endl;
+//     tmap << combinedMaps << std::endl;
+// #endif
 
 
 
@@ -800,8 +745,7 @@ void ThickTracker::execute() {
 	OpalData::getInstance()->setInPrepState(false);
 	selectDT();
 
-    trackParticles_m(
-        mapBeamLine);
+    track_m(mapBeamLine);
     
 //-------------------------
 
@@ -811,9 +755,9 @@ void ThickTracker::execute() {
 }
 
 
-void ThickTracker::trackParticles_m(
-    const std::list<structMapTracking>& mapBeamLine) {
-IpplTimings::startTimer(mapTracking_m);
+void ThickTracker::track_m(const lstruct_t& mapBeamLine)
+{
+    IpplTimings::startTimer(mapTracking_m);
 #ifdef PHIL_WRITE
     std::ofstream tplot;
     tplot.open("tunePlot.txt", std::ios::app);
@@ -823,76 +767,59 @@ IpplTimings::startTimer(mapTracking_m);
     for (int i=0; i<6; i++){
         tFMatrix[i][i]=1.;
     }
+    
+    FMatrix<double, 1, 4> dispInitialVal;
+    for (int i=2; i<4; i++){
+        dispInitialVal[0][i]=0;
+    }
+    
+    map_t transferMap;
 
     int sliceidx=0;
-    cfMatrix_t N, invN;
-    FVector<double, 6> particle;
-    map_t scMap;
+//     cfMatrix_t N, invN;
+//     map_t scMap;
     dumpStats(sliceidx, true, true);
 
     fMatrix_t refSigma;
-    double betagamma=itsBunch_m->getInitialBeta() * itsBunch_m->getInitialGamma();
     refSigma = (itsBunch_m->getSigmaMatrix());
 //     setNMatrix(refSigma, N, invN);
 
+    double position = zstart_m;
+    
+    
+    this->advanceDispersion_m(tFMatrix, dispInitialVal, position);
+    
     //(1) Loop Beamline
     for(auto mapBeamLineit=mapBeamLine.begin(); mapBeamLineit != mapBeamLine.end(); ++mapBeamLineit) {
-
         //(2) Loop Slices
-        for (std::size_t slice=0; slice < mapBeamLineit->nSlices; slice++){
+        for (std::size_t slice=0; slice < mapBeamLineit->nSlices; slice++) {
 
             //apply SC on map
-            if (considerSC_m){
-                double stepSC = mapBeamLineit->stepSize;
-                scMap = getSpaceChargeMap(stepSC);
-
-                tplot << "sPos=  " << mapBeamLineit->elementPos + mapBeamLineit->stepSize * (slice+1) << std::endl;
-                tFMatrix = scMap.linearTerms() * mapBeamLineit->elementMap.linearTerms() * tFMatrix;
-                mapAnalyser_m.linTAnalyze(tFMatrix);
-            }
-
-            //(3) Loop Particles
-            for (unsigned int partidx=0; partidx< itsBunch_m->getLocalNum(); ++partidx){
-
-
-
-                for (int d = 0; d < 3; ++d) {
-                    particle[2 * d] = itsBunch_m->R[partidx](d);
-                    particle[2 *d + 1] = itsBunch_m->P[partidx](d);
-                }
-
-                //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+pParticle*pParticle)/betagamma)
-                                -1./itsBunch_m->getInitialBeta();
-
-                //Apply Space charges
-                if (considerSC_m){
-                    particle = scMap * particle;
-                }
-
-                //Apply element map
-                particle = mapBeamLineit->elementMap * particle;
-
-                //Units back
-                particle[1]*=betagamma;
-                particle[3]*=betagamma;
-
-                pParticle =  std::sqrt( ((particle[5] + 1./itsBunch_m->getInitialBeta()) * betagamma)
-                             *((particle[5] + 1./itsBunch_m->getInitialBeta()) * betagamma) -1);
-
-                particle[5] = std::sqrt(pParticle*pParticle - particle[1]*particle[1] - particle[3]*particle[3]);
-
-                itsBunch_m->set_part(particle, partidx);
-
-
-
-            }
+//             if (considerSC_m){
+//                 double stepSC = mapBeamLineit->stepSize;
+// //                 scMap = getSpaceChargeMap(stepSC);
+// 
+//                 tplot << "sPos=  " << mapBeamLineit->elementPos + mapBeamLineit->stepSize * (slice+1) << std::endl;
+//                 tFMatrix = scMap.linearTerms() * mapBeamLineit->elementMap.linearTerms() * tFMatrix;
+//                 mapAnalyser_m.linTAnalyze(tFMatrix);
+//             }
+
+            this->advanceParticles_m(mapBeamLineit);
+            
+            this->concatenateMaps_m(mapBeamLineit->elementMap, transferMap);
+            
+            tFMatrix= transferMap.linearTerms();
+
+//             if (!considerSC_m) {
+//                 tplot << "sPos=  " << position + mapBeamLineit->stepSize << std::endl;
+//                 //linTAnalyze(tFMatrix);
+//             }
+            
+            position += mapBeamLineit->stepSize;
+            
+            this->advanceDispersion_m(tFMatrix, dispInitialVal, position);
+            
+            
             bool const psDump = true;   //((itsBunch_m->getGlobalTrackStep() % Options::psDumpFreq) + 1 == Options::psDumpFreq);
             bool const statDump = true; //((itsBunch_m->getGlobalTrackStep() % Options::statDumpFreq) + 1 == Options::statDumpFreq);
 
@@ -902,11 +829,65 @@ IpplTimings::startTimer(mapTracking_m);
             refSigma = itsBunch_m->getSigmaMatrix();
 
             //printPhaseShift(refSigma ,mapBeamLineit->elementMap.linearTerms(), N);
-            sliceidx ++;
+            sliceidx++;
         }
     }
     IpplTimings::stopTimer(mapTracking_m);
+}
+
+
+void ThickTracker::advanceParticles_m(lstruct_t::const_iterator& mit) {
+    
+    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);
+        }
+
+        //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();
+        
+//         //Apply Space charges
+//         if (considerSC_m){
+//             particle = scMap * particle;
+//         }
+
+        //Apply element map
+        particle = mit->elementMap * particle;
+
+        //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);
+
+        particle[5] = std::sqrt(pParticle * pParticle -
+                                particle[1] * particle[1] -
+                                particle[3] * particle[3] );
+        
+        itsBunch_m->set_part(particle, ip);
+    }
+}
+
+
+void ThickTracker::concatenateMaps_m(const map_t& x, map_t& y) {
+    IpplTimings::startTimer(mapCombination_m);
+    y = x * y;
+    IpplTimings::stopTimer(mapCombination_m);
 }
 
 
@@ -919,7 +900,10 @@ IpplTimings::startTimer(mapTracking_m);
  *  +
  *  \begin{pmatrix} R_{16} \\ R_{26} \end{pmatrix}\f]
  */
-void ThickTracker::trackDisp(fMatrix_t tempMatrix, FMatrix<double, 1, 4> initialVal, double pos) {
+void ThickTracker::advanceDispersion_m(fMatrix_t tempMatrix,
+                                       FMatrix<double, 1, 4> initialVal,
+                                       double pos)
+{
 
 #ifdef PHIL_WRITE
     std::ofstream disp;