diff --git a/src/Distribution/CMakeLists.txt b/src/Distribution/CMakeLists.txt
index 9f6797cf9171a261fd0be89694a3036cdd43a31d..74eadb79ff74f20ed5fbda63318355db0b450050 100644
--- a/src/Distribution/CMakeLists.txt
+++ b/src/Distribution/CMakeLists.txt
@@ -16,6 +16,7 @@ set (HDRS
LaserProfile.h
MapGenerator.h
matrix_vector_operation.h
+ rdm.h
SigmaGenerator.h
)
diff --git a/src/Distribution/SigmaGenerator.cpp b/src/Distribution/SigmaGenerator.cpp
index 372546937a26cda5137c59ea699796c6f987e206..4a7e551b598d246bd60c88f388ff942d06f45d6d 100644
--- a/src/Distribution/SigmaGenerator.cpp
+++ b/src/Distribution/SigmaGenerator.cpp
@@ -323,10 +323,10 @@ bool SigmaGenerator::match(double accuracy,
while (error_m > accuracy && !stop) {
// decouple transfer matrix and compute (inverse) tranformation matrix
- decouple(Mturn,R,invR);
+ vector_type eigen = decouple(Mturn, R,invR);
// construct new initial sigma-matrix
- newSigma = updateInitialSigma(Mturn, R, invR);
+ newSigma = updateInitialSigma(Mturn, eigen, R, invR);
// compute new sigma matrices for all angles (except for initial sigma)
updateSigma(Mscs,Mcycs);
@@ -334,6 +334,15 @@ bool SigmaGenerator::match(double accuracy,
// compute error
error_m = L2ErrorNorm(sigmas_m[0],newSigma);
+ //***
+ std::cout << "Error: " << error_m << std::endl;
+ std::cout << std::abs(sigmas_m[0](0, 0) - sigmas_m[nSteps_m - 1](0, 0)) << std::endl;
+ std::cout << "L1: " << L1ErrorNorm(sigmas_m[0], sigmas_m[nSteps_m - 1]) << std::endl;
+ std::cout << "L2: " << L2ErrorNorm(sigmas_m[0], sigmas_m[nSteps_m - 1]) << std::endl;
+ std::cout << "LInf: " << LInfErrorNorm(sigmas_m[0], sigmas_m[nSteps_m - 1]) << std::endl;
+ std::cin.get();
+ //***
+
// write new initial sigma-matrix into vector
sigmas_m[0] = weight*newSigma + (1.0-weight)*sigmas_m[0];
@@ -427,148 +436,181 @@ bool SigmaGenerator::match(double accuracy,
}
-void SigmaGenerator::eigsolve_m(const matrix_type& Mturn,
- sparse_matrix_type& R)
-{
- typedef gsl_matrix* gsl_matrix_t;
- typedef gsl_vector_complex* gsl_cvector_t;
- typedef gsl_matrix_complex* gsl_cmatrix_t;
- typedef gsl_eigen_nonsymmv_workspace* gsl_wspace_t;
- typedef boost::numeric::ublas::vector<complex_t> complex_vector_type;
-
- gsl_cvector_t evals = gsl_vector_complex_alloc(6);
- gsl_cmatrix_t evecs = gsl_matrix_complex_alloc(6, 6);
- gsl_wspace_t wspace = gsl_eigen_nonsymmv_alloc(6);
- gsl_matrix_t M = gsl_matrix_alloc(6, 6);
-
- // go to GSL
- for (unsigned int i = 0; i < 6; ++i){
- for (unsigned int j = 0; j < 6; ++j) {
- gsl_matrix_set(M, i, j, Mturn(i,j));
- }
- }
-
- /*int status = */gsl_eigen_nonsymmv(M, evals, evecs, wspace);
-
-// if ( !status )
-// throw OpalException("SigmaGenerator::eigsolve_m()",
-// "Couldn't perform eigendecomposition!");
-
- /*status = *///gsl_eigen_nonsymmv_sort(evals, evecs, GSL_EIGEN_SORT_ABS_ASC);
-
-// if ( !status )
+// void SigmaGenerator::eigsolve_m(const matrix_type& Mturn,
+// sparse_matrix_type& R)
+// {
+// typedef gsl_matrix* gsl_matrix_t;
+// typedef gsl_vector_complex* gsl_cvector_t;
+// typedef gsl_matrix_complex* gsl_cmatrix_t;
+// typedef gsl_eigen_nonsymmv_workspace* gsl_wspace_t;
+// typedef boost::numeric::ublas::vector<complex_t> complex_vector_type;
+//
+// gsl_cvector_t evals = gsl_vector_complex_alloc(6);
+// gsl_cmatrix_t evecs = gsl_matrix_complex_alloc(6, 6);
+// gsl_wspace_t wspace = gsl_eigen_nonsymmv_alloc(6);
+// gsl_matrix_t M = gsl_matrix_alloc(6, 6);
+//
+// // go to GSL
+// for (unsigned int i = 0; i < 6; ++i){
+// for (unsigned int j = 0; j < 6; ++j) {
+// gsl_matrix_set(M, i, j, Mturn(i,j));
+// }
+// }
+//
+// /*int status = */gsl_eigen_nonsymmv(M, evals, evecs, wspace);
+//
+// // if ( !status )
+// // throw OpalException("SigmaGenerator::eigsolve_m()",
+// // "Couldn't perform eigendecomposition!");
+//
+// /*status = *///gsl_eigen_nonsymmv_sort(evals, evecs, GSL_EIGEN_SORT_ABS_ASC);
+//
+// // if ( !status )
+// // throw OpalException("SigmaGenerator::eigsolve_m()",
+// // "Couldn't sort eigenvalues and eigenvectors!");
+//
+// // go to UBLAS
+// for( unsigned int i = 0; i < 6; i++){
+// gsl_vector_complex_view evec_i = gsl_matrix_complex_column(evecs, i);
+//
+// for(unsigned int j = 0;j < 6; j++){
+// gsl_complex zgsl = gsl_vector_complex_get(&evec_i.vector, j);
+// complex_t z(GSL_REAL(zgsl), GSL_IMAG(zgsl));
+// R(i,j) = z;
+// }
+// }
+//
+// // Sorting the Eigenvectors
+// // This is an arbitrary threshold that has worked for me. (We should fix this)
+// double threshold = 10e-12;
+// bool isZdirection = false;
+// std::vector<complex_vector_type> zVectors{};
+// std::vector<complex_vector_type> xyVectors{};
+//
+// for(unsigned int i = 0; i < 6; i++){
+// complex_t z = R(i,0);
+// if(std::abs(z) < threshold) z = 0.;
+// if(z == 0.) isZdirection = true;
+// complex_vector_type v(6);
+// if(isZdirection){
+// for(unsigned int j = 0;j < 6; j++){
+// complex_t z = R(i,j);
+// v(j) = z;
+// }
+// zVectors.push_back(v);
+// }
+// else{
+// for(unsigned int j = 0; j < 6; j++){
+// complex_t z = R(i,j);
+// v(j) = z;
+// }
+// xyVectors.push_back(v);
+// }
+// isZdirection = false;
+// }
+//
+// //if z-direction not found, then the system does not behave as expected
+// if(zVectors.size() != 2)
// throw OpalException("SigmaGenerator::eigsolve_m()",
-// "Couldn't sort eigenvalues and eigenvectors!");
-
- // go to UBLAS
- for( unsigned int i = 0; i < 6; i++){
- gsl_vector_complex_view evec_i = gsl_matrix_complex_column(evecs, i);
-
- for(unsigned int j = 0;j < 6; j++){
- gsl_complex zgsl = gsl_vector_complex_get(&evec_i.vector, j);
- complex_t z(GSL_REAL(zgsl), GSL_IMAG(zgsl));
- R(i,j) = z;
- }
- }
-
- // Sorting the Eigenvectors
- // This is an arbitrary threshold that has worked for me. (We should fix this)
- double threshold = 10e-12;
- bool isZdirection = false;
- std::vector<complex_vector_type> zVectors{};
- std::vector<complex_vector_type> xyVectors{};
-
- for(unsigned int i = 0; i < 6; i++){
- complex_t z = R(i,0);
- if(std::abs(z) < threshold) z = 0.;
- if(z == 0.) isZdirection = true;
- complex_vector_type v(6);
- if(isZdirection){
- for(unsigned int j = 0;j < 6; j++){
- complex_t z = R(i,j);
- v(j) = z;
- }
- zVectors.push_back(v);
- }
- else{
- for(unsigned int j = 0; j < 6; j++){
- complex_t z = R(i,j);
- v(j) = z;
- }
- xyVectors.push_back(v);
- }
- isZdirection = false;
- }
-
- //if z-direction not found, then the system does not behave as expected
- if(zVectors.size() != 2)
- throw OpalException("SigmaGenerator::eigsolve_m()",
- "Couldn't find the vertical Eigenvectors.");
+// "Couldn't find the vertical Eigenvectors.");
+//
+// // Norms the Eigenvectors
+// for(unsigned int i = 0; i < 4; i++){
+// double norm{0};
+// for(unsigned int j = 0; j < 6; j++) norm += std::norm(xyVectors[i](j));
+// for(unsigned int j = 0; j < 6; j++) xyVectors[i](j) /= std::sqrt(norm);
+// }
+// for(unsigned int i = 0; i < 2; i++){
+// double norm{0};
+// for(unsigned int j = 0; j < 6; j++) norm += std::norm(zVectors[i](j));
+// for(unsigned int j = 0; j < 6; j++) zVectors[i](j) /= std::sqrt(norm);
+// }
+//
+// for(double i = 0; i < 6; i++){
+// R(i,0) = xyVectors[0](i);
+// R(i,1) = xyVectors[1](i);
+// R(i,2) = zVectors[0](i);
+// R(i,3) = zVectors[1](i);
+// R(i,4) = xyVectors[2](i);
+// R(i,5) = xyVectors[3](i);
+// }
+//
+// gsl_vector_complex_free(evals);
+// gsl_matrix_complex_free(evecs);
+// gsl_eigen_nonsymmv_free(wspace);
+// gsl_matrix_free(M);
+// }
- // Norms the Eigenvectors
- for(unsigned int i = 0; i < 4; i++){
- double norm{0};
- for(unsigned int j = 0; j < 6; j++) norm += std::norm(xyVectors[i](j));
- for(unsigned int j = 0; j < 6; j++) xyVectors[i](j) /= std::sqrt(norm);
- }
- for(unsigned int i = 0; i < 2; i++){
- double norm{0};
- for(unsigned int j = 0; j < 6; j++) norm += std::norm(zVectors[i](j));
- for(unsigned int j = 0; j < 6; j++) zVectors[i](j) /= std::sqrt(norm);
- }
- for(double i = 0; i < 6; i++){
- R(i,0) = xyVectors[0](i);
- R(i,1) = xyVectors[1](i);
- R(i,2) = zVectors[0](i);
- R(i,3) = zVectors[1](i);
- R(i,4) = xyVectors[2](i);
- R(i,5) = xyVectors[3](i);
- }
-
- gsl_vector_complex_free(evals);
- gsl_matrix_complex_free(evecs);
- gsl_eigen_nonsymmv_free(wspace);
- gsl_matrix_free(M);
-}
+// bool SigmaGenerator::invertMatrix_m(const sparse_matrix_type& R,
+// sparse_matrix_type& invR)
+// {
+// typedef boost::numeric::ublas::permutation_matrix<unsigned int> pmatrix_t;
+//
+// //creates a working copy of R
+// cmatrix_type A(R);
+//
+// //permutation matrix for the LU-factorization
+// pmatrix_t pm(A.size1());
+//
+// //LU-factorization
+// int res = lu_factorize(A,pm);
+//
+// if( res != 0)
+// return false;
+//
+// // create identity matrix of invR
+// invR.assign(boost::numeric::ublas::identity_matrix<complex_t>(A.size1()));
+//
+// // backsubstitute to get the inverse
+// boost::numeric::ublas::lu_substitute(A, pm, invR);
+//
+// return true;
+// }
-bool SigmaGenerator::invertMatrix_m(const sparse_matrix_type& R,
- sparse_matrix_type& invR)
+/*void*/
+typename SigmaGenerator::vector_type
+SigmaGenerator::decouple(const matrix_type& Mturn,
+ sparse_matrix_type& R,
+ sparse_matrix_type& invR)
{
- typedef boost::numeric::ublas::permutation_matrix<unsigned int> pmatrix_t;
-
- //creates a working copy of R
- cmatrix_type A(R);
-
- //permutation matrix for the LU-factorization
- pmatrix_t pm(A.size1());
+ // copy one turn matrix
+ matrix_type M(Mturn);
- //LU-factorization
- int res = lu_factorize(A,pm);
+ // reduce 6x6 matrix to 4x4 matrix
+ reduce<matrix_type>(M);
- if( res != 0)
- return false;
+ // compute symplex part
+ matrix_type Ms = rdm_m.symplex(M);
- // create identity matrix of invR
- invR.assign(boost::numeric::ublas::identity_matrix<complex_t>(A.size1()));
+ // diagonalize and compute transformation matrices
+ rdm_m.diagonalize(Ms,R,invR);
- // backsubstitute to get the inverse
- boost::numeric::ublas::lu_substitute(A, pm, invR);
-
- return true;
-}
+ /*
+ * formula (38) in paper of Dr. Christian Baumgarten:
+ * Geometrical method of decoupling
+ *
+ * [0, alpha, 0, 0;
+ * F_{d} = -beta, 0, 0, 0;
+ * 0, 0, 0, gamma;
+ * 0, 0, -delta, 0]
+ *
+ *
+ */
+ vector_type eigen(4);
+ eigen(0) = Ms(0,1); // alpha
+ eigen(1) = - Ms(1,0); // beta
+ eigen(2) = Ms(2,3); // gamma
+ eigen(3) = - Ms(3,2); // delta
+ return eigen;
-void SigmaGenerator::decouple(const matrix_type& Mturn,
- sparse_matrix_type& R,
- sparse_matrix_type& invR)
-{
- this->eigsolve_m(Mturn, R);
+// this->eigsolve_m(Mturn, R);
- if ( !this->invertMatrix_m(R, invR) )
- throw OpalException("SigmaGenerator::decouple()",
- "Couldn't compute inverse matrix!");
+// if ( !this->invertMatrix_m(R, invR) )
+// throw OpalException("SigmaGenerator::decouple()",
+// "Couldn't compute inverse matrix!");
}
@@ -744,61 +786,211 @@ void SigmaGenerator::initialize(double nuz, double ravg)
+// typename SigmaGenerator::matrix_type
+// SigmaGenerator::updateInitialSigma(const matrix_type& /*M*/,
+// sparse_matrix_type& R,
+// sparse_matrix_type& invR)
+// {
typename SigmaGenerator::matrix_type
-SigmaGenerator::updateInitialSigma(const matrix_type& /*M*/,
+SigmaGenerator::updateInitialSigma(const matrix_type& M,
+ const vector_type& eigen,
sparse_matrix_type& R,
sparse_matrix_type& invR)
{
+ /*
+ * This function is based on the paper of Dr. Christian Baumgarten:
+ * Transverse-Longitudinal Coupling by Space Charge in Cyclotrons (2012)
+ */
+
/*
* Function input:
* - M: one turn transfer matrix
+ * - eigen = {alpha, beta, gamma, delta}
* - R: transformation matrix (in paper: E)
* - invR: inverse transformation matrix (in paper: E^{-1}
*/
+// * with diagonal matrix D (stores eigenvalues of sigma*S (emittances apart from +- i),
+// * skew-symmetric matrix (formula (13)), and tranformation matrices E, E^{-1}
+// */
+
+ // normalize emittances
+ double invbg = 1.0 / (beta_m * gamma_m);
+ double ex = emittance_m[0] * invbg;
+ double ey = emittance_m[1] * invbg;
+ double ez = emittance_m[2] * invbg;
+
+ std::cout << "Eigen: " << eigen << std::endl; std::cin.get();
- /* formula (18):
- * sigma = -E*D*E^{-1}*S
- * with diagonal matrix D (stores eigenvalues of sigma*S (emittances apart from +- i),
- * skew-symmetric matrix (formula (13)), and tranformation matrices E, E^{-1}
+
+ // alpha^2-beta*gamma = 1
+
+ /* 0 eigen(0) 0 0
+ * eigen(1) 0 0 0
+ * 0 0 0 eigen(2)
+ * 0 0 eigen(3) 0
+ *
+ * M = cos(mux)*[1, 0; 0, 1] + sin(mux)*[alpha, beta; -gamma, -alpha], Book, p. 242
+ *
+ * -----------------------------------------------------------------------------------
+ * X-DIRECTION and L-DIRECTION
+ * -----------------------------------------------------------------------------------
+ * --> eigen(0) = sin(mux)*betax
+ * --> eigen(1) = -gammax*sin(mux)
+ *
+ * What is sin(mux)? --> alphax = 0 --> -alphax^2+betax*gammax = betax*gammax = 1
+ *
+ * Convention: betax > 0
+ *
+ * 1) betax = 1/gammax
+ * 2) eig0 = sin(mux)*betax
+ * 3) eig1 = -gammax*sin(mux)
+ *
+ * eig0 = sin(mux)/gammax
+ * eig1 = -gammax*sin(mux) <--> 1/gammax = -sin(mux)/eig1
+ *
+ * eig0 = -sin(mux)^2/eig1 --> -sin(mux)^2 = eig0*eig1 --> sin(mux) = sqrt(-eig0*eig1)
+ * --> gammax = -eig1/sin(mux)
+ * --> betax = eig0/sin(mux)
*/
- cmatrix_type S = boost::numeric::ublas::zero_matrix<complex_t>(6,6);
+
+ // x-direction
+ //double alphax = 0.0;
+ double betax = std::sqrt(std::fabs(eigen(0) / eigen(1)));
+ double gammax = 1.0 / betax;
+
+ // l-direction
+ //double alphal = 0.0;
+ double betal = std::sqrt(std::fabs(eigen(2) / eigen(3)));
+ double gammal = 1.0 / betal;
+
+ /*
+ * -----------------------------------------------------------------------------------
+ * Z-DIRECTION
+ * -----------------------------------------------------------------------------------
+ *
+ * m22 m23
+ * m32 m33
+ *
+ * m22 = cos(muz) + alpha*sin(muz)
+ * m33 = cos(muz) - alpha*sin(muz)
+ *
+ * --> cos(muz) = 0.5*(m22 + m33)
+ * sin(muz) = sign(m32)*sqrt(1-cos(muz)^2)
+ *
+ * m22-m33 = 2*alpha*sin(muz) --> alpha = 0.5*(m22-m33)/sin(muz)
+ *
+ * m23 = betaz*sin(muz) --> betaz = m23/sin(muz)
+ * m32 = -gammaz*sin(muz) --> gammaz = -m32/sin(muz)
+ */
+
+ double cosz = 0.5 * (M(2,2) + M(3,3));
+
+ double sign = (std::signbit(M(2,3))) ? double(-1) : double(1);
+
+ double invsinz = sign / std::sqrt(std::fabs( 1.0 - cosz * cosz));
+
+ double alphaz = 0.5 * (M(2,2) - M(3,3)) * invsinz;
+ double betaz = M(2,3) * invsinz;
+ double gammaz = - M(3,2) * invsinz;
+
+ // Convention beta>0
+ if (std::signbit(betaz)) // singbit = true if beta<0, else false
+ betaz *= -1.0;
+
+ // diagonal matrix with eigenvalues
+ matrix_type D = boost::numeric::ublas::zero_matrix<double>(6,6);
+ // x-direction
+ D(0,1) = betax * ex;
+ D(1,0) = - gammax * ex;
+ // z-direction
+ D(2,2) = alphaz * ey;
+ D(3,3) = - alphaz * ey;
+ D(2,3) = betaz * ey;
+ D(3,2) = - gammaz * ey;
+ // l-direction
+ D(4,5) = betal * ez;
+ D(5,4) = - gammal * ez;
+
+ // expand 4x4 transformation matrices to 6x6
+ expand<sparse_matrix_type>(R);
+ expand<sparse_matrix_type>(invR);
+
+ // symplectic matrix
+ sparse_matrix_type S(6,6,6);
S(0,1) = S(2,3) = S(4,5) = 1;
S(1,0) = S(3,2) = S(5,4) = -1;
- // Build new D-Matrix
- // Section 2.4 Eq. 18 in M. Frey Semester thesis
- // D = diag(i*emx,-i*emx,i*emy,-i*emy,i*emz, -i*emz)
+ matrix_type sigma = matt_boost::gemmm<matrix_type>(-invR,D,R);
+ sigma = boost::numeric::ublas::prod(sigma,S);
+ if (write_m) {
+ std::string energy = float2string(E_m);
- cmatrix_type D = boost::numeric::ublas::zero_matrix<complex_t>(6,6);
- double invbg = 1.0 / (beta_m * gamma_m);
- complex_t im(0,1);
- for(unsigned int i = 0; i < 3; ++i){
- D(2*i, 2*i) = emittance_m[i] * invbg * im;
- D(2*i+1, 2*i+1) = -emittance_m[i] * invbg * im;
- }
-
- // Computing of new Sigma
- // sigma = -R*D*R^{-1}*S
- cmatrix_type csigma(6, 6);
- csigma = boost::numeric::ublas::prod(invR, S);
- csigma = boost::numeric::ublas::prod(D, csigma);
- csigma = boost::numeric::ublas::prod(-R, csigma);
-
- matrix_type sigma(6,6);
- for (unsigned int i = 0; i < 6; ++i){
- for (unsigned int j = 0; j < 6; ++j){
- sigma(i,j) = csigma(i,j).real();
- }
- }
+ std::string fname = Util::combineFilePath({
+ OpalData::getInstance()->getAuxiliaryOutputDirectory(),
+ "maps",
+ "InitialSigmaPerAngleForEnergy" + energy + "MeV.dat"
+ });
- for (unsigned int i = 0; i < 6; ++i) {
- if(sigma(i,i) < 0.)
- sigma(i,i) *= -1.0;
+ std::ofstream writeInit(fname, std::ios::app);
+ writeMatrix(writeInit, sigma);
+ writeInit.close();
}
return sigma;
+
+
+// /*
+// * Function input:
+// * - M: one turn transfer matrix
+// * - R: transformation matrix (in paper: E)
+// * - invR: inverse transformation matrix (in paper: E^{-1}
+// */
+//
+// /* formula (18):
+// * sigma = -E*D*E^{-1}*S
+// * with diagonal matrix D (stores eigenvalues of sigma*S (emittances apart from +- i),
+// * skew-symmetric matrix (formula (13)), and tranformation matrices E, E^{-1}
+// */
+//
+// cmatrix_type S = boost::numeric::ublas::zero_matrix<complex_t>(6,6);
+// S(0,1) = S(2,3) = S(4,5) = 1;
+// S(1,0) = S(3,2) = S(5,4) = -1;
+//
+// // Build new D-Matrix
+// // Section 2.4 Eq. 18 in M. Frey Semester thesis
+// // D = diag(i*emx,-i*emx,i*emy,-i*emy,i*emz, -i*emz)
+//
+//
+// cmatrix_type D = boost::numeric::ublas::zero_matrix<complex_t>(6,6);
+// double invbg = 1.0 / (beta_m * gamma_m);
+// complex_t im(0,1);
+// for(unsigned int i = 0; i < 3; ++i){
+// D(2*i, 2*i) = emittance_m[i] * invbg * im;
+// D(2*i+1, 2*i+1) = -emittance_m[i] * invbg * im;
+// }
+//
+// // Computing of new Sigma
+// // sigma = -R*D*R^{-1}*S
+// cmatrix_type csigma(6, 6);
+// csigma = boost::numeric::ublas::prod(invR, S);
+// csigma = boost::numeric::ublas::prod(D, csigma);
+// csigma = boost::numeric::ublas::prod(-R, csigma);
+//
+// matrix_type sigma(6,6);
+// for (unsigned int i = 0; i < 6; ++i){
+// for (unsigned int j = 0; j < 6; ++j){
+// sigma(i,j) = csigma(i,j).real();
+// }
+// }
+//
+// for (unsigned int i = 0; i < 6; ++i) {
+// if(sigma(i,i) < 0.)
+// sigma(i,i) *= -1.0;
+// }
+//
+// return sigma;
}
@@ -940,4 +1132,4 @@ void SigmaGenerator::writeMatrix(std::ofstream& os, const matrix_type& m) {
os << m(i, j) << " ";
}
os << std::endl;
-}
\ No newline at end of file
+}
diff --git a/src/Distribution/SigmaGenerator.h b/src/Distribution/SigmaGenerator.h
index eb9c5ba1716595937e3e370c50a823080aa7c984..b091ae87b1396ee2aa3d92fd6593cc47817d8e47 100644
--- a/src/Distribution/SigmaGenerator.h
+++ b/src/Distribution/SigmaGenerator.h
@@ -44,6 +44,8 @@
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/vector.hpp>
+#include "rdm.h"
+
/// @brief This class computes the matched distribution
class SigmaGenerator
@@ -58,7 +60,7 @@ public:
typedef boost::numeric::ublas::matrix<complex_t> cmatrix_type;
/// Type for storing the sparse maps
- typedef boost::numeric::ublas::compressed_matrix<complex_t,
+ typedef boost::numeric::ublas::compressed_matrix<double /*complex_t*/,
boost::numeric::ublas::row_major
> sparse_matrix_type;
/// Type for storing vectors
@@ -114,15 +116,15 @@ public:
* @param Mturn is a 6x6 dimensional one turn transfer matrix
* @param R is the 6x6 dimensional transformation matrix (gets computed)
*/
- void eigsolve_m(const matrix_type& Mturn, sparse_matrix_type& R);
+// void eigsolve_m(const matrix_type& Mturn, sparse_matrix_type& R);
/*!
* @param R is the 6x6 dimensional transformation matrix
* @param invR is the 6x6 dimensional inverse transformation (gets computed)
* @return true if success
*/
- bool invertMatrix_m(const sparse_matrix_type& R,
- sparse_matrix_type& invR);
+// bool invertMatrix_m(const sparse_matrix_type& R,
+// sparse_matrix_type& invR);
/// Block diagonalizes the symplex part of the one turn transfer matrix
/*! It computes the transformation matrix <b>R</b> and its inverse <b>invR</b>.
@@ -131,7 +133,8 @@ public:
* @param R is the 6x6 dimensional transformation matrix (gets computed)
* @param invR is the 6x6 dimensional inverse transformation (gets computed)
*/
- void decouple(const matrix_type& Mturn, sparse_matrix_type& R, sparse_matrix_type& invR);
+ /*void*/
+ vector_type decouple(const matrix_type& Mturn, sparse_matrix_type& R, sparse_matrix_type& invR);
/// Checks if the sigma-matrix is an eigenellipse and returns the L2 error.
/*!
@@ -240,7 +243,12 @@ private:
* @param R is the transformation matrix
* @param invR is the inverse transformation matrix
*/
+// matrix_type updateInitialSigma(const matrix_type&,
+// sparse_matrix_type&,
+// sparse_matrix_type&);
+
matrix_type updateInitialSigma(const matrix_type&,
+ const vector_type&,
sparse_matrix_type&,
sparse_matrix_type&);
@@ -290,6 +298,16 @@ private:
const container_type& ds_turn);
void writeMatrix(std::ofstream&, const matrix_type&);
+
+ /// <b>RDM</b>-class member used for decoupling
+ RDM<double, unsigned int> rdm_m;
+
+
+ template<class matrix>
+ void reduce(matrix&);
+
+ template<class matrix>
+ void expand(matrix&);
};
@@ -326,4 +344,71 @@ std::array<double,3> SigmaGenerator::getEmittances() const
}};
}
+
+template<class matrix>
+void SigmaGenerator::reduce(matrix& M) {
+ /* The 6x6 matrix gets reduced to a 4x4 matrix in the following way:
+ *
+ * a11 a12 a13 a14 a15 a16
+ * a21 a22 a23 a24 a25 a26 a11 a12 a15 a16
+ * a31 a32 a33 a34 a35 a36 --> a21 a22 a25 a26
+ * a41 a42 a43 a44 a45 a46 a51 a52 a55 a56
+ * a51 a52 a53 a54 a55 a56 a61 a62 a65 a66
+ * a61 a62 a63 a64 a65 a66
+ */
+
+ // copy x- and l-direction to a 4x4 matrix_type
+ matrix_type M4x4(4,4);
+ for (unsigned int i = 0; i < 2; ++i) {
+ // upper left 2x2 [a11,a12;a21,a22]
+ M4x4(i,0) = M(i,0);
+ M4x4(i,1) = M(i,1);
+ // lower left 2x2 [a51,a52;a61,a62]
+ M4x4(i + 2,0) = M(i + 4,0);
+ M4x4(i + 2,1) = M(i + 4,1);
+ // upper right 2x2 [a15,a16;a25,a26]
+ M4x4(i,2) = M(i,4);
+ M4x4(i,3) = M(i,5);
+ // lower right 2x2 [a55,a56;a65,a66]
+ M4x4(i + 2,2) = M(i + 4,4);
+ M4x4(i + 2,3) = M(i + 4,5);
+ }
+
+ M.resize(4,4,false);
+ M.swap(M4x4);
+}
+
+template<class matrix>
+void SigmaGenerator::expand(matrix& M) {
+ /* The 4x4 matrix gets expanded to a 6x6 matrix in the following way:
+ *
+ * a11 a12 0 0 a13 a14
+ * a11 a12 a13 a14 a21 a22 0 0 a23 a24
+ * a21 a22 a23 a24 --> 0 0 1 0 0 0
+ * a31 a32 a33 a34 0 0 0 1 0 0
+ * a41 a42 a43 a44 a31 a32 0 0 a33 a34
+ * a41 a42 0 0 a43 a44
+ */
+
+ matrix M6x6 = boost::numeric::ublas::identity_matrix<double>(6,6);
+
+ for (unsigned int i = 0; i < 2; ++i) {
+ // upper left 2x2 [a11,a12;a21,a22]
+ M6x6(i,0) = M(i,0);
+ M6x6(i,1) = M(i,1);
+ // lower left 2x2 [a31,a32;a41,a42]
+ M6x6(i + 4,0) = M(i + 2,0);
+ M6x6(i + 4,1) = M(i + 2,1);
+ // upper right 2x2 [a13,a14;a23,a24]
+ M6x6(i,4) = M(i,2);
+ M6x6(i,5) = M(i,3);
+ // lower right 2x2 [a22,a34;a43,a44]
+ M6x6(i + 4,4) = M(i + 2,2);
+ M6x6(i + 4,5) = M(i + 2,3);
+ }
+
+ // exchange
+ M.swap(M6x6);
+}
+
#endif
\ No newline at end of file
diff --git a/src/Distribution/rdm.h b/src/Distribution/rdm.h
new file mode 100644
index 0000000000000000000000000000000000000000..bbca9bd993dcd16f329918b9de33eda26a919b6b
--- /dev/null
+++ b/src/Distribution/rdm.h
@@ -0,0 +1,304 @@
+/**
+ * @file rdm.h
+ * They're ordered after the paper of Dr. C. Baumgarten: "Use of real Dirac matrices in two-dimensional coupled linear optics".
+ * The diagonalizing method is based on the paper "Geometrical method of decoupling" (2012) of Dr. C. Baumgarten.
+ * The template parameter Value_type is the datatype of the matrices and variables.
+ *
+ * @author Matthias Frey
+ * @version 1.0
+ */
+#ifndef RDM_H
+#define RDM_H
+
+#include <cmath>
+#include <iostream>
+#include <stdexcept>
+
+#include <boost/numeric/ublas/matrix_sparse.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+
+#include "matrix_vector_operation.h"
+
+// Remark: Comments starting with "///" are doxygen comments, comments with "//" are normal comments.
+
+/// @brief Real Dirac Matrix class
+template<typename Value_type, typename Size_type>
+class RDM
+{
+ public:
+ /// Type of variables
+ typedef Value_type value_type;
+ /// Type for specifying sizes
+ typedef Size_type size_type;
+ /// Sparse matrix type definition
+ typedef boost::numeric::ublas::compressed_matrix<value_type,boost::numeric::ublas::row_major> sparse_matrix_type;
+ /// Dense matrix type definition
+ typedef boost::numeric::ublas::matrix<value_type> matrix_type;
+ /// Dense vector type definition
+ typedef boost::numeric::ublas::vector<value_type> vector_type;
+
+ /// Default constructor (sets only NumOfRDMs and DimOfRDMs)
+ RDM();
+
+ /// Returns the i-th Real Dirac matrix
+ /*!
+ * @param i specifying the matrix (has to be in the range from 0 to 15)
+ */
+ sparse_matrix_type getRDM(short);
+
+ /// Decomposes a real-valued 4x4 matrix into a linear combination and returns a vector containing the coefficients
+ /*!
+ * @param M an arbitrary real-valued 4x4 matrix
+ */
+ vector_type decompose(const matrix_type&);
+
+ /// Takes a vector of coefficients, evaluates the linear combination of RDMs with these coefficients and returns a 4x4 matrix
+ /*!
+ * @param coeffs is a vector of coefficients (at most length NumOfRDMs)
+ */
+ matrix_type combine(const vector_type&);
+
+ /// Brings a 4x4 symplex matrix into Hamilton form and computes the transformation matrix and its inverse
+ /*!
+ * @param Ms is a 4x4 symplex matrix
+ * @param R is the 4x4 transformation matrix (gets computed)
+ * @param invR is the 4x4 inverse transformation matrix (gets computed)
+ */
+ void diagonalize(matrix_type&, sparse_matrix_type&, sparse_matrix_type&);
+
+ /// Returns the symplex part of a 4x4 real-valued matrix
+ /*!
+ * @param M 4x4 real-valued matrix
+ */
+ matrix_type symplex(const matrix_type&);
+
+ /// Returns the cosymplex part of a 4x4 real-valued matrix
+ /*!
+ * @param M 4x4 real-valued matrix
+ */
+ matrix_type cosymplex(const matrix_type&);
+
+ /// The number of real Dirac matrices
+ short NumOfRDMs;
+ /// The matrix dimension (4x4)
+ short DimOfRDMs;
+
+ private:
+ /// Applies a rotation to the matrix M by a given angle
+ /*!
+ * @param M is the matrix to be transformed
+ * @param i is the i-th RDM used for transformation
+ * @param phi is the angle of rotation
+ * @param Rtot is a reference to the current transformation matrix
+ * @param invRtot is a reference to the inverse of the current transformation matrix
+ */
+ void transform(matrix_type&, short, value_type, sparse_matrix_type&, sparse_matrix_type&);
+};
+
+// -----------------------------------------------------------------------------------------------------------------------
+// PUBLIC MEMBER FUNCTIONS
+// -----------------------------------------------------------------------------------------------------------------------
+
+template<typename Value_type, typename Size_type>
+RDM<Value_type, Size_type>::RDM() : NumOfRDMs(16), DimOfRDMs(4) {};
+
+template<typename Value_type, typename Size_type>
+typename RDM<Value_type, Size_type>::sparse_matrix_type RDM<Value_type, Size_type>::getRDM(short i) {
+ sparse_matrix_type rdm(4,4,4); // #nrows, #ncols, #non-zeros
+ switch(i) {
+ case 0: rdm(0,1) = rdm(2,3) = 1; rdm(1,0) = rdm(3,2) = -1; break;
+ case 1: rdm(0,1) = rdm(1,0) = -1; rdm(2,3) = rdm(3,2) = 1; break;
+ case 2: rdm(0,3) = rdm(1,2) = rdm(2,1) = rdm(3,0) = 1; break;
+ case 3: rdm(0,0) = rdm(2,2) = -1; rdm(1,1) = rdm(3,3) = 1; break;
+ case 4: rdm(0,0) = rdm(3,3) = -1; rdm(1,1) = rdm(2,2) = 1; break;
+ case 5: rdm(0,2) = rdm(2,0) = 1; rdm(1,3) = rdm(3,1) = -1; break;
+ case 6: rdm(0,1) = rdm(1,0) = rdm(2,3) = rdm(3,2) = 1; break;
+ case 7: rdm(0,3) = rdm(2,1) = 1; rdm(1,2) = rdm(3,0) = -1; break;
+ case 8: rdm(0,1) = rdm(3,2) = 1; rdm(1,0) = rdm(2,3) = -1; break;
+ case 9: rdm(0,2) = rdm(1,3) = -1; rdm(2,0) = rdm(3,1) = 1; break;
+ case 10: rdm(0,2) = rdm(3,1) = 1; rdm(1,3) = rdm(2,0) = -1; break;
+ case 11: rdm(0,2) = rdm(1,3) = rdm(2,0) = rdm(3,1) = -1; break;
+ case 12: rdm(0,0) = rdm(1,1) = -1; rdm(2,2) = rdm(3,3) = 1; break;
+ case 13: rdm(0,3) = rdm(3,0) = -1; rdm(1,2) = rdm(2,1) = 1; break;
+ case 14: rdm(0,3) = rdm(1,2) = -1; rdm(2,1) = rdm(3,0) = 1; break;
+ case 15: rdm(0,0) = rdm(1,1) = rdm(2,2) = rdm(3,3) = 1; break;
+ default: throw std::out_of_range("Error in RDM::generate: 0 <= i <= 15"); break;
+ }
+ return rdm;
+}
+
+template<typename Value_type, typename Size_type>
+typename RDM<Value_type, Size_type>::vector_type RDM<Value_type, Size_type>::decompose(const matrix_type& M) {
+ /*
+ * Formula (11) from paper:
+ * Geometrical method of decoupling
+ */
+ if (M.size1() != 4 && M.size1() != M.size2())
+ throw std::length_error("Error in RDM::decompose: Wrong matrix dimensions.");
+
+ vector_type coeffs(16);
+
+ matrix_type left(4,4), right(4,4), rdm2(4,4);
+
+ value_type inv32 = 1.0 / 32.0;
+
+ // do it with iterators
+ for (short i = 0; i < 16; ++i) {
+ sparse_matrix_type rdm = getRDM(i);
+ left = boost::numeric::ublas::prod(M,rdm);
+ right = boost::numeric::ublas::prod(rdm,M);
+ rdm2 = boost::numeric::ublas::prod(rdm,rdm);
+ left *= inv32;
+ right *= inv32;
+ left += right;
+ coeffs(i) = matt_boost::trace(rdm2) * matt_boost::trace(left);
+ }
+ return coeffs;
+}
+
+template<typename Value_type, typename Size_type>
+typename RDM<Value_type, Size_type>::matrix_type RDM<Value_type, Size_type>::combine(const vector_type& coeffs) {
+ if (coeffs.size() > 16)
+ throw std::length_error("Error in RDM::combine: Wrong size of coefficient vector.");
+
+ // initialize a 4x4 zero matrix
+ matrix_type M = boost::numeric::ublas::zero_matrix<value_type>(4,4);
+
+ // evaluate linear combination
+ for (short i = 0; i < 16; ++i)
+ M += coeffs(i)*getRDM(i);
+
+ return M;
+}
+
+template<typename Value_type, typename Size_type>
+void RDM<Value_type, Size_type>::diagonalize(matrix_type& Ms, sparse_matrix_type& R, sparse_matrix_type& invR) {
+
+ // R and invR store the total transformation
+ R = boost::numeric::ublas::identity_matrix<value_type>(4);
+ invR = R;
+
+ vector_type P(3), E(3), B(3), b;
+ value_type mr, mg, mb, eps;
+
+ // Lambda function to compute vectors E, P, B and scalar eps (it takes the current Ms as reference argument (--> [&])
+ auto mult = [&](short i) {
+ /*
+ * For computing E, P, B, eps according to formula (C4) from paper:
+ * Geometrical method of decoupling
+ */
+ matrix_type tmp = boost::numeric::ublas::prod(Ms,getRDM(i))+boost::numeric::ublas::prod(getRDM(i),Ms);
+ return 0.125*matt_boost::trace(tmp);
+ };
+
+ // 1. Transformation with \gamma_{0}
+ P(0) = mult(1); P(1) = mult(2); P(2) = mult(3);
+ E(0) = mult(4); E(1) = mult(5); E(2) = mult(6);
+ B(0) = -mult(7); B(1) = -mult(8); B(2) = -mult(9);
+ mr = boost::numeric::ublas::inner_prod(E,B); // formula (31), paper: Geometrical method of decoupling
+ mg = boost::numeric::ublas::inner_prod(B,P); // formula (31), paper: Geometrical method of decoupling
+
+ transform(Ms,short(0),0.5 * std::atan2(mg,mr),R,invR);
+
+ // 2. Transformation with \gamma_{7}
+ eps = -mult(0);
+ P(0) = mult(1); P(1) = mult(2); P(2) = mult(3);
+ E(0) = mult(4); E(1) = mult(5); E(2) = mult(6);
+ B(0) = -mult(7); B(1) = -mult(8); B(2) = -mult(9);
+ b = eps*B+matt_boost::cross_prod(E,P); // formula (32), paper: Geometrical method of decoupling
+
+ transform(Ms,short(7),0.5 * std::atan2(b(2),b(1)),R,invR);
+
+ // 3. Transformation with \gamma_{9}
+ eps = -mult(0);
+ P(0) = mult(1); P(1) = mult(2); P(2) = mult(3);
+ E(0) = mult(4); E(1) = mult(5); E(2) = mult(6);
+ B(0) = -mult(7); B(1) = -mult(8); B(2) = -mult(9);
+ b = eps*B+matt_boost::cross_prod(E,P);
+
+ transform(Ms,short(9),-0.5 * std::atan2(b(0),b(1)),R,invR);
+
+ // 4. Transformation with \gamma_{2}
+ eps = -mult(0);
+ P(0) = mult(1); P(1) = mult(2); P(2) = mult(3);
+ E(0) = mult(4); E(1) = mult(5); E(2) = mult(6);
+ B(0) = -mult(7); B(1) = -mult(8); B(2) = -mult(9);
+ mr = boost::numeric::ublas::inner_prod(E,B);
+ b = eps*B+matt_boost::cross_prod(E,P);
+
+ // Transformation distinction made according to function "rdm_Decouple_F" in rdm.c of Dr. Christian Baumgarten
+
+ if (std::fabs(mr) < std::fabs(b(1))) {
+ transform(Ms,short(2),0.5 * std::atanh(mr / b(1)),R,invR);
+ } else {
+ transform(Ms,short(2),0.5 * std::atanh(b(1) / mr),R,invR);
+ }
+
+ eps = -mult(0);
+ P(0) = mult(1); P(1) = mult(2); P(2) = mult(3);
+ E(0) = mult(4); E(1) = mult(5); E(2) = mult(6);
+ B(0) = -mult(7); B(1) = -mult(8); B(2) = -mult(9);
+
+ // formula (31), paper: Geometrical method of decoupling
+ mr = boost::numeric::ublas::inner_prod(E,B);
+ mg = boost::numeric::ublas::inner_prod(B,P);
+ mb = boost::numeric::ublas::inner_prod(E,P);
+
+ value_type P2 = boost::numeric::ublas::inner_prod(P,P);
+ value_type E2 = boost::numeric::ublas::inner_prod(E,E);
+
+ // 5. Transform with \gamma_{0}
+ transform(Ms,short(0),0.25 * std::atan2(mb,0.5 * (E2 - P2)),R,invR);
+
+ // 6. Transformation with \gamma_{8}
+ P(0) = mult(1); P(2) = mult(3);
+
+ transform(Ms,short(8),-0.5 * std::atan2(P(2),P(0)),R,invR);
+}
+
+template<typename Value_type, typename Size_type>
+typename RDM<Value_type, Size_type>::matrix_type RDM<Value_type, Size_type>::symplex(const matrix_type& M) {
+ /*
+ * formula(16), p. 3
+ */
+ sparse_matrix_type rdm0 = getRDM(0);
+ return 0.5 * (M + matt_boost::gemmm<matrix_type>(rdm0,boost::numeric::ublas::trans(M),rdm0));
+}
+
+
+template<typename Value_type, typename Size_type>
+typename RDM<Value_type, Size_type>::matrix_type RDM<Value_type, Size_type>::cosymplex(const matrix_type& M) {
+ /*
+ * formula (16), p, 3
+ */
+ sparse_matrix_type rdm0 = getRDM(0);
+ return 0.5 * (M - matt_boost::gemmm<matrix_type>(rdm0,boost::numeric::ublas::trans(M),rdm0));
+}
+
+// -----------------------------------------------------------------------------------------------------------------------
+// PRIVATE MEMBER FUNCTIONS
+// -----------------------------------------------------------------------------------------------------------------------
+
+template<typename Value_type, typename Size_type>
+void RDM<Value_type, Size_type>::transform(matrix_type& M, short i, value_type phi, sparse_matrix_type& Rtot, sparse_matrix_type& invRtot) {
+
+ if (phi) { // if phi == 0 --> nothing happens, since R and invR would be identity_matrix matrix
+ sparse_matrix_type R(4,4), invR(4,4);
+ sparse_matrix_type I = boost::numeric::ublas::identity_matrix<value_type>(4);
+
+ if (i < 7 && i != 0 && i < 11 && i != 14) {
+ R = I * std::cosh(phi) + getRDM(i) * std::sinh(phi);
+ invR = I * std::cosh(phi) - getRDM(i) * std::sinh(phi);
+ } else {
+ R = I * std::cos(phi) + getRDM(i) * std::sin(phi);
+ invR = I * std::cos(phi) - getRDM(i) * std::sin(phi);
+ }
+ // update matrices
+ M = matt_boost::gemmm<matrix_type>(R,M,invR);
+ Rtot = boost::numeric::ublas::prod(R,Rtot);
+ invRtot = boost::numeric::ublas::prod(invRtot,invR);
+ }
+}
+
+#endif
\ No newline at end of file