25 #include <gsl/gsl_math.h>
26 #include <gsl/gsl_eigen.h>
27 #include <gsl/gsl_linalg.h>
31 : mapAnalysis_m(
IpplTimings::getTimer(
"mapAnalysis"))
32 , bunchAnalysis_m(
IpplTimings::getTimer(
"bunchAnalysis"))
38 tplot.open (
"tunePlot.txt",std::ios::app);
39 tplot << std::setprecision(16);
45 fMatrix_t blocktMatrix, scalingMatrix, invScalingMatrix;
52 for (
int i=0; i <6 ; i++){
55 for (
int j=0; j <6 ; j+=2){
56 temp += 2*(eigenVecM[j][i] *
std::conj(eigenVecM[j+1][i])).
imag();
63 for (
int j=0; j <6 ; j++){
72 for (
int i=0; i <6 ; i++){
73 for (
int j=0; j <6 ; j++){
74 eigenVecMT[i][j]= eigenVecM[j][i];
81 std::array<int,6> pairs;
88 for (
int i=0; i<6;i++) pairs[i]=-1;
91 for (
int i=0; i <6 ; i++){
92 for (
int j=0; j <6 ; j++){
93 if (K[i][j]>1-prec && K[i][j]<1+prec){
109 for (
int i=0; i<6; i++) {
110 pairsIt =
std::find (pairs.begin(), pairs.end(), i);
111 if (pairsIt != pairs.end()){
115 for (
int j=0; i!=j && j <6 ; j++){
116 double diff =
std::abs( eigenValM[i][i] - eigenValM[j][j]);
126 for (
int i=0; i < 6 && idx < 6; i++){
127 pairsIt =
std::find (pairs.begin(), pairs.end(), i);
129 if (pairsIt != pairs.end()){
142 for (
int i=0; i <6 ; i++){
143 eigenValM[i][i]=tempValM[pairs[i]][pairs[i]];
144 for (
int j=0; j <6 ; j++){
145 eigenVecM[j][i]=tempM[j][pairs[i]];
156 for (
int i=0; i<6; i++){
166 for (
int i = 0; i < 3; i++){
175 double lenEigenVec = 0;
176 for (
int j = 0; j < 3; j++){
177 lenEigenVec +=
std::abs(eigenVecM[i][j]);
184 tplot<<
"3: ("<<betaTunes3[i]<<
",0)"<<
std::endl;
199 for (
int i=0; i <6 ; i=i+2){
200 skewMatrix[0+i][1+i]=1.;
201 skewMatrix[1+i][0+i]=-1.;
204 sigMatrix=sigMatrix*skewMatrix;
206 cfMatrix_t eigenValM, eigenVecM, invEigenVecM;
213 for (
int i=0; i<6; i++){
214 for (
int j =0; j<6; j++){
215 sigmaBlockM[i][j]=cblocksigM[i][j].real();
220 cfMatrix_t teigenValM, teigenVecM, tinvEigenVecM;
221 eigenDecomp_m(sigmaBlockM, teigenValM, teigenVecM, tinvEigenVecM);
223 for (
int i=0; i<6; i++){
238 for (
int i = 0; i < 6; i++) {
239 for (
int j = 0; j < 6; j++) {
246 gsl_matrix_view m = gsl_matrix_view_array(data, 6, 6);
247 gsl_vector_complex *eval = gsl_vector_complex_alloc(6);
248 gsl_matrix_complex *evec = gsl_matrix_complex_alloc(6, 6);
249 gsl_matrix_complex *eveci = gsl_matrix_complex_alloc(6, 6);
250 gsl_permutation * p = gsl_permutation_alloc(6);
251 gsl_eigen_nonsymmv_workspace * w = gsl_eigen_nonsymmv_alloc(6);
254 gsl_eigen_nonsymmv(&m.matrix, eval, evec, w);
255 gsl_eigen_nonsymmv_free(w);
256 gsl_eigen_nonsymmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
258 for (
int i = 0; i < 6; ++i) {
259 eigenVal[i][i] = std::complex<double>(
260 GSL_REAL(gsl_vector_complex_get(eval, i)),
261 GSL_IMAG(gsl_vector_complex_get(eval, i)));
262 for (
int j = 0; j < 6; ++j) {
263 eigenVec[i][j] = std::complex<double>(
264 GSL_REAL(gsl_matrix_complex_get(evec, i, j)),
265 GSL_IMAG(gsl_matrix_complex_get(evec, i, j)));
270 gsl_linalg_complex_LU_decomp(evec, p, &s);
271 gsl_linalg_complex_LU_invert(evec, p, eveci);
274 for (
int i = 0; i < 6; ++i) {
275 for (
int j = 0; j < 6; ++j) {
276 invEigenVec[i][j] = std::complex<double>(
277 GSL_REAL(gsl_matrix_complex_get(eveci, i, j)),
278 GSL_IMAG(gsl_matrix_complex_get(eveci, i, j)));
284 gsl_vector_complex_free(eval);
285 gsl_matrix_complex_free(evec);
286 gsl_matrix_complex_free(eveci);
296 cfMatrix_t cM, qMatrix, invqMatrix, nMatrix, invnMatrix, rMatrix;
298 for (
int i=0; i<6; i++){
299 for (
int j =0; j<6; j++){
300 cM[i][j]=std::complex<double>(M[i][j],0);
304 for (
int i=0; i <6 ; i=i+2){
305 qMatrix[0+i][0+i]=std::complex<double>(1.,0);
306 qMatrix[0+i][1+i]=std::complex<double>(0,1.);
307 qMatrix[1+i][0+i]=std::complex<double>(1.,0);
308 qMatrix[1+i][1+i]=std::complex<double>(0,-1);
310 invqMatrix[0+i][0+i]=std::complex<double>(1.,0);
311 invqMatrix[0+i][1+i]=std::complex<double>(1.,0);
312 invqMatrix[1+i][0+i]=std::complex<double>(0.,-1.);
313 invqMatrix[1+i][1+i]=std::complex<double>(0,1.);
318 nMatrix=eigenVecM*qMatrix;
319 invnMatrix= invqMatrix* invEigenVecM;
322 rMatrix= invnMatrix * cM * nMatrix;
333 for (
int i=0; i<6; i++){
334 for (
int j =0; j<6; j++){
335 ctM[i][j]=std::complex<double>(tM[i][j],0);
344 std::array<double, 3> phi;
346 for (
int i = 0; i < 3; i++){
353 for (
int i=0; i<6; i++){
354 for (
int j =0; j<6; j++){
355 cR[i][j]=std::complex<double>(R1[i][j],0);
356 N1[i][j]=oldN[i][j].real();
364 cfMatrix_t eigenValM, eigenVecM, invEigenVecM, eigenVecMT;
369 cfMatrix_t cM, qMatrix, invqMatrix, nMatrix, invnMatrix, rMatrix;
376 for (
int i=0; i<6; i++){
377 for (
int j =0; j<6; j++){
378 cM[i][j]=std::complex<double>(M[i][j],0);
382 for (
int i=0; i <6 ; i=i+2){
383 qMatrix[0+i][0+i]=std::complex<double>(1.,0);
384 qMatrix[0+i][1+i]=std::complex<double>(0,1.);
385 qMatrix[1+i][0+i]=std::complex<double>(1.,0);
386 qMatrix[1+i][1+i]=std::complex<double>(0,-1);
388 invqMatrix[0+i][0+i]=std::complex<double>(1.,0);
389 invqMatrix[0+i][1+i]=std::complex<double>(1.,0);
390 invqMatrix[1+i][0+i]=std::complex<double>(0.,-1.);
391 invqMatrix[1+i][1+i]=std::complex<double>(0,1.);
398 invN= invqMatrix* invEigenVecM;
405 for (
int i = 0; i < 3; i++){
407 R[2*i+1][2*i+1] = R[2*i][2*i];
409 R[2*i+1][2*i] = -R[2*i][2*i+1];
418 for (
int i = 0; i < 3; i++){
428 for (
int i=0; i<6; i++){
429 for (
int j =0; j<6; j++){
430 M[i][j]=cM[i][j].real();
438 for (
int i=0; i<6; i++){
439 for (
int j =0; j<6; j++){
440 M[i][j]=cM[i][j].imag();
448 for (
int i=0; i<6; i++){
449 for (
int j =0; j<6; j++){
450 cM[i][j]=std::complex<double>(M[i][j],0);
459 gsl_set_error_handler_off();
461 gsl_matrix_complex *m = gsl_matrix_complex_alloc(6, 6);
462 gsl_matrix_complex *invm = gsl_matrix_complex_alloc(6, 6);
463 gsl_permutation * p = gsl_permutation_alloc(6);
469 for (
int i = 0; i < 6; ++i) {
470 for (
int j = 0; j < 6; ++j) {
472 gsl_matrix_complex_set( m, i, j, temp);
477 int eigenDecompStatus = gsl_linalg_complex_LU_decomp(m, p, &s);
478 if (eigenDecompStatus != 0){
479 std::cout<<
"This actually works" <<
std::endl;
484 int invertStatus = gsl_linalg_complex_LU_invert(m, p, invm);
486 if ( invertStatus ) {
492 if (invertStatus != 0){
493 std::cout<<
"This actually works" <<
std::endl;
500 for (
int i = 0; i < 6; ++i) {
501 for (
int j = 0; j < 6; ++j) {
502 invM[i][j] = std::complex<double>(
503 GSL_REAL(gsl_matrix_complex_get(invm, i, j)),
504 GSL_IMAG(gsl_matrix_complex_get(invm, i, j)));
509 gsl_matrix_complex_free(m);
510 gsl_matrix_complex_free(invm);
511 gsl_permutation_free(p);
IpplTimings::TimerRef bunchAnalysis_m
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
PETE_TUTree< FnArcCos, typename T::PETE_Expr_t > acos(const PETE_Expr< T > &l)
constexpr double e
The value of .
IpplTimings::TimerRef mapAnalysis_m
void setNMatrix_m(fMatrix_t &M, cfMatrix_t &N, cfMatrix_t &invN)
A templated representation for vectors.
PETE_TUTree< FnFabs, typename T::PETE_Expr_t > fabs(const PETE_Expr< T > &l)
Tps< T > sin(const Tps< T > &x)
Sine.
fMatrix_t realPartOfMatrix_m(cfMatrix_t cM)
fMatrix_t imagPartOfMatrix_m(cfMatrix_t cM)
Tps< T > log(const Tps< T > &x)
Natural logarithm.
FLieGenerator< T, N > real(const FLieGenerator< std::complex< T >, N > &)
Take real part of a complex generator.
fMatrix_t createSkewMatrix_m()
void rearrangeEigen_m(cfMatrix_t &EigenVal, cfMatrix_t &EigenVec)
cfMatrix_t complexTypeMatrix_m(fMatrix_t M)
void linSigAnalyze(fMatrix_t &sigMatrix)
constexpr double pi
The value of .
m_complex conj(const m_complex &c)
void printPhaseShift_m(fMatrix_t &Sigma, fMatrix_t tM, cfMatrix_t &oldN)
PETE_TUTree< FnArcTan, typename T::PETE_Expr_t > atan(const PETE_Expr< T > &l)
static void startTimer(TimerRef t)
void eigenDecomp_m(const fMatrix_t &M, cfMatrix_t &eigenVal, cfMatrix_t &eigenVec, cfMatrix_t &invEigenVec)
fMatrix_t createRotMatrix_m(std::array< double, 3 > phi)
const T * find(const T table[], const std::string &name)
Look up name.
Tps< T > sqrt(const Tps< T > &x)
Square root.
cfMatrix_t invertMatrix_m(const cfMatrix_t &M)
Tps< T > cos(const Tps< T > &x)
Cosine.
cfMatrix_t getBlockDiagonal_m(const fMatrix_t &M, cfMatrix_t &eigenVecM, cfMatrix_t &invEigenVecM)
PETE_TUTree< FnArcSin, typename T::PETE_Expr_t > asin(const PETE_Expr< T > &l)
void linTAnalyze(const fMatrix_t &tMatrix)
std::string::iterator iterator
static void stopTimer(TimerRef t)
FLieGenerator< T, N > imag(const FLieGenerator< std::complex< T >, N > &)
Take imaginary part of a complex generator.
Inform & endl(Inform &inf)