SigmaGenerator.h

Go to the documentation of this file.

//
// Class: SigmaGenerator
// The SigmaGenerator class uses the class <b>ClosedOrbitFinder</b> to get the parameters(inverse bending radius, path length
// field index and tunes) to initialize the sigma matrix.
// The main function of this class is <b>match(double, unsigned int)</b>, where it iteratively tries to find a matched
// distribution for given
// emittances, energy and current. The computation stops when the L2-norm is smaller than a user-defined tolerance. \n
// In default mode it prints all space charge maps, cyclotron maps and second moment matrices. The orbit properties, i.e.
// tunes, average radius, orbit radius, inverse bending radius, path length, field index and frequency error, are printed
// as well.
//
// Copyright (c) 2014, 2018, Matthias Frey, Cristopher Cortes, ETH Zürich
// All rights reserved
//
// Implemented as part of the Semester thesis by Matthias Frey
// "Matched Distributions in Cyclotrons"
//
// Some adaptations done as part of the Bachelor thesis by Cristopher Cortes
// "Limitations of a linear transfer map method for finding matched distributions in high intensity cyclotrons"
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//

#ifndef SIGMAGENERATOR_H
#define SIGMAGENERATOR_H

#include <array>
#include <fstream>
#include <functional>
#include <string>
#include <vector>

#include <boost/numeric/ublas/matrix.hpp>

#include "FixedAlgebra/FTps.h"
#include "Physics/Physics.h"
#include "Physics/Units.h"

#include "Distribution/RealDiracMatrix.h"

class Cyclotron;

class SigmaGenerator
{
public:
    typedef RealDiracMatrix::matrix_t matrix_t;
    typedef RealDiracMatrix::sparse_matrix_t sparse_matrix_t;
    typedef RealDiracMatrix::vector_t vector_t;
    typedef std::vector<double> container_t;
    typedef FTps<double,2*3> Series;
    typedef FVps<double,2*3> Map;
    typedef std::function<Series(double,double,double)> Hamiltonian;
    typedef std::function<Series(double,double,double)> SpaceCharge;


    SigmaGenerator(double I, double ex, double ey, double ez,
                   double E, double m, double q, const Cyclotron* cycl,
                   unsigned int N, unsigned int Nsectors, unsigned int truncOrder, bool write = true);


    bool match(double accuracy, unsigned int maxit, unsigned int maxitOrbit,
               Cyclotron* cycl, double denergy, double rguess, bool full);


    vector_t decouple(const matrix_t& Mturn, sparse_matrix_t& R, sparse_matrix_t& invR);


    double isEigenEllipse(const matrix_t& Mturn, const matrix_t& sigma);

    matrix_t& getSigma();

    unsigned int getIterations() const;

    double getError() const;

    std::array<double,3> getEmittances() const;

    const double& getInjectionRadius() const {
        return rinit_m;
    }

    const double& getInjectionMomentum() const {
        return prinit_m;
    }

private:
    double I_m;
    std::array<double,3> emittance_m;
    double wo_m;
    double E_m;
    double gamma_m;
    double gamma2_m;
    double nh_m;
    double beta_m;
    double mass_m;
    double q_m;
    unsigned int niterations_m;
    bool converged_m;
    double Emin_m;
    double Emax_m;
    unsigned int N_m;
    unsigned int nSectors_m;
    unsigned int nStepsPerSector_m;

    unsigned int nSteps_m;

    double error_m;

    unsigned int truncOrder_m;

    bool write_m;

    matrix_t sigma_m;

    std::vector<matrix_t> sigmas_m;

    Hamiltonian H_m;

    SpaceCharge Hsc_m;

    Series x_m, px_m, z_m, pz_m, l_m, delta_m;

    double rinit_m;
    double prinit_m;

    void initialize(double, double);


    matrix_t updateInitialSigma(const matrix_t&,
                                   const vector_t&,
                                   sparse_matrix_t&,
                                   sparse_matrix_t&);


    void updateSigma(const std::vector<matrix_t>&,
                     const std::vector<matrix_t>&);


    double L2ErrorNorm(const matrix_t&, const matrix_t&);


    std::string float2string(double val);


    void writeOrbitOutput_m(const container_t& r_turn,
                            const container_t& peo,
                            const container_t& h_turn,
                            const container_t& fidx_turn,
                            const container_t& ds_turn);

    void writeMatrix(std::ofstream&, const matrix_t&);

    RealDiracMatrix rdm_m;


    template<class matrix>
    void reduce(matrix&);

    template<class matrix>
    void expand(matrix&);
};


inline
typename SigmaGenerator::matrix_t&
SigmaGenerator::getSigma()
{
    return sigma_m;
}


inline
unsigned int SigmaGenerator::getIterations() const
{
    return (converged_m) ? niterations_m : 0;
}


inline
double SigmaGenerator::getError() const
{
    return error_m;
}


inline
std::array<double,3> SigmaGenerator::getEmittances() const
{
    double bgam = gamma_m*beta_m;
    return std::array<double,3>{{
            emittance_m[0] / Physics::pi / bgam * Units::m2mm * Units::rad2mrad,
        emittance_m[1] / Physics::pi / bgam * Units::m2mm * Units::rad2mrad,
        emittance_m[2] / Physics::pi / bgam * Units::m2mm * Units::rad2mrad
    }};
}


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_t
    matrix_t 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