SigmaGenerator.h 48 KB
Newer Older
1 2 3
/**
 * @file SigmaGenerator.h
 * The SigmaGenerator class uses the class <b>ClosedOrbitFinder</b> to get the parameters (inverse bending radius, path length
frey_m's avatar
frey_m committed
4
 * field index and tunes) to initialize the sigma matrix.
5 6
 * The main function of this class is <b>match(value_type, size_type)</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
7 8 9 10
 * 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.
 *
frey_m's avatar
frey_m committed
11 12
 * @author Matthias Frey, Cristopher Cortes
 * @version 1.1
13
 */
14 15 16 17 18
#ifndef SIGMAGENERATOR_H
#define SIGMAGENERATOR_H

#include <array>
#include <cmath>
19
#include <fstream>
adelmann's avatar
adelmann committed
20
#include <functional>
21
#include <iomanip>
22 23 24 25
#include <iterator>
#include <limits>
#include <list>
#include <numeric>
26
#include <sstream>
adelmann's avatar
adelmann committed
27
#include <string>
28
#include <utility>
29 30
#include <vector>

31
#include "Physics/Physics.h"
32
#include "Utilities/Options.h"
Andreas Adelmann's avatar
Andreas Adelmann committed
33 34
#include "Utilities/Options.h"
#include "Utilities/OpalException.h"
35

36 37 38 39
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/vector.hpp>

40 41 42
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
frey_m's avatar
frey_m committed
43
#include <boost/numeric/ublas/io.hpp>
44 45 46 47 48

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

49
#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
50
#if BOOST_VERSION >= 106000
gsell's avatar
gsell committed
51
#include <boost/numeric/odeint/integrate/check_adapter.hpp>
52
#endif
53 54 55 56

#include "matrix_vector_operation.h"
#include "ClosedOrbitFinder.h"
#include "MapGenerator.h"
adelmann's avatar
adelmann committed
57
#include "Harmonics.h"
58

59

frey_m's avatar
frey_m committed
60 61
extern Inform *gmsg;

62
/// @brief This class computes the matched distribution
63 64 65
template<typename Value_type, typename Size_type>
class SigmaGenerator
{
frey_m's avatar
frey_m committed
66 67 68 69 70 71 72 73 74
public:
    /// Type of variables
    typedef Value_type value_type;
    /// Type of constant variables
    typedef const value_type const_value_type;
    /// Type for specifying sizes
    typedef Size_type size_type;
    /// Type for storing maps
    typedef boost::numeric::ublas::matrix<value_type> matrix_type;
kraus's avatar
kraus committed
75

76
    typedef std::complex<value_type> complex_t;
kraus's avatar
kraus committed
77

78 79
    /// Type for storing complex matrices
    typedef boost::numeric::ublas::matrix<complex_t> cmatrix_type;
kraus's avatar
kraus committed
80

frey_m's avatar
frey_m committed
81
    /// Type for storing the sparse maps
82
    typedef boost::numeric::ublas::compressed_matrix<complex_t,
frey_m's avatar
frey_m committed
83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
                                                     boost::numeric::ublas::row_major
                                                     > sparse_matrix_type;
    /// Type for storing vectors
    typedef boost::numeric::ublas::vector<value_type> vector_type;
    /// Container for storing the properties for each angle
    typedef std::vector<value_type> container_type;
    /// Type of the truncated powere series
    typedef FTps<value_type,2*3> Series;
    /// Type of a map
    typedef FVps<value_type,2*3> Map;
    /// Type of the Hamiltonian for the cyclotron
    typedef std::function<Series(value_type,value_type,value_type)> Hamiltonian;
    /// Type of the Hamiltonian for the space charge
    typedef std::function<Series(value_type,value_type,value_type)> SpaceCharge;

    /// Constructs an object of type SigmaGenerator
    /*!
     * @param I specifies the current for which a matched distribution should be found, \f$ [I] = A \f$
     * @param ex is the emittance in x-direction (horizontal), \f$ \left[\varepsilon_{x}\right] = \pi\ mm\ mrad  \f$
     * @param ey is the emittance in y-direction (longitudinal), \f$ \left[\varepsilon_{y}\right] = \pi\ mm\ mrad  \f$
     * @param ez is the emittance in z-direction (vertical), \f$ \left[\varepsilon_{z}\right] = \pi\ mm\ mrad  \f$
     * @param E is the energy, \f$ \left[E\right] = MeV \f$
     * @param m is the mass of the particles \f$ \left[m\right] = \frac{MeV}{c^{2}} \f$
106
     * @param cycl is the cyclotron element
frey_m's avatar
frey_m committed
107 108
     * @param N is the number of integration steps (closed orbit computation). That's why its also the number
     *    of maps (for each integration step a map)
109
     * @param Nsectors is the number of sectors that the field map is averaged over (closed orbit computation)
frey_m's avatar
frey_m committed
110 111 112 113
     * @param truncOrder is the truncation order for power series of the Hamiltonian
     * @param write is a boolean (default: true). If true all maps of all iterations are stored, otherwise not.
     */
    SigmaGenerator(value_type I, value_type ex, value_type ey, value_type ez,
114
                   value_type E, value_type m, const Cyclotron* cycl,
115
                   size_type N, size_type Nsectors, size_type truncOrder, bool write = true);
frey_m's avatar
frey_m committed
116 117 118 119 120 121 122 123 124

    /// Searches for a matched distribution.
    /*!
     * Returns "true" if a matched distribution within the accuracy could be found, returns "false" otherwise.
     * @param accuracy is used for computing the equilibrium orbit (to a certain accuracy)
     * @param maxit is the maximum number of iterations (the program stops if within this value no stationary
     *    distribution was found)
     * @param maxitOrbit is the maximum number of iterations for finding closed orbit
     * @param angle defines the start of the sector (one can choose any angle between 0° and 360°)
frey_m's avatar
frey_m committed
125
     * @param denergy energy step size for closed orbit finder [MeV]
126
     * @param rguess value of radius for closed orbit finder
frey_m's avatar
frey_m committed
127 128 129 130
     * @param type specifies the magnetic field format (e.g. CARBONCYCL)
     * @param harmonic is a boolean. If "true" the harmonics are used instead of the closed orbit finder.
     * @param full match over full turn not just single sector
     */
131
    bool match(value_type accuracy, size_type maxit, size_type maxitOrbit,
frey_m's avatar
frey_m committed
132
               Cyclotron* cycl, value_type denergy, value_type rguess, bool harmonic, bool full);
kraus's avatar
kraus committed
133

134 135 136 137
    /*!
     * Eigenvalue / eigenvector solver
     * @param Mturn is a 6x6 dimensional one turn transfer matrix
     * @param R is the 6x6 dimensional transformation matrix (gets computed)
frey_m's avatar
frey_m committed
138
     */
139
    void eigsolve_m(const matrix_type& Mturn, sparse_matrix_type& R);
kraus's avatar
kraus committed
140

frey_m's avatar
frey_m committed
141
    /*!
142 143 144 145 146 147 148 149 150
     * @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);

    /// 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>.
kraus's avatar
kraus committed
151
     *
frey_m's avatar
frey_m committed
152
     * @param Mturn is a 6x6 dimensional one turn transfer matrix
153 154
     * @param R is the 6x6 dimensional transformation matrix (gets computed)
     * @param invR is the 6x6 dimensional inverse transformation (gets computed)
frey_m's avatar
frey_m committed
155
     */
156
    void decouple(const matrix_type& Mturn, sparse_matrix_type& R, sparse_matrix_type& invR);
frey_m's avatar
frey_m committed
157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176

    /// Checks if the sigma-matrix is an eigenellipse and returns the L2 error.
    /*!
     * The idea of this function is taken from Dr. Christian Baumgarten's program.
     * @param Mturn is the one turn transfer matrix
     * @param sigma is the sigma matrix to be tested
     */
    value_type isEigenEllipse(const matrix_type& Mturn, const matrix_type& sigma);

    /// Returns the converged stationary distribution
    matrix_type& getSigma();

    /// Returns the number of iterations needed for convergence (if not converged, it returns zero)
    size_type getIterations() const;

    /// Returns the error (if the program didn't converged, one can call this function to check the error)
    value_type getError() const;

    /// Returns the emittances (ex,ey,ez) in \f$ \pi\ mm\ mrad \f$ for which the code converged (since the whole simulation is done on normalized emittances)
    std::array<value_type,3> getEmittances() const;
kraus's avatar
kraus committed
177

frey_m's avatar
frey_m committed
178 179 180
    const double& getInjectionRadius() const {
        return rinit_m;
    }
kraus's avatar
kraus committed
181

frey_m's avatar
frey_m committed
182 183 184
    const double& getInjectionMomentum() const {
        return prinit_m;
    }
185 186

    private:
frey_m's avatar
frey_m committed
187 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 213 214
    /// Stores the value of the current, \f$ \left[I\right] = A \f$
    value_type I_m;
    /// Stores the desired emittances, \f$ \left[\varepsilon_{x}\right] = \left[\varepsilon_{y}\right] = \left[\varepsilon_{z}\right] = mm \ mrad \f$
    std::array<value_type,3> emittance_m;
    /// Is the orbital frequency, \f$ \left[\omega_{o}\right] = \frac{1}{s} \f$
    value_type wo_m;
    /// Stores the user-define energy, \f$ \left[E\right] = MeV \f$
    value_type E_m;
    /// Relativistic factor (which can be computed out ot the kinetic energy and rest mass (potential energy)), \f$ \left[\gamma\right] = 1 \f$
    value_type gamma_m;
    /// Relativistic factor squared
    value_type gamma2_m;
    /// Harmonic number, \f$ \left[N_{h}\right] = 1 \f$
    value_type nh_m;
    /// Velocity (c/v), \f$ \left[\beta\right] = 1 \f$
    value_type beta_m;
    /// Is the mass of the particles, \f$ \left[m\right] = \frac{MeV}{c^{2}} \f$
    value_type m_m;
    /// Is the number of iterations needed for convergence
    size_type niterations_m;
    /// Is true if converged, false otherwise
    bool converged_m;
    /// Minimum energy needed in cyclotron, \f$ \left[E_{min}\right] = MeV \f$
    value_type Emin_m;
    /// Maximum energy reached in cyclotron, \f$ \left[E_{max}\right] = MeV \f$
    value_type Emax_m;
    /// Number of integration steps for closed orbit computation
    size_type N_m;
215 216
    /// Number of (symmetric) sectors the field is averaged over
    size_type nSectors_m;
frey_m's avatar
frey_m committed
217 218
    /// Number of integration steps per sector (--> also: number of maps per sector)
    size_type nStepsPerSector_m;
kraus's avatar
kraus committed
219

frey_m's avatar
frey_m committed
220 221
    /// Number of integration steps in total
    size_type nSteps_m;
kraus's avatar
kraus committed
222

frey_m's avatar
frey_m committed
223 224
    /// Error of computation
    value_type error_m;
kraus's avatar
kraus committed
225

frey_m's avatar
frey_m committed
226 227
    /// Truncation order of the power series for the Hamiltonian (cyclotron and space charge)
    size_type truncOrder_m;
kraus's avatar
kraus committed
228

frey_m's avatar
frey_m committed
229 230
    /// Decides for writing output (default: true)
    bool write_m;
kraus's avatar
kraus committed
231

frey_m's avatar
frey_m committed
232 233
    /// Stores the stationary distribution (sigma matrix)
    matrix_type sigma_m;
234

frey_m's avatar
frey_m committed
235 236
    /// Vector storing the second moment matrix for each angle
    std::vector<matrix_type> sigmas_m;
237

frey_m's avatar
frey_m committed
238 239
    /// Stores the Hamiltonian of the cyclotron
    Hamiltonian H_m;
kraus's avatar
kraus committed
240

frey_m's avatar
frey_m committed
241 242
    /// Stores the Hamiltonian for the space charge
    SpaceCharge Hsc_m;
kraus's avatar
kraus committed
243

frey_m's avatar
frey_m committed
244 245
    /// All variables x, px, y, py, z, delta
    Series x_m, px_m, y_m, py_m, z_m, delta_m;
kraus's avatar
kraus committed
246

frey_m's avatar
frey_m committed
247 248 249 250 251 252
    double rinit_m;
    double prinit_m;

    /*! Initializes a first guess of the sigma matrix with the assumption of
     * a spherical symmetric beam (ex = ey = ez). For each angle split the
     * same initial guess is taken.
kraus's avatar
kraus committed
253
     *
frey_m's avatar
frey_m committed
254 255 256 257
     * @param nuz is the vertical tune
     * @param ravg is the average radius of the closed orbit
     */
    void initialize(value_type, value_type);
kraus's avatar
kraus committed
258

frey_m's avatar
frey_m committed
259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282
    /// Computes the new initial sigma matrix
    /*!
     * @param M is the 6x6 one turn transfer matrix
     * @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&);

    /// Computes new sigma matrices (one for each angle)
    /*!
     * Mscs is a vector of all space charge maps
     * Mcycs is a vector of all cyclotron maps
     */
    void updateSigma(const std::vector<matrix_type>&,
                     const std::vector<matrix_type>&);

    /// Returns the L2-error norm between the old and new sigma-matrix
    /*!
     * @param oldS is the old sigma matrix
     * @param newS is the new sigma matrix
     */
    value_type L2ErrorNorm(const matrix_type&, const matrix_type&);
kraus's avatar
kraus committed
283 284


frey_m's avatar
frey_m committed
285 286 287 288 289 290 291 292 293 294 295 296
    /// Returns the Lp-error norm between the old and new sigma-matrix
    /*!
     * @param oldS is the old sigma matrix
     * @param newS is the new sigma matrix
     */
    value_type L1ErrorNorm(const matrix_type&, const matrix_type&);

    /// Transforms a floating point value to a string
    /*!
     * @param val is the floating point value which is transformed to a string
     */
    std::string float2string(value_type val);
kraus's avatar
kraus committed
297

frey_m's avatar
frey_m committed
298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315
    /// Called within SigmaGenerator::match().
    /*!
     * @param tunes
     * @param ravg is the average radius [m]
     * @param r_turn is the radius [m]
     * @param peo is the momentum
     * @param h_turn is the inverse bending radius
     * @param fidx_turn is the field index
     * @param ds_turn is the path length element
     */
    void writeOrbitOutput_m(const std::pair<value_type,value_type>& tunes,
                            const value_type& ravg,
                            const value_type& freqError,
                            const container_type& r_turn,
                            const container_type& peo,
                            const container_type& h_turn,
                            const container_type&  fidx_turn,
                            const container_type&  ds_turn);
316 317 318 319 320 321 322
};

// -----------------------------------------------------------------------------------------------------------------------
// PUBLIC MEMBER FUNCTIONS
// -----------------------------------------------------------------------------------------------------------------------

template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
323 324 325 326 327 328
SigmaGenerator<Value_type, Size_type>::SigmaGenerator(value_type I,
                                                      value_type ex,
                                                      value_type ey,
                                                      value_type ez,
                                                      value_type E,
                                                      value_type m,
329
                                                      const Cyclotron* cycl,
frey_m's avatar
frey_m committed
330
                                                      size_type N,
331
                                                      size_type Nsectors,
frey_m's avatar
frey_m committed
332 333
                                                      size_type truncOrder,
                                                      bool write)
334
    : I_m(I)
335
    , wo_m(cycl->getRfFrequ(0)*1E6/cycl->getCyclHarm()*2.0*Physics::pi)
336 337 338
    , E_m(E)
    , gamma_m(E/m+1.0)
    , gamma2_m(gamma_m*gamma_m)
339
    , nh_m(cycl->getCyclHarm())
340 341 342 343
    , beta_m(std::sqrt(1.0-1.0/gamma2_m))
    , m_m(m)
    , niterations_m(0)
    , converged_m(false)
344 345
    , Emin_m(cycl->getFMLowE())
    , Emax_m(cycl->getFMHighE())
346
    , N_m(N)
347
    , nSectors_m(Nsectors)
348
    , nStepsPerSector_m(N/cycl->getSymmetry())
349 350 351 352 353 354 355
    , nSteps_m(N)
    , error_m(std::numeric_limits<value_type>::max())
    , truncOrder_m(truncOrder)
    , write_m(write)
    , sigmas_m(nStepsPerSector_m)
    , rinit_m(0.0)
    , prinit_m(0.0)
356
{
357 358
    // set emittances (initialization like that due to old compiler version)
    // [ex] = [ey] = [ez] = pi*mm*mrad --> [emittance] = mm mrad
359 360 361
    emittance_m[0] = ex * Physics::pi;
    emittance_m[1] = ey * Physics::pi;
    emittance_m[2] = ez * Physics::pi;
362 363

    // minimum beta*gamma
364
    value_type minGamma = Emin_m / m_m + 1.0;
365 366 367 368 369 370 371
    value_type bgam = std::sqrt(minGamma * minGamma - 1.0);

    // normalized emittance (--> multiply with beta*gamma)
    // [emittance] = mm mrad
    emittance_m[0] *= bgam;
    emittance_m[1] *= bgam;
    emittance_m[2] *= bgam;
kraus's avatar
kraus committed
372

adelmann's avatar
adelmann committed
373 374
    // Define the Hamiltonian
    Series::setGlobalTruncOrder(truncOrder_m);
kraus's avatar
kraus committed
375

adelmann's avatar
adelmann committed
376 377 378 379 380 381 382
    // infinitesimal elements
    x_m = Series::makeVariable(0);
    px_m = Series::makeVariable(1);
    y_m = Series::makeVariable(2);
    py_m = Series::makeVariable(3);
    z_m = Series::makeVariable(4);
    delta_m = Series::makeVariable(5);
kraus's avatar
kraus committed
383

adelmann's avatar
adelmann committed
384 385 386 387 388
    H_m = [&](value_type h, value_type kx, value_type ky) {
        return 0.5*px_m*px_m + 0.5*kx*x_m*x_m - h*x_m*delta_m +
               0.5*py_m*py_m + 0.5*ky*y_m*y_m +
               0.5*delta_m*delta_m/gamma2_m;
    };
kraus's avatar
kraus committed
389

adelmann's avatar
adelmann committed
390 391
    Hsc_m = [&](value_type sigx, value_type sigy, value_type sigz) {
        // convert m from MeV/c^2 to eV*s^{2}/m^{2}
392
        value_type m = m_m * 1.0e6 / (Physics::c * Physics::c);
kraus's avatar
kraus committed
393

adelmann's avatar
adelmann committed
394
        // formula (57)
395 396 397
        value_type lam = 2.0 * Physics::pi*Physics::c / (wo_m * nh_m); // wavelength, [lam] = m
        value_type K3 = 3.0 * /* physics::q0 */ 1.0 * I_m * lam / (20.0 * std::sqrt(5.0) * Physics::pi * Physics::epsilon_0 * m *
                        Physics::c * Physics::c * Physics::c * beta_m * beta_m * gamma_m * gamma2_m);            // [K3] = m
kraus's avatar
kraus committed
398

adelmann's avatar
adelmann committed
399
        value_type milli = 1.0e-3;
kraus's avatar
kraus committed
400

adelmann's avatar
adelmann committed
401
        // formula (30), (31)
snuverink_j's avatar
snuverink_j committed
402
        // [sigma(0,0)] = mm^{2} --> [sx] = [sy] = [sz] = mm
kraus's avatar
kraus committed
403
        // multiply with 0.001 to get meter --> [sx] = [sy] = [sz] = m
adelmann's avatar
adelmann committed
404 405 406
        value_type sx = std::sqrt(std::fabs(sigx)) * milli;
        value_type sy = std::sqrt(std::fabs(sigy)) * milli;
        value_type sz = std::sqrt(std::fabs(sigz)) * milli;
kraus's avatar
kraus committed
407

adelmann's avatar
adelmann committed
408
        value_type tmp = sx * sy;                                           // [tmp] = m^{2}
kraus's avatar
kraus committed
409

adelmann's avatar
adelmann committed
410 411
        value_type f = std::sqrt(tmp) / (3.0 * gamma_m * sz);               // [f] = 1
        value_type kxy = K3 * std::fabs(1.0 - f) / ((sx + sy) * sz); // [kxy] = 1/m
kraus's avatar
kraus committed
412

adelmann's avatar
adelmann committed
413 414 415
        value_type Kx = kxy / sx;
        value_type Ky = kxy / sy;
        value_type Kz = K3 * f / (tmp * sz);
kraus's avatar
kraus committed
416

adelmann's avatar
adelmann committed
417 418 419 420
        return -0.5 * Kx * x_m * x_m
               -0.5 * Ky * y_m * y_m
               -0.5 * Kz * z_m * z_m * gamma2_m;
     };
421 422 423
}

template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
424 425 426
  bool SigmaGenerator<Value_type, Size_type>::match(value_type accuracy,
                                                    size_type maxit,
                                                    size_type maxitOrbit,
427
                                                    Cyclotron* cycl,
frey_m's avatar
frey_m committed
428
                                                    value_type denergy,
frey_m's avatar
frey_m committed
429
                                                    value_type rguess,
430
                                                    bool harmonic, bool full)
431 432
{
    /* compute the equilibrium orbit for energy E_
433 434 435 436 437
     * and get the the following properties:
     * - inverse bending radius h
     * - step sizes of path ds
     * - tune nuz
     */
cortes_c's avatar
cortes_c committed
438

439
    try {
440 441
        if ( !full )
            nSteps_m = nStepsPerSector_m;
Cris_cortes's avatar
Cris_cortes committed
442

443 444 445 446 447 448
        // object for space charge map and cyclotron map
        MapGenerator<value_type,
                     size_type,
                     Series,
                     Map,
                     Hamiltonian,
449
                     SpaceCharge> mapgen(nSteps_m);
450 451

        // compute cyclotron map and space charge map for each angle and store them into a vector
452
        std::vector<matrix_type> Mcycs(nSteps_m), Mscs(nSteps_m);
453

454
        container_type h(nSteps_m), r(nSteps_m), ds(nSteps_m), fidx(nSteps_m);
455 456 457 458 459
        value_type ravg = 0.0, const_ds = 0.0;
        std::pair<value_type,value_type> tunes;

        if (!harmonic) {
            ClosedOrbitFinder<value_type, size_type,
460
                boost::numeric::odeint::runge_kutta4<container_type> > cof(m_m, N_m, cycl, false, nSectors_m);
461

frey_m's avatar
frey_m committed
462
            if ( !cof.findOrbit(accuracy, maxitOrbit, E_m, denergy, rguess) ) {
463 464 465 466
                throw OpalException("SigmaGenerator::match()",
                                    "Closed orbit finder didn't converge.");
            }

467
            cof.computeOrbitProperties(E_m);
468 469

            // properties of one turn
470 471
            tunes = cof.getTunes();
            ravg = cof.getAverageRadius();                   // average radius
kraus's avatar
kraus committed
472

473
            value_type angle = cycl->getPHIinit();
frey_m's avatar
frey_m committed
474
            container_type h_turn = cof.getInverseBendingRadius(angle);
475
            container_type r_turn = cof.getOrbit(angle);
frey_m's avatar
frey_m committed
476 477
            container_type ds_turn = cof.getPathLength(angle);
            container_type fidx_turn = cof.getFieldIndex(angle);
478

479
            container_type peo = cof.getMomentum(angle);
kraus's avatar
kraus committed
480

481 482 483 484 485 486 487
            // write properties to file
            if (write_m)
                writeOrbitOutput_m(tunes, ravg, cof.getFrequencyError(),
                                   r_turn, peo, h_turn, fidx_turn, ds_turn);

            // write to terminal
            *gmsg << "* ----------------------------" << endl
488
                  << "* Closed orbit info:" << endl
489 490 491 492
                  << "*" << endl
                  << "* average radius: " << ravg << " [m]" << endl
                  << "* initial radius: " << r_turn[0] << " [m]" << endl
                  << "* initial momentum: " << peo[0] << " [Beta Gamma]" << endl
493
                  << "* frequency error: " << cof.getFrequencyError()*100 <<" [ % ] "<< endl
494 495 496 497
                  << "* horizontal tune: " << tunes.first << endl
                  << "* vertical tune: " << tunes.second << endl
                  << "* ----------------------------" << endl << endl;

498 499 500 501 502
            // copy properties
            std::copy_n(r_turn.begin(), nSteps_m, r.begin());
            std::copy_n(h_turn.begin(), nSteps_m, h.begin());
            std::copy_n(fidx_turn.begin(), nSteps_m, fidx.begin());
            std::copy_n(ds_turn.begin(), nSteps_m, ds.begin());
kraus's avatar
kraus committed
503

504 505
            rinit_m = r[0];
            prinit_m = peo[0];
kraus's avatar
kraus committed
506

507 508 509 510
        } else {
            *gmsg << "Not yet supported." << endl;
            return false;
        }
kraus's avatar
kraus committed
511

512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535
        // initialize sigma matrices (for each angle one) (first guess)
        initialize(tunes.second,ravg);

        // for writing
        std::ofstream writeMturn, writeMcyc, writeMsc;

        if (write_m) {

            std::string energy = float2string(E_m);

            writeMturn.open("data/OneTurnMapForEnergy"+energy+"MeV.dat",std::ios::app);
            writeMsc.open("data/SpaceChargeMapPerAngleForEnergy"+energy+"MeV.dat",std::ios::app);
            writeMcyc.open("data/CyclotronMapPerAngleForEnergy"+energy+"MeV.dat",std::ios::app);

            writeMturn << "--------------------------------" << std::endl;
            writeMturn << "Iteration: 0 " << std::endl;
            writeMturn << "--------------------------------" << std::endl;

            writeMsc << "--------------------------------" << std::endl;
            writeMsc << "Iteration: 0 " << std::endl;
            writeMsc << "--------------------------------" << std::endl;
        }

        // calculate only for a single sector (a nSector_-th) of the whole cyclotron
536
        for (size_type i = 0; i < nSteps_m; ++i) {
537
            if (!harmonic) {
frey_m's avatar
frey_m committed
538 539 540
                Mcycs[i] = mapgen.generateMap(H_m(h[i],
                                                  h[i]*h[i]+fidx[i],
                                                  -fidx[i]),
541 542
                                              ds[i],truncOrder_m);

frey_m's avatar
frey_m committed
543 544 545 546
                Mscs[i]  = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),
                                                    sigmas_m[i](2,2),
                                                    sigmas_m[i](4,4)),
                                              ds[i],truncOrder_m);
547
            } else {
frey_m's avatar
frey_m committed
548 549 550 551
                Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),
                                                   sigmas_m[i](2,2),
                                                   sigmas_m[i](4,4)),
                                             const_ds,truncOrder_m);
552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567
            }

            if (write_m) {
                writeMcyc << Mcycs[i] << std::endl;
                writeMsc << Mscs[i] << std::endl;
            }
        }

        // one turn matrix
        mapgen.combine(Mscs,Mcycs);
        matrix_type Mturn = mapgen.getMap();

        if (write_m)
            writeMturn << Mturn << std::endl;

        // (inverse) transformation matrix
frey_m's avatar
frey_m committed
568
        sparse_matrix_type R(6, 6), invR(6, 6);
569 570 571 572 573 574 575

        // new initial sigma matrix
        matrix_type newSigma(6,6);

        // for exiting loop
        bool stop = false;

frey_m's avatar
frey_m committed
576
        value_type weight = 0.5;
577 578 579

        while (error_m > accuracy && !stop) {
            // decouple transfer matrix and compute (inverse) tranformation matrix
580
            decouple(Mturn,R,invR);
581 582

            // construct new initial sigma-matrix
583
            newSigma = updateInitialSigma(Mturn, R, invR);
584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600

            // compute new sigma matrices for all angles (except for initial sigma)
            updateSigma(Mscs,Mcycs);

            // compute error
            error_m = L2ErrorNorm(sigmas_m[0],newSigma);

            // write new initial sigma-matrix into vector
            sigmas_m[0] = weight*newSigma + (1.0-weight)*sigmas_m[0];

            if (write_m) {
                writeMsc << "--------------------------------" << std::endl;
                writeMsc << "Iteration: " << niterations_m + 1 << std::endl;
                writeMsc << "--------------------------------" << std::endl;
            }

            // compute new space charge maps
601
            for (size_type i = 0; i < nSteps_m; ++i) {
602
                if (!harmonic) {
frey_m's avatar
frey_m committed
603 604 605 606
                    Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),
                                                       sigmas_m[i](2,2),
                                                       sigmas_m[i](4,4)),
                                                 ds[i],truncOrder_m);
607
                } else {
frey_m's avatar
frey_m committed
608 609 610 611
                    Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),
                                                       sigmas_m[i](2,2),
                                                       sigmas_m[i](4,4)),
                                                 const_ds,truncOrder_m);
612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647
                }

                if (write_m)
                    writeMsc << Mscs[i] << std::endl;
            }

            // construct new one turn transfer matrix M
            mapgen.combine(Mscs,Mcycs);
            Mturn = mapgen.getMap();

            if (write_m) {
                writeMturn << "--------------------------------" << std::endl;
                writeMturn << "Iteration: " << niterations_m + 1 << std::endl;
                writeMturn << "--------------------------------" << std::endl;
                writeMturn << Mturn << std::endl;
            }

            // check if number of iterations has maxit exceeded.
            stop = (niterations_m++ > maxit);
        }

        // store converged sigma-matrix
        sigma_m.resize(6,6,false);
        sigma_m.swap(newSigma);

        // returns if the sigma matrix has converged
        converged_m = error_m < accuracy;

        // Close files
        if (write_m) {
            writeMturn.close();
            writeMsc.close();
            writeMcyc.close();
        }

    } catch(const std::exception& e) {
cortes_c's avatar
cortes_c committed
648
        std::cerr << e.what() << std::endl;
649
    }
kraus's avatar
kraus committed
650

651 652 653
    if ( converged_m && write_m ) {
        // write tunes
        std::ofstream writeSigmaMatched("data/MatchedDistributions.dat", std::ios::app);
kraus's avatar
kraus committed
654

655
        std::array<double,3> emit = this->getEmittances();
kraus's avatar
kraus committed
656

657 658 659 660 661 662 663 664 665 666 667 668 669 670
        writeSigmaMatched << "* Converged (Ex, Ey, Ez) = (" << emit[0]
                          << ", " << emit[1] << ", " << emit[2]
                          << ") pi mm mrad for E= " << E_m << " (MeV)"
                          << std::endl << "* Sigma-Matrix " << std::endl;

        for(unsigned int i = 0; i < sigma_m.size1(); ++ i) {
            writeSigmaMatched << std::setprecision(4)  << std::setw(11)
                              << sigma_m(i,0);
            for(unsigned int j = 1; j < sigma_m.size2(); ++ j) {
                writeSigmaMatched << " & " <<  std::setprecision(4)
                                  << std::setw(11) << sigma_m(i,j);
            }
            writeSigmaMatched << " \\\\" << std::endl;
        }
kraus's avatar
kraus committed
671

672 673
        writeSigmaMatched << std::endl;
        writeSigmaMatched.close();
674
    }
675 676

    return converged_m;
677
}
678

679

680
template<typename Value_type, typename Size_type>
681 682
void SigmaGenerator<Value_type, Size_type>::eigsolve_m(const matrix_type& Mturn,
                                                       sparse_matrix_type& R)
frey_m's avatar
frey_m committed
683
{
684 685 686 687
    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;
688
    typedef boost::numeric::ublas::vector<complex_t> complex_vector_type;
kraus's avatar
kraus committed
689

690 691 692 693
    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);
kraus's avatar
kraus committed
694

695 696 697 698 699 700
    // go to GSL
    for (size_type i = 0; i < 6; ++i){
        for (size_type j = 0; j < 6; ++j) {
            gsl_matrix_set(M, i, j, Mturn(i,j));
        }
    }
kraus's avatar
kraus committed
701

frey_m's avatar
frey_m committed
702
    /*int status = */gsl_eigen_nonsymmv(M, evals, evecs, wspace);
kraus's avatar
kraus committed
703

frey_m's avatar
frey_m committed
704 705 706
//     if ( !status )
//         throw OpalException("SigmaGenerator::eigsolve_m()",
//                             "Couldn't perform eigendecomposition!");
kraus's avatar
kraus committed
707

708
    /*status = *///gsl_eigen_nonsymmv_sort(evals, evecs, GSL_EIGEN_SORT_ABS_ASC);
kraus's avatar
kraus committed
709

frey_m's avatar
frey_m committed
710 711 712
//     if ( !status )
//         throw OpalException("SigmaGenerator::eigsolve_m()",
//                             "Couldn't sort eigenvalues and eigenvectors!");
kraus's avatar
kraus committed
713

714 715 716
    // go to UBLAS
    for( size_type i = 0; i < 6; i++){
        gsl_vector_complex_view evec_i = gsl_matrix_complex_column(evecs, i);
kraus's avatar
kraus committed
717

718 719 720 721 722 723
        for(size_type 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;
        }
    }
kraus's avatar
kraus committed
724

725 726 727 728 729 730
    // Sorting the Eigenvectors
    // This is an arbitrary threshold that has worked for me. (We should fix this)
    value_type threshold = 10e-12;
    bool isZdirection = false;
    std::vector<complex_vector_type> zVectors{};
    std::vector<complex_vector_type> xyVectors{};
kraus's avatar
kraus committed
731

732 733
    for(size_type i = 0; i < 6; i++){
        complex_t z = R(i,0);
kraus's avatar
kraus committed
734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751
        if(std::abs(z) < threshold) z = 0.;
        if(z == 0.) isZdirection = true;
        complex_vector_type v(6);
        if(isZdirection){
            for(size_type j = 0;j < 6; j++){
                complex_t z = R(i,j);
                v(j) = z;
            }
            zVectors.push_back(v);
        }
        else{
            for(size_type j = 0; j < 6; j++){
                complex_t z = R(i,j);
                v(j) = z;
            }
            xyVectors.push_back(v);
        }
        isZdirection = false;
752
    }
kraus's avatar
kraus committed
753

754 755 756 757
    //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.");
kraus's avatar
kraus committed
758

759 760 761 762
    // Norms the Eigenvectors
     for(size_type i = 0; i < 4; i++){
        value_type norm{0};
        for(size_type j = 0; j < 6; j++) norm += std::norm(xyVectors[i](j));
kraus's avatar
kraus committed
763
        for(size_type j = 0; j < 6; j++) xyVectors[i](j) /= std::sqrt(norm);
764 765 766 767
    }
    for(size_type i = 0; i < 2; i++){
        value_type norm{0};
        for(size_type j = 0; j < 6; j++) norm += std::norm(zVectors[i](j));
kraus's avatar
kraus committed
768
        for(size_type j = 0; j < 6; j++) zVectors[i](j) /= std::sqrt(norm);
769
    }
kraus's avatar
kraus committed
770

771 772 773 774 775 776
    for(value_type 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);
kraus's avatar
kraus committed
777
        R(i,5) = xyVectors[3](i);
778
    }
kraus's avatar
kraus committed
779

780 781 782 783 784
    gsl_vector_complex_free(evals);
    gsl_matrix_complex_free(evecs);
    gsl_eigen_nonsymmv_free(wspace);
    gsl_matrix_free(M);
}
785 786


787 788 789 790 791
template<typename Value_type, typename Size_type>
bool SigmaGenerator<Value_type, Size_type>::invertMatrix_m(const sparse_matrix_type& R,
                                                           sparse_matrix_type& invR)
{
    typedef boost::numeric::ublas::permutation_matrix<size_type> pmatrix_t;
kraus's avatar
kraus committed
792

793 794
    //creates a working copy of R
    cmatrix_type A(R);
kraus's avatar
kraus committed
795

796 797
    //permutation matrix for the LU-factorization
    pmatrix_t pm(A.size1());
kraus's avatar
kraus committed
798

799 800
    //LU-factorization
    int res = lu_factorize(A,pm);
kraus's avatar
kraus committed
801

802 803
    if( res != 0)
        return false;
kraus's avatar
kraus committed
804

805 806
    // create identity matrix of invR
    invR.assign(boost::numeric::ublas::identity_matrix<complex_t>(A.size1()));
kraus's avatar
kraus committed
807

808 809
    // backsubstitute to get the inverse
    boost::numeric::ublas::lu_substitute(A, pm, invR);
kraus's avatar
kraus committed
810

811 812
    return true;
}
813 814


815 816 817 818
template<typename Value_type, typename Size_type>
void SigmaGenerator<Value_type, Size_type>::decouple(const matrix_type& Mturn,
                                                     sparse_matrix_type& R,
                                                     sparse_matrix_type& invR)
kraus's avatar
kraus committed
819
{
820
    this->eigsolve_m(Mturn, R);
kraus's avatar
kraus committed
821

822 823 824
    if ( !this->invertMatrix_m(R, invR) )
        throw OpalException("SigmaGenerator::decouple()",
                            "Couldn't compute inverse matrix!");
825
}
826

827
template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
828 829 830 831
typename SigmaGenerator<Value_type, Size_type>::value_type
SigmaGenerator<Value_type, Size_type>::isEigenEllipse(const matrix_type& Mturn,
                                                      const matrix_type& sigma)
{
832
    // compute sigma matrix after one turn
frey_m's avatar
frey_m committed
833 834 835
    matrix_type newSigma = matt_boost::gemmm<matrix_type>(Mturn,
                                                          sigma,
                                                          boost::numeric::ublas::trans(Mturn));
836 837 838

    // return L2 error
    return L2ErrorNorm(sigma,newSigma);
839 840 841
}

template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
842 843 844
inline typename SigmaGenerator<Value_type, Size_type>::matrix_type&
SigmaGenerator<Value_type, Size_type>::getSigma()
{
845
    return sigma_m;
846 847 848
}

template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
849 850 851
inline typename SigmaGenerator<Value_type, Size_type>::size_type
SigmaGenerator<Value_type, Size_type>::getIterations() const
{
852
    return (converged_m) ? niterations_m : size_type(0);
853 854
}

adelmann's avatar
adelmann committed
855
template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
856 857 858
inline typename SigmaGenerator<Value_type, Size_type>::value_type
SigmaGenerator<Value_type, Size_type>::getError() const
{
adelmann's avatar
adelmann committed
859 860 861 862
    return error_m;
}

template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
863 864 865
inline std::array<Value_type,3>
SigmaGenerator<Value_type, Size_type>::getEmittances() const
{
adelmann's avatar
adelmann committed
866
    value_type bgam = gamma_m*beta_m;
gsell's avatar
gsell committed
867
    return std::array<value_type,3>{{
frey_m's avatar
frey_m committed
868 869 870 871
        emittance_m[0]/Physics::pi/bgam,
        emittance_m[1]/Physics::pi/bgam,
        emittance_m[2]/Physics::pi/bgam
    }};
adelmann's avatar
adelmann committed
872 873
}

874 875 876 877 878
// -----------------------------------------------------------------------------------------------------------------------
// PRIVATE MEMBER FUNCTIONS
// -----------------------------------------------------------------------------------------------------------------------

template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
879 880 881
void SigmaGenerator<Value_type, Size_type>::initialize(value_type nuz,
                                                       value_type ravg)
{
882
    /*
frey_m's avatar
frey_m committed
883 884
     * The initialization is based on the analytical solution of
     * using a spherical symmetric beam in the paper:
885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901
     * Transverse-longitudinal coupling by space charge in cyclotrons
     * by Dr. Christian Baumgarten
     * (formulas: (46), (56), (57))
     */


    /* Units:
     * ----------------------------------------------
     * [wo] = 1/s
     * [nh] = 1
     * [q0] = e
     * [I] = A
     * [eps0] = (A*s)^{2}/(N*m^{2})
     * [E0] = MeV/(c^{2}) (with speed of light c)
     * [beta] = 1
     * [gamma] = 1
     * [m] = kg
kraus's avatar
kraus committed
902
     *
903 904 905 906 907 908 909 910 911 912
     * [lam] = m
     * [K3] = m
     * [alpha] = 10^{3}/(pi*mrad)
     * ----------------------------------------------
     */

    // helper constants
    value_type invbg = 1.0 / (beta_m * gamma_m);
    value_type micro = 1.0e-6;
    value_type mega = 1.0e6;
913
    //value_type kilo = 1.0e3;
914 915

    // convert mass m_m from MeV/c^2 to eV*s^{2}/m^{2}
916
    value_type m = m_m * mega/(Physics::c * Physics::c);        // [m] = eV*s^{2}/m^{2}, [m_m] = MeV/c^2
917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939

    /* Emittance [ex] = [ey] = [ez] = mm mrad (emittance_m are normalized emittances
     * (i.e. emittance multiplied with beta*gamma)
     */
    value_type ex = emittance_m[0] * invbg;                        // [ex] = mm mrad
    value_type ey = emittance_m[1] * invbg;                        // [ey] = mm mrad
    value_type ez = emittance_m[2] * invbg;                        // [ez] = mm mrad

    // convert normalized emittances: mm mrad --> m rad (mm mrad: millimeter milliradian)
    ex *= micro;
    ey *= micro;
    ez *= micro;

    // initial guess of emittance, [e] = m rad
    value_type e = std::cbrt(ex * ey * ez);             // cbrt computes cubic root (C++11) <cmath>

    // cyclotron radius [rcyc] = m
    value_type rcyc = ravg / beta_m;

    // "average" inverse bending radius
    value_type h = 1.0 / ravg;            // [h] = 1/m

    // formula (57)
940 941 942
    value_type lam = 2.0 * Physics::pi * Physics::c / (wo_m * nh_m); // wavelength, [lam] = m
    value_type K3 = 3.0 * /* physics::q0 */ 1.0 * I_m * lam / (20.0 * std::sqrt(5.0) * Physics::pi * Physics::epsilon_0 * m *
                    Physics::c * Physics::c * Physics::c * beta_m * beta_m * gamma2_m * gamma_m);               // [K3] = m
kraus's avatar
kraus committed
943

944
    value_type alpha = /* physics::q0 */ 1.0 * Physics::mu_0 * I_m / (5.0 * std::sqrt(10.0) * m * Physics::c *
945
                       gamma_m * nh_m) * std::sqrt(rcyc * rcyc * rcyc / (e * e * e));                           // [alpha] = 1/rad --> [alpha] = 1
kraus's avatar
kraus committed
946

947 948 949 950 951 952 953 954 955
    value_type sig0 = std::sqrt(2.0 * rcyc * e) / gamma_m;                                                      // [sig0] = m*sqrt(rad) --> [sig0] = m

    // formula (56)
    value_type sig;                                     // [sig] = m
    if (alpha >= 2.5) {
        sig = sig0 * std::cbrt(1.0 + alpha);            // cbrt computes cubic root (C++11) <cmath>
    } else if (alpha >= 0) {
        sig = sig0 * (1 + alpha * (0.25 - 0.03125 * alpha));
    } else {
frey_m's avatar
frey_m committed
956 957
        throw OpalException("SigmaGenerator::initialize()",
                            "Negative alpha value: " + std::to_string(alpha) + " < 0");
958
    }
959 960 961 962

    // K = Kx = Ky = Kz
    value_type K = K3 * gamma_m / (3.0 * sig * sig * sig);   // formula (46), [K] = 1/m^{2}
    value_type kx = h * h * gamma2_m;                        // formula (46) (assumption of an isochronous cyclotron), [kx] = 1/m^{2}
kraus's avatar
kraus committed
963

964 965
    value_type a = 0.5 * kx - K;    // formula (22) (with K = Kx = Kz), [a] = 1/m^{2}
    value_type b = K * K;           // formula (22) (with K = Kx = Kz and kx = h^2*gamma^2), [b] = 1/m^{4}
kraus's avatar
kraus committed
966 967


968 969
    // b must be positive, otherwise no real-valued frequency
    if (b < 0)
frey_m's avatar
frey_m committed
970 971
        throw OpalException("SigmaGenerator::initialize()",
                            "Negative value --> No real-valued frequency.");
972 973 974

    value_type tmp = a * a - b;           // [tmp] = 1/m^{4}
    if (tmp < 0)
frey_m's avatar
frey_m committed
975 976
        throw OpalException("SigmaGenerator::initialize()",
                            "Square root of negative number.");
kraus's avatar
kraus committed
977

978
    tmp = std::sqrt(tmp);               // [tmp] = 1/m^{2}
kraus's avatar
kraus committed
979

980
    if (a < tmp)
frey_m's avatar
frey_m committed
981 982
        throw OpalException("Error in SigmaGenerator::initialize()",
                            "Square root of negative number.");
kraus's avatar
kraus committed
983

adelmann's avatar
adelmann committed
984
    if (h * h * nuz * nuz <= K)
frey_m's avatar
frey_m committed
985 986
        throw OpalException("SigmaGenerator::initialize()",
                            "h^{2} * nu_{z}^{2} <= K (i.e. square root of negative number)");
987 988 989 990 991 992 993

    value_type Omega = std<