upadting

src/Distribution/Distribution.cpp

@@ -1362,19 +1362,18 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
           Units of the Sigma Matrix:
 
           spatial: mm
-          moment:  rad
+          moment:  mrad
 
         */
 
         if (Options::cloTuneOnly)
             throw OpalException("Do only CLO and tune calculation","... ");
 
-
         auto sigma = siggen->getSigma();
-
+        // convertion from mm mrad to m rad 
         for (unsigned int i = 0; i < 6; ++ i)
             for (unsigned int j = 0; j < 6; ++ j) sigma(i, j) *= 1e-6;
-        
+                
         for (unsigned int i = 0; i < 3; ++ i) {
             if ( sigma(2 * i, 2 * i) < 0 || sigma(2 * i + 1, 2 * i + 1) < 0 )
                 throw OpalException("Distribution::CreateMatchedGaussDistribution()",
@@ -1387,8 +1386,11 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
         sigmaP_m[2] = std::sqrt(sigma(3, 3))*beta*gamma;
         sigmaR_m[1] = std::sqrt(sigma(4, 4));
         // delta = |p - p_0|/|p_0| = |p - p_0|/beta*gamma --> l'^2 = (delta*beta*gamma)^2 - (x'^2 + z'^2) 
-        double lprime2 = sigma(5, 5)*beta*gamma - sigmaP_m[0]*sigmaP_m[0] - sigmaP_m[2]*sigmaP_m[2];
-        sigmaP_m[1] = std::sqrt(std::abs(lprime2));
+        //double lprime2 = sigma(5, 5)*(beta*beta*gamma*gamma) - sigmaP_m[0]*sigmaP_m[0] - sigmaP_m[2]*sigmaP_m[2];
+        //sigmaP_m[1] = std::sqrt(std::abs(lprime2));
+        double pl2 = (std::sqrt(sigma(5,5)) + 1)*(std::sqrt(sigma(5,5)) + 1)*beta*gamma*beta*gamma - sigmaP_m[0]*sigmaP_m[0] - sigmaP_m[2]*sigmaP_m[2];
+        double pl = std::sqrt(pl2);
+        sigmaP_m[1] = gamma*(pl - beta*gamma);
 
         correlationMatrix_m(1, 0) = sigma(0, 1) / (sqrt(sigma(0, 0) * sigma(1, 1)));
         correlationMatrix_m(3, 2) = sigma(2, 3) / (sqrt(sigma(2, 2) * sigma(3, 3)));
@@ -2502,7 +2504,7 @@ void Distribution::generateGaussZ(size_t numberOfParticles) {
     }
     //Sets again the standard GSL error handler on
     gsl_set_error_handler(NULL);
-
+    //Just to be sure
     if (errcode == GSL_EDOM) {
         throw OpalException("Distribution::GenerateGaussZ",
                             "gsl_linalg_cholesky_decomp failed");
@@ -3490,11 +3492,11 @@ void Distribution::setAttributes() {
     // itsAttr[Attrib::Distribution::E2]
     //     = Attributes::makeReal("E2", "If E2<Eb, we compute the tunes from the beams energy Eb to E2 with dE=0.25 MeV ", 0.0);
     itsAttr[Attrib::Distribution::RESIDUUM]
-        = Attributes::makeReal("RESIDUUM", "Residuum for the closed orbit finder and sigma matrix generator ", 1e-8);
+        = Attributes::makeReal("RESIDUUM", "Residuum for the closed orbit finder and sigma matrix generator ", 1e-9);
     itsAttr[Attrib::Distribution::MAXSTEPSCO]
         = Attributes::makeReal("MAXSTEPSCO", "Maximum steps used to find closed orbit ", 100);
     itsAttr[Attrib::Distribution::MAXSTEPSSI]
-        = Attributes::makeReal("MAXSTEPSSI", "Maximum steps used to find matched distribution ",2000);
+        = Attributes::makeReal("MAXSTEPSSI", "Maximum steps used to find matched distribution ",1000);
     itsAttr[Attrib::Distribution::ORDERMAPS]
         = Attributes::makeReal("ORDERMAPS", "Order used in the field expansion ", 7);
     itsAttr[Attrib::Distribution::MAGSYM]

src/Distribution/MapGenerator.h

@@ -147,7 +147,7 @@ MapGenerator<Value_type,
     
     for (size_type i = 0; i < 6; ++i) {
         for (size_type j = 0; j < 6; ++j) {
-            map(i,j) = matrix[i][j];
+	    map(i,j) = matrix[i][j];
         }
     }
     

src/Distribution/SigmaGenerator.h

@@ -316,6 +316,11 @@ class SigmaGenerator
          * Mcycs is a vector of all cyclotron maps
          */
         void updateSigma(const std::vector<matrix_type>&, const std::vector<matrix_type>&);
+	/// Writes the one turn matrix 
+	/*
+	 * @param One turn Matrix and iteration number
+	 */
+	void writeOneTurn(const matrix_type&, const value_type);
 	
 	/// Writes the Matched Distributions
 	/*
@@ -470,9 +475,9 @@ template<typename Value_type, typename Size_type>
      * - step sizes of path ds
      * - tune nuz
      */
-    // object for space cyclotron map
-    maxit = 150;
+    
     try{
+        // object for space cyclotron map
         MapGenerator<value_type,size_type,Series,Map,Hamiltonian,SpaceCharge> mapgen(nStepsPerSector_m);
 
 	container_type h(nStepsPerSector_m), r(nStepsPerSector_m), ds(nStepsPerSector_m), fidx(nStepsPerSector_m);
@@ -498,7 +503,7 @@ template<typename Value_type, typename Size_type>
 	    tunes = cof.getTunes();
 	    ravg = cof.getAverageRadius();                   // average radius
 
-	    container_type peo = cof.getMomentum(angle);
+	    container_type peo = cof.getMomentum(angle);    //radial Momentum
 
 	    // write to terminal
 	    *gmsg << "* ----------------------------" << endl
@@ -518,12 +523,12 @@ template<typename Value_type, typename Size_type>
 	    size_type start = deg_step * angle;
 
 	    // copy properties of the length of one sector (--> nStepsPerSector_m)
-	    nStepsPerSector_m = ds_turn.size();
+	    
 	    std::copy_n(r_turn.begin()+start,nStepsPerSector_m, r.begin());
 	    std::copy_n(h_turn.begin()+start,nStepsPerSector_m, h.begin());
 	    std::copy_n(fidx_turn.begin()+start,nStepsPerSector_m, fidx.begin());
 	    std::copy_n(ds_turn.begin()+start,nStepsPerSector_m, ds.begin());
-
+	    
 	    if(write_m){
 	        writeOrbitOutput_m(tunes,ravg, cof.getFrequencyError(),r_turn,peo,h_turn,fidx_turn,ds_turn);
 	    }
@@ -533,23 +538,27 @@ template<typename Value_type, typename Size_type>
 	    return false;
 	}
 
-	if(Options::cloTuneOnly)
-	  throw OpalException("Do only CLO and tune calculation","... ");
-
-	// computes the cyclotron maps for each angle and stores them into a vector
-	std::vector<matrix_type> Mcycs(ds.size());
-    
-	// calculate only for a single sector (a nSector_-th) of the whole cyclotron
-	for (size_type i = 0; i < ds.size(); ++i) 
-	    Mcycs[i] = mapgen.generateMap(H_m(h[i],h[i]*h[i]+fidx[i],-fidx[i]),ds[i],truncOrder_m);
+	if(Options::cloTuneOnly){
+	  //throw OpalException("Do only CLO and tune calculation","... ");
+	  *gmsg << "Doing just close orbit and tune calculation only" << endl;
+	}
+	else{
+	    // computes the cyclotron maps for each angle and stores them into a vector
+	    std::vector<matrix_type> Mcycs(nStepsPerSector_m);
+	  
+	    // calculate only for a single sector (a nSector_-th) of the whole cyclotron
+	    for (size_type i = 0; i < nStepsPerSector_m; ++i) 
+	        Mcycs[i] = mapgen.generateMap(H_m(h[i],h[i]*h[i]+fidx[i],-fidx[i]),ds[i],truncOrder_m);
     
-	// search for a match distribution for each permutation 
-	for(size_type permutation = 0; permutation < 1; permutation++){
-	    findMatchDistribution(accuracy,h,fidx,ds,tunes,ravg,maxit,Mcycs,harmonic,2);
+	    // search for a match distribution for each permutation
+	    // Set to search for the permutation which converged to a matched distribution upon tracking
+	    // Note that 4 permutations converge to the same matched distribution
+	    for(size_type permutation = 0; permutation < 1; permutation++){
+	        findMatchDistribution(accuracy,h,fidx,ds,tunes,ravg,maxit,Mcycs,harmonic,2);
+	    }
+	    std::cout << "OPAL> * Number of matched distributions: " <<matchedSigmas_m.size() << std::endl;
+	    std::cout << "OPAL> * Last matched distribution information: "<<std::endl;
 	}
-	std::cout << "OPAL> * Number of matched distributions: " <<matchedSigmas_m.size() << std::endl;
-	std::cout << "OPAL> * Last matched distribution information: "<<std::endl;
-	
     }catch(const std::exception& e) {
         std::cerr << e.what() << std::endl;
     }
@@ -557,7 +566,7 @@ template<typename Value_type, typename Size_type>
     if ( 0 < matchedSigmas_m.size()) return (0 < matchedSigmas_m.size());
     //Returns an identity matrix, if  no matched distribution was found
     else{
-        *gmsg << "No matched distribution found, returning unity matrix instead.";
+        *gmsg << "No matched distribution found, returning unity matrix instead or CLOTUNEONLY option"<<endl;
         matrix_type One(6,6);
 	for(size_type i = 0; i < 6; i++){
 	  for(size_type j = 0; j < 6; j++) One(i,j) = 0;
@@ -778,7 +787,7 @@ template<typename Value_type, typename Size_type>
 
     try {
         // stores computed space charge map for each angle
-        std::vector<matrix_type> Mscs(ds.size());
+        std::vector<matrix_type> Mscs(nStepsPerSector_m);
 
 	// (inverse) transformation matrix
         complex_matrix_type R(6,6), invR(6,6);
@@ -789,13 +798,14 @@ template<typename Value_type, typename Size_type>
 
 	// for exiting the loop if maxit is reached
         bool stop = false;
+	// checks if the eigenvector matrix E and E^-1 could be built
+	bool isBuild = false;
 
 	// initialize sigma matrices (for each angle one) (first guess)
         initialize(tunes.second,ravg);
 
         // calculate only for a single sector (a nSector_-th) of the whole cyclotron
-	// for (size_type i = 0; i < nStepsPerSector_m; ++i) {
-        for (size_type i = 0; i < ds.size(); ++i) {
+        for (size_type i = 0; i < nStepsPerSector_m; ++i) {
             if (!harmonic) {
                 Mscs[i]  = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),sigmas_m[i](2,2),sigmas_m[i](4,4)),ds[i],truncOrder_m);
              } else {
@@ -807,16 +817,19 @@ template<typename Value_type, typename Size_type>
          mapgen.combine(Mscs,Mcycs);
          matrix_type Mturn = mapgen.getMap();
 
+	 writeOneTurn(Mturn, niterations_m);
+
          while (error_m > accuracy && !stop) {
 	     /*
 	      * Version 2.0 decouple was substituded by setEigenVectors
 	      * and updateInitialSigma was changed according to Eq. 18 Sec. 2.4 
 	      * Frey Semester thesis
 	      */
-	     // Returns the Eigenvalues of the Twiss-Matrix,
+	     // Sets the Eigenvalues of the Twiss-Matrix,
 	     // orders them and computes R and invR
-	     setEigenVectors(Mturn,R,invR,permutation);
-	 	
+	     // Also checks if the Eigenvector matrix E is invertible
+	     isBuild = setEigenVectors(Mturn,R,invR,permutation);
+	     
 	     // computes the new sigma matrix
 	     newSigma = updateInitialSigma(R,invR);
 		
@@ -830,7 +843,7 @@ template<typename Value_type, typename Size_type>
 	     sigmas_m[0] = weight*newSigma + (1.0-weight)*sigmas_m[0];
 
 	     // compute new space charge maps
-	     for (size_type i = 0; i < ds.size(); ++i) {
+	     for (size_type i = 0; i < nStepsPerSector_m; ++i) {
 	         if (!harmonic) {
 		     Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),sigmas_m[i](2,2),sigmas_m[i](4,4)),ds[i],truncOrder_m);
 		 } else {
@@ -841,9 +854,11 @@ template<typename Value_type, typename Size_type>
 	     // construct new one turn transfer matrix M
 	     mapgen.combine(Mscs,Mcycs);
 	     Mturn = mapgen.getMap();
-
-	     // check if number of iterations has maxit exceeded.
-	     stop = (niterations_m++ > maxit);
+	     
+	     writeOneTurn(Mturn, niterations_m);
+	     
+	     // check if number of iterations has maxit exceeded or if the eigenvector matrix could be built
+	     stop = (niterations_m++ > maxit) || !isBuild;
 
 	     // check for convergence
 	     converged_m = error_m < accuracy;
@@ -916,20 +931,17 @@ template<typename Value_type, typename Size_type>
     gsl_matrix_complex_free(eigenVectors);
     gsl_matrix_free(MTurngsl);
 
-    // Sorts the Eigenvectors such that
-    // Re(eig_x) < Re(eig_y) < Re(eig_z)
-
+    // Sorts the Eigenvectors
     std::vector<size_type> permutation = permutations_m[perm];
-
+    value_type threshold = 10e-12;
+    std::vector<complex_vector_type> zVectors{};
+    std::vector<complex_vector_type> xyVectors{};
     bool isZdirection = false;
-    std::vector<complex_vector_type> zVectors;
-    std::vector<complex_vector_type> xyVectors; 
-
+    
     for(size_type i = 0; i < 6; i++){
         complex_number_type z = R(i,0);
-	if(std::abs(z) < 1e-13) z = 0.;
+	if(std::abs(z) < threshold) z = 0.;
 	if(z == 0.) isZdirection = true;
-	
 	complex_vector_type v(6);
 	if(isZdirection){
 	    for(size_type j = 0;j < matrixDimension; j++){
@@ -945,7 +957,17 @@ template<typename Value_type, typename Size_type>
 	    }
 	    xyVectors.push_back(v);
 	}
-        isZdirection = false;
+	isZdirection = false;
+    }
+    //z-direction not found, then the system does not have the expected form
+    if(zVectors.size() != 2){
+        *gmsg << "Could not found the vertical eigenvectors of the system"<<endl;
+        complex_matrix_type Copy = R;
+        for(size_type i = 0; i < 6; i++){
+	    for(size_type j = 0; j < 6; j++) R(i,j) = Copy(j,i);
+        }
+        InvertMatrix(R,invR);
+        return false;
     }
     // Norms the Eigenvectors
      for(size_type i = 0; i < 4; i++){
@@ -1057,8 +1079,28 @@ typename SigmaGenerator<Value_type, Size_type>::value_type SigmaGenerator<Value_
     return std::sqrt(std::inner_product(diff.data().begin(),diff.data().end(),diff.data().begin(),0.0));
 }
 
+template<typename Value_type, typename Size_type>  void SigmaGenerator<Value_type, Size_type>::writeOneTurn(const matrix_type& M, const value_type iteration){
+
+    std::stringstream ss;
+    ss << iteration;
+    std::string it = ss.str();
+
+    std::string energy = float2string(E_m);
+    std::ofstream oneTurn("data/oneTurn"+energy+"MeV.txt", std::ios::app);
+    oneTurn<<"Iteration: "+it<<std::endl;
+    for(size_type i = 0; i < 6; i++){
+        for(size_type j = 0; j < 6; j++){
+	    oneTurn<<std::left<<std::setprecision(5)
+		   <<std::setw(12)<<M(i,j);
+        }
+        oneTurn<<std::endl;
+    }
+    oneTurn<<std::endl;
+    oneTurn.close();
+}
+
 template<typename Value_type, typename Size_type> void SigmaGenerator<Value_type, Size_type>::writeMatched(const matrix_type& match, const size_type permutation,const std::vector<matrix_type>& Mcycs, const std::vector<matrix_type>& Mscs, const matrix_type& Mcyc){
-    
+
     std::string energy = float2string(E_m);
     std::ofstream writeSigmas("data/Matches"+energy+"MeV.txt", std::ios::app);
     writeSigmas<<std::left<<std::setprecision(5)