MapAnalyser.cpp 8.99 KB
Newer Older
snuverink_j's avatar
snuverink_j committed
1
//
2
// Source file of the MapAnalyser class, analyses linear maps from OPAL-map.
snuverink_j's avatar
snuverink_j committed
3
//
4 5
// Copyright (c) 2018, Philippe Ganz, ETH Zürich
// All rights reserved
snuverink_j's avatar
snuverink_j committed
6
//
7
// OPAL is licensed under GNU GPL version 3.
snuverink_j's avatar
snuverink_j committed
8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

#include "MapAnalyser.h"

#include <fstream>

#include "Physics/Physics.h"

#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>


MapAnalyser::MapAnalyser()
    : mapAnalysis_m(IpplTimings::getTimer("mapAnalysis"))
    , bunchAnalysis_m(IpplTimings::getTimer("bunchAnalysis"))
{ }

25
//Returns the eigenDecompositon for the fMatrix M = eigenVec * eigenVal * invEigenVec
snuverink_j's avatar
snuverink_j committed
26 27
void MapAnalyser::eigenDecomp_m(const fMatrix_t& M, cfMatrix_t& eigenVal, cfMatrix_t& eigenVec, cfMatrix_t& invEigenVec){

28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
    double data[36];
    int idx, s;
    for (int i = 0; i < 6; i++) {
        for (int j = 0; j < 6; j++) {
            idx = i * 6 + j;
            data[idx] = M[i][j];
        }
    }

    gsl_matrix_view m = gsl_matrix_view_array(data, 6, 6);
    gsl_vector_complex *eval = gsl_vector_complex_alloc(6);
    gsl_matrix_complex *evec = gsl_matrix_complex_alloc(6, 6);
    gsl_matrix_complex *eveci = gsl_matrix_complex_alloc(6, 6);
    gsl_permutation * p = gsl_permutation_alloc(6);
    gsl_eigen_nonsymmv_workspace * w = gsl_eigen_nonsymmv_alloc(6);

    //get Eigenvalues and Eigenvectors
    gsl_eigen_nonsymmv(&m.matrix, eval, evec, w);
    gsl_eigen_nonsymmv_free(w);
    gsl_eigen_nonsymmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);

    for (int i = 0; i < 6; ++i) {
        eigenVal[i][i] = std::complex<double>(
                                              GSL_REAL(gsl_vector_complex_get(eval, i)),
                                              GSL_IMAG(gsl_vector_complex_get(eval, i)));
        for (int j = 0; j < 6; ++j) {
            eigenVec[i][j] = std::complex<double>(
                                                  GSL_REAL(gsl_matrix_complex_get(evec, i, j)),
                                                  GSL_IMAG(gsl_matrix_complex_get(evec, i, j)));
        }
    }

    //invert Eigenvectormatrix
    gsl_linalg_complex_LU_decomp(evec, p, &s);
    gsl_linalg_complex_LU_invert(evec, p, eveci);

    //Create invEigenVecMatrix
    for (int i = 0; i < 6; ++i) {
        for (int j = 0; j < 6; ++j) {
            invEigenVec[i][j] = std::complex<double>(
                                                     GSL_REAL(gsl_matrix_complex_get(eveci, i, j)),
                                                     GSL_IMAG(gsl_matrix_complex_get(eveci, i, j)));
        }
    }

    //free space
    gsl_vector_complex_free(eval);
    gsl_matrix_complex_free(evec);
    gsl_matrix_complex_free(eveci);
snuverink_j's avatar
snuverink_j committed
77 78 79 80 81
}


//Transforms the Matirx to a block diagonal rotation Matrix
MapAnalyser::cfMatrix_t MapAnalyser::getBlockDiagonal_m(const fMatrix_t& M,
82
                                                        cfMatrix_t& eigenVecM, cfMatrix_t& invEigenVecM){
snuverink_j's avatar
snuverink_j committed
83

84
    cfMatrix_t cM, qMatrix, invqMatrix, nMatrix, invnMatrix, rMatrix;
snuverink_j's avatar
snuverink_j committed
85

86 87 88
    for (int i = 0; i < 6; i++){
        for (int j = 0; j < 6; j++){
            cM[i][j]=std::complex<double>(M[i][j], 0);
89 90
        }
    }
snuverink_j's avatar
snuverink_j committed
91

92 93 94 95 96
    for (int i = 0; i < 6; i = i+2){
        qMatrix[  i][  i] = std::complex<double>(1., 0);
        qMatrix[  i][1+i] = std::complex<double>(0, 1.);
        qMatrix[1+i][  i] = std::complex<double>(1., 0);
        qMatrix[1+i][1+i] = std::complex<double>(0, -1);
snuverink_j's avatar
snuverink_j committed
97

98 99 100 101
        invqMatrix[  i][  i] = std::complex<double>(1., 0);
        invqMatrix[  i][1+i] = std::complex<double>(1., 0);
        invqMatrix[1+i][  i] = std::complex<double>(0., -1.);
        invqMatrix[1+i][1+i] = std::complex<double>(0, 1.);
102
    }
103 104
    qMatrix /= std::sqrt(2.);
    invqMatrix /= std::sqrt(2);
snuverink_j's avatar
snuverink_j committed
105

106 107
    nMatrix = eigenVecM*qMatrix;
    invnMatrix = invqMatrix* invEigenVecM;
snuverink_j's avatar
snuverink_j committed
108 109


110
    rMatrix = invnMatrix * cM * nMatrix;
snuverink_j's avatar
snuverink_j committed
111

112
    return rMatrix;
snuverink_j's avatar
snuverink_j committed
113 114 115 116 117 118 119
}


void MapAnalyser::printPhaseShift_m(fMatrix_t& Sigma, fMatrix_t tM, cfMatrix_t& oldN){
    cfMatrix_t N1, cinvN, cR, ctM, N2;
    fMatrix_t R1, S, sigmaS;

120 121 122
    for (int i = 0; i < 6; i++){
        for (int j = 0; j < 6; j++){
            ctM[i][j] = std::complex<double>(tM[i][j], 0);
123 124
        }
    }
snuverink_j's avatar
snuverink_j committed
125 126 127 128 129 130 131 132 133 134 135 136

    S = createSkewMatrix_m();
    sigmaS = Sigma*S;

    setNMatrix_m(sigmaS, N2, cinvN);

    std::array<double, 3> phi;

    for (int i = 0; i < 3; i++){
        phi[i] = std::atan(oldN[2*i+1][i].real()/oldN[2*i][2*i].real());
    }

137
    R1 = createRotMatrix_m(phi);
snuverink_j's avatar
snuverink_j committed
138

139 140 141 142
    for (int i = 0; i < 6; i++){
        for (int j = 0; j < 6; j++){
            cR[i][j] = std::complex<double>(R1[i][j], 0);
            N1[i][j] = oldN[i][j].real();
snuverink_j's avatar
snuverink_j committed
143 144 145 146 147 148 149 150 151 152 153 154 155
        }
    }
}


void MapAnalyser::setNMatrix_m(fMatrix_t& M, cfMatrix_t& N, cfMatrix_t& invN){

    cfMatrix_t eigenValM, eigenVecM, invEigenVecM, eigenVecMT;

    eigenDecomp_m(M, eigenValM, eigenVecM, invEigenVecM);

    cfMatrix_t cM, qMatrix, invqMatrix, nMatrix, invnMatrix, rMatrix;

156 157 158
    //std::ofstream tmap;
    //tmap.open ("TransferMap.txt",std::ios::app);
    //tmap << std::setprecision(16);
snuverink_j's avatar
snuverink_j committed
159 160


161 162 163
    for (int i = 0; i < 6; i++){
        for (int j = 0; j < 6; j++){
            cM[i][j] = std::complex<double>(M[i][j], 0);
snuverink_j's avatar
snuverink_j committed
164
        }
165
    }
snuverink_j's avatar
snuverink_j committed
166

167 168 169 170 171
    for (int i = 0; i < 6; i = i+2){
        qMatrix[  i][  i] = std::complex<double>(1., 0);
        qMatrix[  i][1+i] = std::complex<double>(0, 1.);
        qMatrix[1+i][  i] = std::complex<double>(1., 0);
        qMatrix[1+i][1+i] = std::complex<double>(0, -1);
snuverink_j's avatar
snuverink_j committed
172

173 174 175 176
        invqMatrix[  i][  i] = std::complex<double>(1., 0);
        invqMatrix[  i][1+i] = std::complex<double>(1., 0);
        invqMatrix[1+i][  i] = std::complex<double>(0., -1.);
        invqMatrix[1+i][1+i] = std::complex<double>(0, 1.);
177
    }
178 179
    qMatrix    /= std::sqrt(2);
    invqMatrix /= std::sqrt(2);
snuverink_j's avatar
snuverink_j committed
180 181


182 183
    N = eigenVecM*qMatrix;
    invN =  invqMatrix* invEigenVecM;
snuverink_j's avatar
snuverink_j committed
184 185 186
}


187
MapAnalyser::fMatrix_t MapAnalyser::createRotMatrix_m(std::array<double, 3> phi){
snuverink_j's avatar
snuverink_j committed
188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212
    fMatrix_t R;

    for (int i = 0; i < 3; i++){
        R[2*i][2*i] = std::cos(phi[1]);
        R[2*i+1][2*i+1] = R[2*i][2*i];
        R[2*i][2*i+1] = std::sin(phi[1]);
        R[2*i+1][2*i] = -R[2*i][2*i+1];
    }
    return R;
}


MapAnalyser::fMatrix_t MapAnalyser::createSkewMatrix_m(){
    fMatrix_t S;

    for (int i = 0; i < 3; i++){
        S[2*i][2*i+1] = 1;
        S[2*i+1][2*i] = -1;
    }
    return S;
}


MapAnalyser::fMatrix_t MapAnalyser::realPartOfMatrix_m(cfMatrix_t cM){
    fMatrix_t M;
213 214 215
    for (int i = 0; i < 6; i++){
        for (int j = 0; j < 6; j++){
            M[i][j] = cM[i][j].real();
snuverink_j's avatar
snuverink_j committed
216
        }
217
    }
snuverink_j's avatar
snuverink_j committed
218 219 220 221 222
    return M;
}

MapAnalyser::fMatrix_t MapAnalyser::imagPartOfMatrix_m(cfMatrix_t cM){
    fMatrix_t M;
223 224 225
    for (int i = 0; i < 6; i++){
        for (int j = 0; j < 6; j++){
            M[i][j] = cM[i][j].imag();
snuverink_j's avatar
snuverink_j committed
226
        }
227
    }
snuverink_j's avatar
snuverink_j committed
228 229 230 231 232
    return M;
}

MapAnalyser::cfMatrix_t MapAnalyser::complexTypeMatrix_m(fMatrix_t M){
    cfMatrix_t cM;
233 234 235
    for (int i = 0; i < 6; i++){
        for (int j = 0; j < 6; j++){
            cM[i][j] = std::complex<double>(M[i][j], 0);
snuverink_j's avatar
snuverink_j committed
236
        }
237
    }
snuverink_j's avatar
snuverink_j committed
238 239 240 241
    return cM;
}

MapAnalyser::cfMatrix_t MapAnalyser::invertMatrix_m(const cfMatrix_t& M){
242

snuverink_j's avatar
snuverink_j committed
243 244 245 246 247 248 249 250 251 252 253 254
    gsl_set_error_handler_off();
    //gsl_vector_complex *m = gsl_vector_complex_alloc(6);
    gsl_matrix_complex *m = gsl_matrix_complex_alloc(6, 6);
    gsl_matrix_complex *invm = gsl_matrix_complex_alloc(6, 6);
    gsl_permutation * p = gsl_permutation_alloc(6);
    gsl_complex temp;
    int s;


    //Create invEigenVecMatrix
    for (int i = 0; i < 6; ++i) {
        for (int j = 0; j < 6; ++j) {
255
            GSL_SET_COMPLEX(&temp, std::real(M[i][j]), std::imag(M[i][j]));
snuverink_j's avatar
snuverink_j committed
256 257 258 259 260 261 262 263 264 265
            gsl_matrix_complex_set( m, i, j, temp);
        }
    }

    //invert Eigenvectormatrix
    int eigenDecompStatus = gsl_linalg_complex_LU_decomp(m, p, &s);
    if (eigenDecompStatus != 0){
        std::cout<< "This actually works" << std::endl;
        //gsl_set_error_handler (NULL);

266 267
    }

snuverink_j's avatar
snuverink_j committed
268
    int invertStatus = gsl_linalg_complex_LU_invert(m, p, invm);
269

snuverink_j's avatar
snuverink_j committed
270 271 272 273
    if ( invertStatus ) {
        std::cout << "Error" << std::endl;
        std::exit(1);
    }
274 275


snuverink_j's avatar
snuverink_j committed
276 277 278 279
    if (invertStatus != 0){
        std::cout<< "This actually works" << std::endl;
        //gsl_set_error_handler (NULL);

280
    }
snuverink_j's avatar
snuverink_j committed
281 282 283 284 285 286

    cfMatrix_t invM;
    //Create invEigenVecMatrix
    for (int i = 0; i < 6; ++i) {
        for (int j = 0; j < 6; ++j) {
            invM[i][j] = std::complex<double>(
287 288
                                              GSL_REAL(gsl_matrix_complex_get(invm, i, j)),
                                              GSL_IMAG(gsl_matrix_complex_get(invm, i, j)));
snuverink_j's avatar
snuverink_j committed
289 290 291 292 293 294 295 296 297 298 299 300
        }
    }

    //free space
    gsl_matrix_complex_free(m);
    gsl_matrix_complex_free(invm);
    gsl_permutation_free(p);


    return invM;
}

301 302 303 304
void MapAnalyser::normalizeEigen_m(cfMatrix_t& eigenVecM, cfMatrix_t& invEigenVecM) {
    //normalize eigen Vectors
    for (int i = 0; i < 6; i++){
        double temp = 0;
snuverink_j's avatar
snuverink_j committed
305

306 307 308 309
        for (int j = 0; j < 6; j += 2){
            temp += 2 * (eigenVecM[j][i] * std::conj(eigenVecM[j+1][i])).imag();
        }
        temp = std::fabs(temp);
snuverink_j's avatar
snuverink_j committed
310

311 312 313 314 315 316
        if (temp > 1e-10){
            for (int j = 0; j < 6; j++){
                eigenVecM[j][i] /= std::sqrt(temp);
                invEigenVecM[j][i] /= std::sqrt(temp);
            }
        }
317
    }
318
}