MapAnalyser.cpp 15.7 KB
Newer Older
snuverink_j's avatar
snuverink_j committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: MapAnalyser
//   Organizes the function for a linear map analysis from
//   ThickTracker.
//   Transfer map -> tunes, symplecticity and stability
//   Sigma Matrix -> (not projected) beam emittance
//
// ------------------------------------------------------------------------
//
// $Author: ganz $
//
// ------------------------------------------------------------------------


#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"))
{ }

//Analyzes a TransferMap for the tunes, symplecticity and stability
void MapAnalyser::linTAnalyze(const fMatrix_t& tMatrix){
    std::ofstream tplot;
    tplot.open ("tunePlot.txt",std::ios::app);
    tplot << std::setprecision(16);
40
    //================================
snuverink_j's avatar
snuverink_j committed
41 42 43 44 45 46 47 48

    IpplTimings::startTimer(mapAnalysis_m);


    fMatrix_t blocktMatrix, scalingMatrix, invScalingMatrix;

    //get the EigenValues/-Vectors
    cfMatrix_t eigenValM, eigenVecM, invEigenVecM;
49 50 51 52 53
    eigenDecomp_m(tMatrix, eigenValM, eigenVecM, invEigenVecM);

    //normalize eigen Vectors
    for (int i=0; i <6 ; i++){
        double temp=0;
snuverink_j's avatar
snuverink_j committed
54 55 56 57 58 59 60 61 62 63 64 65 66 67

        for (int j=0; j <6 ; j+=2){
            temp += 2*(eigenVecM[j][i] * std::conj(eigenVecM[j+1][i])).imag();
        }
        temp=std::fabs(temp);

        if (temp > 1e-10){
            //TODO is it neccessary to normalize the eigenValues?
            //eigenValM[i][i] *= std::sqrt(temp);
            for (int j=0; j <6 ; j++){
                eigenVecM[j][i] /= std::sqrt(temp);
                invEigenVecM[j][i] /= std::sqrt(temp);
            }
        }
68 69 70 71 72 73 74 75 76
    }

    //TODO pack transposition in a neat function
    cfMatrix_t eigenVecMT;
    for (int i=0; i <6 ; i++){
        for (int j=0; j <6 ; j++){
            eigenVecMT[i][j]= eigenVecM[j][i];
        }
    }
snuverink_j's avatar
snuverink_j committed
77

78
    double prec =0.01;
snuverink_j's avatar
snuverink_j committed
79

80 81
    int idx=0;
    std::array<int,6> pairs;
snuverink_j's avatar
snuverink_j committed
82

83 84
    fMatrix_t skewMatrix = createSkewMatrix_m();
    cfMatrix_t cSkewMatrix = complexTypeMatrix_m(skewMatrix);
snuverink_j's avatar
snuverink_j committed
85

86
    fMatrix_t K = imagPartOfMatrix_m(eigenVecMT *cSkewMatrix* eigenVecM);
snuverink_j's avatar
snuverink_j committed
87

88
    for (int i=0; i<6;i++) pairs[i]=-1;
snuverink_j's avatar
snuverink_j committed
89

90 91
    //Printer
    for (int i=0; i <6 ; i++){
snuverink_j's avatar
snuverink_j committed
92 93 94 95 96 97 98 99 100
        for (int j=0; j <6 ; j++){
            if (K[i][j]>1-prec && K[i][j]<1+prec){
                pairs[idx]=i;
                pairs[idx+1]=j;
                idx+=2;

                if (idx==6) break;
            }
        }
101 102
    }

snuverink_j's avatar
snuverink_j committed
103

104 105
    //Fill empty elements in pairs
    std::array<int,6>::iterator pairsIt;
snuverink_j's avatar
snuverink_j committed
106

107
    //Compare eigenvalues
snuverink_j's avatar
snuverink_j committed
108 109 110 111
    int negidx=6-1;
    for (int i=0; i<6; i++) {
        pairsIt = std::find (pairs.begin(), pairs.end(), i);
        if (pairsIt != pairs.end()){
112
            continue;
snuverink_j's avatar
snuverink_j committed
113 114 115 116 117 118 119 120 121 122 123
        }

        for (int j=0; i!=j && j <6 ; j++){
            double diff = std::abs( eigenValM[i][i] - eigenValM[j][j]);
            if (diff< 0.001){
                pairs[negidx]=i;
                pairs[negidx-1]=j;
                negidx-=2;
            }
        }
    }
124

snuverink_j's avatar
snuverink_j committed
125 126 127 128 129 130 131 132 133 134 135
    //Fill the paris vector
    for (int i=0; i < 6 && idx < 6; i++){
        pairsIt = std::find (pairs.begin(), pairs.end(), i);

        if (pairsIt != pairs.end()){
            continue;
        } else{
            pairs[idx]=i;
            idx++;
        }
    }
136

137
    // :FIXME: why commented out? To be removed?
138 139 140
    cfMatrix_t tempM = eigenVecM ;
    //cfMatrix_t tempInvM = eigenVecM ;
    cfMatrix_t tempValM = eigenValM ;
snuverink_j's avatar
snuverink_j committed
141 142


143 144
    for (int i=0; i <6 ; i++){
        eigenValM[i][i]=tempValM[pairs[i]][pairs[i]];
snuverink_j's avatar
snuverink_j committed
145 146 147 148 149
        for (int j=0; j <6 ; j++){
            eigenVecM[j][i]=tempM[j][pairs[i]];
            //eigenVecMT[i][j]=eigenVecM[j][i];
            //invEigenVecM[i][j]=tempInvM[pairs[i]][j];
        }
150 151
    }

snuverink_j's avatar
snuverink_j committed
152
    invEigenVecM= invertMatrix_m(eigenVecM);
153
    cfMatrix_t cblocktMatrix = getBlockDiagonal_m(tMatrix, eigenVecM, invEigenVecM);
snuverink_j's avatar
snuverink_j committed
154

155 156
    FVector<std::complex<double>, 3> betaTunes, betaTunes2;
    FVector<double, 3> betaTunes3;
snuverink_j's avatar
snuverink_j committed
157

158
    // :FIXME: why commented out
159
    //rearrangeEigen(eigenValM, eigenVecM);
snuverink_j's avatar
snuverink_j committed
160

161 162 163
    for (int i = 0; i < 3; i++){
        betaTunes[i]=std::log(eigenValM[i*2][i*2])/ (2*Physics::pi * std::complex<double>(0, 1));
        betaTunes[i]= std::asin(cblocktMatrix[i*2][i*2+1])/(std::complex<double>(2*Physics::pi, 0));
snuverink_j's avatar
snuverink_j committed
164

165 166
        betaTunes2[i]= std::acos(cblocktMatrix[i*2][i*2])/(std::complex<double>(2*Physics::pi, 0));
        betaTunes2[i]= std::acos(cblocktMatrix[i*2][i*2].real())/(std::complex<double>(2*Physics::pi, 0));
snuverink_j's avatar
snuverink_j committed
167

168
        betaTunes3[i]= std::acos(eigenValM[i*2][i*2].real()) / (2*Physics::pi);
snuverink_j's avatar
snuverink_j committed
169

170 171 172
        double lenEigenVec = 0;
        for (int j = 0; j < 3; j++){
            lenEigenVec += std::abs(eigenVecM[i][j]);
snuverink_j's avatar
snuverink_j committed
173

174 175 176 177 178 179 180 181 182
        }
        lenEigenVec = std::sqrt(lenEigenVec);
        tplot<<"1: "<<betaTunes[i] <<std::endl;
        tplot<<"2: "<<betaTunes2[i]<< std::endl;
        //tplot<<"3: "<<lenEigenVec<< std::endl;
        tplot<<"3: ("<<betaTunes3[i]<< ",0)"<< std::endl;

        IpplTimings::stopTimer(mapAnalysis_m);
        //TODO: do something with the tunes etc
183
        //:FIXME: are there some plans what to do?
184
    }
snuverink_j's avatar
snuverink_j committed
185 186 187 188 189 190 191

}


//Analyses the Sigma Matrix
void MapAnalyser::linSigAnalyze(fMatrix_t& sigMatrix){
    IpplTimings::startTimer(bunchAnalysis_m);
192
    fMatrix_t sigmaBlockM;
snuverink_j's avatar
snuverink_j committed
193

194 195 196 197 198
    fMatrix_t skewMatrix;
    for (int i=0; i <6 ; i=i+2){
        skewMatrix[0+i][1+i]=1.;
        skewMatrix[1+i][0+i]=-1.;
    }
snuverink_j's avatar
snuverink_j committed
199

200
    sigMatrix=sigMatrix*skewMatrix;
snuverink_j's avatar
snuverink_j committed
201

202
    cfMatrix_t eigenValM, eigenVecM, invEigenVecM;
snuverink_j's avatar
snuverink_j committed
203

204
    eigenDecomp_m(sigMatrix,eigenValM, eigenVecM, invEigenVecM);
snuverink_j's avatar
snuverink_j committed
205

206
    cfMatrix_t cblocksigM = getBlockDiagonal_m(sigMatrix, eigenVecM,invEigenVecM);
snuverink_j's avatar
snuverink_j committed
207 208


209 210 211 212
    for (int i=0; i<6; i++){
        for (int j =0; j<6; j++){
            sigmaBlockM[i][j]=cblocksigM[i][j].real();
        }
snuverink_j's avatar
snuverink_j committed
213

214
    }
snuverink_j's avatar
snuverink_j committed
215

216 217
    cfMatrix_t teigenValM, teigenVecM, tinvEigenVecM;
    eigenDecomp_m(sigmaBlockM, teigenValM, teigenVecM, tinvEigenVecM);
snuverink_j's avatar
snuverink_j committed
218

219
    for (int i=0; i<6; i++){
snuverink_j's avatar
snuverink_j committed
220

221
    }
snuverink_j's avatar
snuverink_j committed
222 223


224 225
    //TODO: Where to go with EigenValues?
    IpplTimings::stopTimer(bunchAnalysis_m);
snuverink_j's avatar
snuverink_j committed
226 227 228 229 230 231
}


//Returns the eigenDecomositon for the fMatrix M = eigenVec * eigenVal * invEigenVec
void MapAnalyser::eigenDecomp_m(const fMatrix_t& M, cfMatrix_t& eigenVal, cfMatrix_t& eigenVec, cfMatrix_t& invEigenVec){

232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282
    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
283 284 285 286 287 288 289


}


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

292
    cfMatrix_t cM, qMatrix, invqMatrix, nMatrix, invnMatrix, rMatrix;
snuverink_j's avatar
snuverink_j committed
293

294 295 296 297 298
    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
299

300 301 302 303 304
    for (int i=0; i <6 ; i=i+2){
        qMatrix[0+i][0+i]=std::complex<double>(1.,0);
        qMatrix[0+i][1+i]=std::complex<double>(0,1.);
        qMatrix[1+i][0+i]=std::complex<double>(1.,0);
        qMatrix[1+i][1+i]=std::complex<double>(0,-1);
snuverink_j's avatar
snuverink_j committed
305

306 307 308 309 310 311 312
        invqMatrix[0+i][0+i]=std::complex<double>(1.,0);
        invqMatrix[0+i][1+i]=std::complex<double>(1.,0);
        invqMatrix[1+i][0+i]=std::complex<double>(0.,-1.);
        invqMatrix[1+i][1+i]=std::complex<double>(0,1.);
    }
    qMatrix/=std::sqrt(2.);
    invqMatrix/=std::sqrt(2);
snuverink_j's avatar
snuverink_j committed
313

314 315
    nMatrix=eigenVecM*qMatrix;
    invnMatrix= invqMatrix* invEigenVecM;
snuverink_j's avatar
snuverink_j committed
316 317


318
    rMatrix= invnMatrix * cM * nMatrix;
snuverink_j's avatar
snuverink_j committed
319

320
    return rMatrix;
snuverink_j's avatar
snuverink_j committed
321 322 323 324 325 326 327 328 329

}


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;

    for (int i=0; i<6; i++){
330 331 332 333
        for (int j =0; j<6; j++){
            ctM[i][j]=std::complex<double>(tM[i][j],0);
        }
    }
snuverink_j's avatar
snuverink_j committed
334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366

    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());
    }


    R1=createRotMatrix_m(phi);

    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();
        }
    }
}


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;

367 368 369
    //std::ofstream tmap;
    //tmap.open ("TransferMap.txt",std::ios::app);
    //tmap << std::setprecision(16);
snuverink_j's avatar
snuverink_j committed
370 371


372 373 374
    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
375
        }
376
    }
snuverink_j's avatar
snuverink_j committed
377

378 379 380 381 382
    for (int i=0; i <6 ; i=i+2){
        qMatrix[0+i][0+i]=std::complex<double>(1.,0);
        qMatrix[0+i][1+i]=std::complex<double>(0,1.);
        qMatrix[1+i][0+i]=std::complex<double>(1.,0);
        qMatrix[1+i][1+i]=std::complex<double>(0,-1);
snuverink_j's avatar
snuverink_j committed
383

384 385 386 387 388 389 390
        invqMatrix[0+i][0+i]=std::complex<double>(1.,0);
        invqMatrix[0+i][1+i]=std::complex<double>(1.,0);
        invqMatrix[1+i][0+i]=std::complex<double>(0.,-1.);
        invqMatrix[1+i][1+i]=std::complex<double>(0,1.);
    }
    qMatrix/=std::sqrt(2.);
    invqMatrix/=std::sqrt(2);
snuverink_j's avatar
snuverink_j committed
391 392


393 394
    N=eigenVecM*qMatrix;
    invN= invqMatrix* invEigenVecM;
snuverink_j's avatar
snuverink_j committed
395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424
}


MapAnalyser::fMatrix_t MapAnalyser::createRotMatrix_m(std::array<double,3> phi){
    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;
    for (int i=0; i<6; i++){
425 426
        for (int j =0; j<6; j++){
            M[i][j]=cM[i][j].real();
snuverink_j's avatar
snuverink_j committed
427
        }
428
    }
snuverink_j's avatar
snuverink_j committed
429 430 431 432 433 434
    return M;
}

MapAnalyser::fMatrix_t MapAnalyser::imagPartOfMatrix_m(cfMatrix_t cM){
    fMatrix_t M;
    for (int i=0; i<6; i++){
435 436
        for (int j =0; j<6; j++){
            M[i][j]=cM[i][j].imag();
snuverink_j's avatar
snuverink_j committed
437
        }
438
    }
snuverink_j's avatar
snuverink_j committed
439 440 441 442 443 444
    return M;
}

MapAnalyser::cfMatrix_t MapAnalyser::complexTypeMatrix_m(fMatrix_t M){
    cfMatrix_t cM;
    for (int i=0; i<6; i++){
445 446
        for (int j =0; j<6; j++){
            cM[i][j]=std::complex<double>(M[i][j],0);
snuverink_j's avatar
snuverink_j committed
447
        }
448
    }
snuverink_j's avatar
snuverink_j committed
449 450 451 452 453
    return cM;
}


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

snuverink_j's avatar
snuverink_j committed
455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477
    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) {
            GSL_SET_COMPLEX(&temp,std::real(M[i][j]),std::imag(M[i][j]));
            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);

478 479
    }

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

snuverink_j's avatar
snuverink_j committed
482 483 484 485
    if ( invertStatus ) {
        std::cout << "Error" << std::endl;
        std::exit(1);
    }
486 487


snuverink_j's avatar
snuverink_j committed
488 489 490 491
    if (invertStatus != 0){
        std::cout<< "This actually works" << std::endl;
        //gsl_set_error_handler (NULL);

492
    }
snuverink_j's avatar
snuverink_j committed
493 494 495 496 497 498

    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>(
499 500
                                              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
501 502 503 504 505 506 507 508 509 510 511 512 513
        }
    }

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


    return invM;

}

514
#if 0
snuverink_j's avatar
snuverink_j committed
515 516
//TODO Work in progress
void MapAnalyser::rearrangeEigen_m(cfMatrix_t& eigenVal, cfMatrix_t& EigenVec){
517
#ifdef PHIL_WRITE
snuverink_j's avatar
snuverink_j committed
518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557
    std::ofstream out;
    out.open("OUT.txt", std::ios::app);
#endif
    double precision = 1e-9;
    int pair = 0;
    std::complex<double> mult;
    std::pair<int, int> pairs[3];
    cfMatrix_t eigVa, eigVe, temp;

    out << "+++++++++++++++++++++++++++++++" << std::endl;
    //out << EigenVec << std::endl;
    for (int i = 0; pair < 3 && i <2*DIM; i++){
        for (int j = 0; pair < 3 && j < i; j++){
            mult=0;
            if (true){//std::abs(EigenVal[i][i].real() - EigenVal[j][j].real()) < precision
                    //&& std::abs(EigenVal[i][i].imag() + EigenVal[j][j].imag() ) < precision ){
                out << eigenVal[i][i] << "    " << eigenVal[j][j] << std::endl;
                out << "/////////////Calc///////////" << std::endl;
                for (int k = 0; k <DIM; k++){

                    out << -1. * EigenVec[2*k+1][i]  << EigenVec[2*k][j] << std::endl;
                    out << EigenVec[2*k][i] << EigenVec[2*k+1][j] << std::endl;


                    mult += -1. * EigenVec[2*k+1][i] * EigenVec[2*k][j];
                    mult += EigenVec[2*k][i] * EigenVec[2*k+1][j];
                    out <<"in between mult:  " <<  mult << std::endl;
                }

                out << mult << std::endl;
                out << "Criteria    " << std::abs(mult.imag()-1.) << std::endl;
                if (std::abs(mult.imag()-1.) < precision){
                    pairs[pair]= std::pair<int,int>(i,j);
                    out << i << "    " << j << std::endl;
                    pair++;
                }

            }

         }
558
    }
snuverink_j's avatar
snuverink_j committed
559
}
560
#endif