SigmaGenerator.h 48.6 KB
Newer Older
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
//
// Class: SigmaGenerator.h
// 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(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
// 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/>.
//

31 32 33 34 35
#ifndef SIGMAGENERATOR_H
#define SIGMAGENERATOR_H

#include <array>
#include <cmath>
36
#include <fstream>
adelmann's avatar
adelmann committed
37
#include <functional>
38
#include <iomanip>
39 40 41 42
#include <iterator>
#include <limits>
#include <list>
#include <numeric>
43
#include <sstream>
adelmann's avatar
adelmann committed
44
#include <string>
45
#include <utility>
46 47
#include <vector>

48
#include "Physics/Physics.h"
49
#include "Utilities/Options.h"
Andreas Adelmann's avatar
Andreas Adelmann committed
50 51
#include "Utilities/Options.h"
#include "Utilities/OpalException.h"
52

53 54
#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>

55 56 57 58
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/vector.hpp>

59 60 61
#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
62
#include <boost/numeric/ublas/io.hpp>
63 64 65 66 67

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

68 69 70
#include "matrix_vector_operation.h"
#include "ClosedOrbitFinder.h"
#include "MapGenerator.h"
adelmann's avatar
adelmann committed
71
#include "Harmonics.h"
72

frey_m's avatar
frey_m committed
73 74
extern Inform *gmsg;

75
/// @brief This class computes the matched distribution
76 77 78
template<typename Value_type, typename Size_type>
class SigmaGenerator
{
frey_m's avatar
frey_m committed
79 80 81 82 83 84 85 86 87
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
88

89
    typedef std::complex<value_type> complex_t;
kraus's avatar
kraus committed
90

91 92
    /// Type for storing complex matrices
    typedef boost::numeric::ublas::matrix<complex_t> cmatrix_type;
kraus's avatar
kraus committed
93

frey_m's avatar
frey_m committed
94
    /// Type for storing the sparse maps
95
    typedef boost::numeric::ublas::compressed_matrix<complex_t,
frey_m's avatar
frey_m committed
96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
                                                     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$
119
     * @param cycl is the cyclotron element
frey_m's avatar
frey_m committed
120 121
     * @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)
122
     * @param Nsectors is the number of sectors that the field map is averaged over (closed orbit computation)
frey_m's avatar
frey_m committed
123 124 125 126
     * @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,
127
                   value_type E, value_type m, const Cyclotron* cycl,
128
                   size_type N, size_type Nsectors, size_type truncOrder, bool write = true);
frey_m's avatar
frey_m committed
129 130 131 132 133 134 135 136 137

    /// 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
138
     * @param denergy energy step size for closed orbit finder [MeV]
139
     * @param rguess value of radius for closed orbit finder
frey_m's avatar
frey_m committed
140 141 142 143
     * @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
     */
144
    bool match(value_type accuracy, size_type maxit, size_type maxitOrbit,
frey_m's avatar
frey_m committed
145
               Cyclotron* cycl, value_type denergy, value_type rguess, bool harmonic, bool full);
kraus's avatar
kraus committed
146

147 148 149 150
    /*!
     * 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
151
     */
152
    void eigsolve_m(const matrix_type& Mturn, sparse_matrix_type& R);
kraus's avatar
kraus committed
153

frey_m's avatar
frey_m committed
154
    /*!
155 156 157 158 159 160 161 162 163
     * @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
164
     *
frey_m's avatar
frey_m committed
165
     * @param Mturn is a 6x6 dimensional one turn transfer matrix
166 167
     * @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
168
     */
169
    void decouple(const matrix_type& Mturn, sparse_matrix_type& R, sparse_matrix_type& invR);
frey_m's avatar
frey_m committed
170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189

    /// 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
190

frey_m's avatar
frey_m committed
191 192 193
    const double& getInjectionRadius() const {
        return rinit_m;
    }
kraus's avatar
kraus committed
194

frey_m's avatar
frey_m committed
195 196 197
    const double& getInjectionMomentum() const {
        return prinit_m;
    }
198 199

    private:
frey_m's avatar
frey_m committed
200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227
    /// 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;
228 229
    /// Number of (symmetric) sectors the field is averaged over
    size_type nSectors_m;
frey_m's avatar
frey_m committed
230 231
    /// Number of integration steps per sector (--> also: number of maps per sector)
    size_type nStepsPerSector_m;
kraus's avatar
kraus committed
232

frey_m's avatar
frey_m committed
233 234
    /// Number of integration steps in total
    size_type nSteps_m;
kraus's avatar
kraus committed
235

frey_m's avatar
frey_m committed
236 237
    /// Error of computation
    value_type error_m;
kraus's avatar
kraus committed
238

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

frey_m's avatar
frey_m committed
242 243
    /// Decides for writing output (default: true)
    bool write_m;
kraus's avatar
kraus committed
244

frey_m's avatar
frey_m committed
245 246
    /// Stores the stationary distribution (sigma matrix)
    matrix_type sigma_m;
247

frey_m's avatar
frey_m committed
248 249
    /// Vector storing the second moment matrix for each angle
    std::vector<matrix_type> sigmas_m;
250

frey_m's avatar
frey_m committed
251 252
    /// Stores the Hamiltonian of the cyclotron
    Hamiltonian H_m;
kraus's avatar
kraus committed
253

frey_m's avatar
frey_m committed
254 255
    /// Stores the Hamiltonian for the space charge
    SpaceCharge Hsc_m;
kraus's avatar
kraus committed
256

frey_m's avatar
frey_m committed
257 258
    /// 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
259

frey_m's avatar
frey_m committed
260 261 262 263 264 265
    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
266
     *
frey_m's avatar
frey_m committed
267 268 269 270
     * @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
271

frey_m's avatar
frey_m committed
272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295
    /// 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
296 297


frey_m's avatar
frey_m committed
298 299 300 301 302 303 304 305 306 307 308 309
    /// 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
310

frey_m's avatar
frey_m committed
311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328
    /// 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);
329 330 331 332 333 334 335
};

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

template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
336 337 338 339 340 341
SigmaGenerator<Value_type, Size_type>::SigmaGenerator(value_type I,
                                                      value_type ex,
                                                      value_type ey,
                                                      value_type ez,
                                                      value_type E,
                                                      value_type m,
342
                                                      const Cyclotron* cycl,
frey_m's avatar
frey_m committed
343
                                                      size_type N,
344
                                                      size_type Nsectors,
frey_m's avatar
frey_m committed
345 346
                                                      size_type truncOrder,
                                                      bool write)
347
    : I_m(I)
348
    , wo_m(cycl->getRfFrequ(0)*1E6/cycl->getCyclHarm()*2.0*Physics::pi)
349 350 351
    , E_m(E)
    , gamma_m(E/m+1.0)
    , gamma2_m(gamma_m*gamma_m)
352
    , nh_m(cycl->getCyclHarm())
353 354 355 356
    , beta_m(std::sqrt(1.0-1.0/gamma2_m))
    , m_m(m)
    , niterations_m(0)
    , converged_m(false)
357 358
    , Emin_m(cycl->getFMLowE())
    , Emax_m(cycl->getFMHighE())
359
    , N_m(N)
360
    , nSectors_m(Nsectors)
361
    , nStepsPerSector_m(N/cycl->getSymmetry())
362 363 364 365 366 367 368
    , 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)
369
{
370 371
    // set emittances (initialization like that due to old compiler version)
    // [ex] = [ey] = [ez] = pi*mm*mrad --> [emittance] = mm mrad
372 373 374
    emittance_m[0] = ex * Physics::pi;
    emittance_m[1] = ey * Physics::pi;
    emittance_m[2] = ez * Physics::pi;
375 376

    // minimum beta*gamma
377
    value_type minGamma = Emin_m / m_m + 1.0;
378 379 380 381 382 383 384
    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
385

adelmann's avatar
adelmann committed
386 387
    // Define the Hamiltonian
    Series::setGlobalTruncOrder(truncOrder_m);
kraus's avatar
kraus committed
388

adelmann's avatar
adelmann committed
389 390 391 392 393 394 395
    // 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
396

adelmann's avatar
adelmann committed
397 398 399 400 401
    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
402

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

adelmann's avatar
adelmann committed
407
        // formula (57)
408 409 410
        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
411

adelmann's avatar
adelmann committed
412
        value_type milli = 1.0e-3;
kraus's avatar
kraus committed
413

adelmann's avatar
adelmann committed
414
        // formula (30), (31)
snuverink_j's avatar
snuverink_j committed
415
        // [sigma(0,0)] = mm^{2} --> [sx] = [sy] = [sz] = mm
kraus's avatar
kraus committed
416
        // multiply with 0.001 to get meter --> [sx] = [sy] = [sz] = m
adelmann's avatar
adelmann committed
417 418 419
        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
420

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

adelmann's avatar
adelmann committed
423 424
        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
425

adelmann's avatar
adelmann committed
426 427 428
        value_type Kx = kxy / sx;
        value_type Ky = kxy / sy;
        value_type Kz = K3 * f / (tmp * sz);
kraus's avatar
kraus committed
429

adelmann's avatar
adelmann committed
430 431 432 433
        return -0.5 * Kx * x_m * x_m
               -0.5 * Ky * y_m * y_m
               -0.5 * Kz * z_m * z_m * gamma2_m;
     };
434 435 436
}

template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
437 438 439
  bool SigmaGenerator<Value_type, Size_type>::match(value_type accuracy,
                                                    size_type maxit,
                                                    size_type maxitOrbit,
440
                                                    Cyclotron* cycl,
frey_m's avatar
frey_m committed
441
                                                    value_type denergy,
frey_m's avatar
frey_m committed
442
                                                    value_type rguess,
443
                                                    bool harmonic, bool full)
444 445
{
    /* compute the equilibrium orbit for energy E_
446 447 448 449 450
     * and get the the following properties:
     * - inverse bending radius h
     * - step sizes of path ds
     * - tune nuz
     */
cortes_c's avatar
cortes_c committed
451

452
    try {
453 454
        if ( !full )
            nSteps_m = nStepsPerSector_m;
Cris_cortes's avatar
Cris_cortes committed
455

456 457 458 459 460 461
        // object for space charge map and cyclotron map
        MapGenerator<value_type,
                     size_type,
                     Series,
                     Map,
                     Hamiltonian,
462
                     SpaceCharge> mapgen(nSteps_m);
463 464

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

467
        container_type h(nSteps_m), r(nSteps_m), ds(nSteps_m), fidx(nSteps_m);
468 469 470 471 472
        value_type ravg = 0.0, const_ds = 0.0;
        std::pair<value_type,value_type> tunes;

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

frey_m's avatar
frey_m committed
475
            if ( !cof.findOrbit(accuracy, maxitOrbit, E_m, denergy, rguess) ) {
476 477 478 479
                throw OpalException("SigmaGenerator::match()",
                                    "Closed orbit finder didn't converge.");
            }

480
            cof.computeOrbitProperties(E_m);
481 482

            // properties of one turn
483 484
            tunes = cof.getTunes();
            ravg = cof.getAverageRadius();                   // average radius
kraus's avatar
kraus committed
485

486
            value_type angle = cycl->getPHIinit();
frey_m's avatar
frey_m committed
487
            container_type h_turn = cof.getInverseBendingRadius(angle);
488
            container_type r_turn = cof.getOrbit(angle);
frey_m's avatar
frey_m committed
489 490
            container_type ds_turn = cof.getPathLength(angle);
            container_type fidx_turn = cof.getFieldIndex(angle);
491

492
            container_type peo = cof.getMomentum(angle);
kraus's avatar
kraus committed
493

494 495 496 497 498 499 500
            // 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
501
                  << "* Closed orbit info:" << endl
502 503 504 505
                  << "*" << endl
                  << "* average radius: " << ravg << " [m]" << endl
                  << "* initial radius: " << r_turn[0] << " [m]" << endl
                  << "* initial momentum: " << peo[0] << " [Beta Gamma]" << endl
506
                  << "* frequency error: " << cof.getFrequencyError()*100 <<" [ % ] "<< endl
507 508 509 510
                  << "* horizontal tune: " << tunes.first << endl
                  << "* vertical tune: " << tunes.second << endl
                  << "* ----------------------------" << endl << endl;

511 512 513 514 515
            // 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
516

517 518
            rinit_m = r[0];
            prinit_m = peo[0];
kraus's avatar
kraus committed
519

520 521 522 523
        } else {
            *gmsg << "Not yet supported." << endl;
            return false;
        }
kraus's avatar
kraus committed
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
        // 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
549
        for (size_type i = 0; i < nSteps_m; ++i) {
550
            if (!harmonic) {
frey_m's avatar
frey_m committed
551 552 553
                Mcycs[i] = mapgen.generateMap(H_m(h[i],
                                                  h[i]*h[i]+fidx[i],
                                                  -fidx[i]),
554 555
                                              ds[i],truncOrder_m);

frey_m's avatar
frey_m committed
556 557 558 559
                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);
560
            } else {
frey_m's avatar
frey_m committed
561 562 563 564
                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);
565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580
            }

            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
581
        sparse_matrix_type R(6, 6), invR(6, 6);
582 583 584 585 586 587 588

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

        // for exiting loop
        bool stop = false;

frey_m's avatar
frey_m committed
589
        value_type weight = 0.5;
590 591 592

        while (error_m > accuracy && !stop) {
            // decouple transfer matrix and compute (inverse) tranformation matrix
593
            decouple(Mturn,R,invR);
594 595

            // construct new initial sigma-matrix
596
            newSigma = updateInitialSigma(Mturn, R, invR);
597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613

            // 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
614
            for (size_type i = 0; i < nSteps_m; ++i) {
615
                if (!harmonic) {
frey_m's avatar
frey_m committed
616 617 618 619
                    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);
620
                } else {
frey_m's avatar
frey_m committed
621 622 623 624
                    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);
625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660
                }

                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
661
        std::cerr << e.what() << std::endl;
662
    }
kraus's avatar
kraus committed
663

664 665 666
    if ( converged_m && write_m ) {
        // write tunes
        std::ofstream writeSigmaMatched("data/MatchedDistributions.dat", std::ios::app);
kraus's avatar
kraus committed
667

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

670 671 672 673 674 675 676 677 678 679 680 681 682 683
        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
684

685 686
        writeSigmaMatched << std::endl;
        writeSigmaMatched.close();
687
    }
688 689

    return converged_m;
690
}
691

692

693
template<typename Value_type, typename Size_type>
694 695
void SigmaGenerator<Value_type, Size_type>::eigsolve_m(const matrix_type& Mturn,
                                                       sparse_matrix_type& R)
frey_m's avatar
frey_m committed
696
{
697 698 699 700
    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;
701
    typedef boost::numeric::ublas::vector<complex_t> complex_vector_type;
kraus's avatar
kraus committed
702

703 704 705 706
    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
707

708 709 710 711 712 713
    // 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
714

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

frey_m's avatar
frey_m committed
717 718 719
//     if ( !status )
//         throw OpalException("SigmaGenerator::eigsolve_m()",
//                             "Couldn't perform eigendecomposition!");
kraus's avatar
kraus committed
720

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

frey_m's avatar
frey_m committed
723 724 725
//     if ( !status )
//         throw OpalException("SigmaGenerator::eigsolve_m()",
//                             "Couldn't sort eigenvalues and eigenvectors!");
kraus's avatar
kraus committed
726

727 728 729
    // 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
730

731 732 733 734 735 736
        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
737

738 739 740 741 742 743
    // 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
744

745 746
    for(size_type i = 0; i < 6; i++){
        complex_t z = R(i,0);
kraus's avatar
kraus committed
747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764
        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;
765
    }
kraus's avatar
kraus committed
766

767 768 769 770
    //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
771

772 773 774 775
    // 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
776
        for(size_type j = 0; j < 6; j++) xyVectors[i](j) /= std::sqrt(norm);
777 778 779 780
    }
    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
781
        for(size_type j = 0; j < 6; j++) zVectors[i](j) /= std::sqrt(norm);
782
    }
kraus's avatar
kraus committed
783

784 785 786 787 788 789
    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
790
        R(i,5) = xyVectors[3](i);
791
    }
kraus's avatar
kraus committed
792

793 794 795 796 797
    gsl_vector_complex_free(evals);
    gsl_matrix_complex_free(evecs);
    gsl_eigen_nonsymmv_free(wspace);
    gsl_matrix_free(M);
}
798 799


800 801 802 803 804
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
805

806 807
    //creates a working copy of R
    cmatrix_type A(R);
kraus's avatar
kraus committed
808

809 810
    //permutation matrix for the LU-factorization
    pmatrix_t pm(A.size1());
kraus's avatar
kraus committed
811

812 813
    //LU-factorization
    int res = lu_factorize(A,pm);
kraus's avatar
kraus committed
814

815 816
    if( res != 0)
        return false;
kraus's avatar
kraus committed
817

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

821 822
    // backsubstitute to get the inverse
    boost::numeric::ublas::lu_substitute(A, pm, invR);
kraus's avatar
kraus committed
823

824 825
    return true;
}
826 827


828 829 830 831
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
832
{
833
    this->eigsolve_m(Mturn, R);
kraus's avatar
kraus committed
834

835 836 837
    if ( !this->invertMatrix_m(R, invR) )
        throw OpalException("SigmaGenerator::decouple()",
                            "Couldn't compute inverse matrix!");
838
}
839

840
template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
841 842 843 844
typename SigmaGenerator<Value_type, Size_type>::value_type
SigmaGenerator<Value_type, Size_type>::isEigenEllipse(const matrix_type& Mturn,
                                                      const matrix_type& sigma)
{
845
    // compute sigma matrix after one turn
frey_m's avatar
frey_m committed
846 847 848
    matrix_type newSigma = matt_boost::gemmm<matrix_type>(Mturn,
                                                          sigma,
                                                          boost::numeric::ublas::trans(Mturn));
849 850 851

    // return L2 error
    return L2ErrorNorm(sigma,newSigma);
852 853 854
}

template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
855 856 857
inline typename SigmaGenerator<Value_type, Size_type>::matrix_type&
SigmaGenerator<Value_type, Size_type>::getSigma()
{
858
    return sigma_m;
859 860 861
}

template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
862 863 864
inline typename SigmaGenerator<Value_type, Size_type>::size_type
SigmaGenerator<Value_type, Size_type>::getIterations() const
{
865
    return (converged_m) ? niterations_m : size_type(0);
866 867
}

adelmann's avatar
adelmann committed
868
template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
869 870 871
inline typename SigmaGenerator<Value_type, Size_type>::value_type
SigmaGenerator<Value_type, Size_type>::getError() const
{
adelmann's avatar
adelmann committed
872 873 874 875
    return error_m;
}

template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
876 877 878
inline std::array<Value_type,3>
SigmaGenerator<Value_type, Size_type>::getEmittances() const
{
adelmann's avatar
adelmann committed
879
    value_type bgam = gamma_m*beta_m;
gsell's avatar
gsell committed
880
    return std::array<value_type,3>{{
frey_m's avatar
frey_m committed
881 882 883 884
        emittance_m[0]/Physics::pi/bgam,
        emittance_m[1]/Physics::pi/bgam,
        emittance_m[2]/Physics::pi/bgam
    }};
adelmann's avatar
adelmann committed
885 886
}

887 888 889 890 891
// -----------------------------------------------------------------------------------------------------------------------
// PRIVATE MEMBER FUNCTIONS
// -----------------------------------------------------------------------------------------------------------------------

template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
892 893 894
void SigmaGenerator<Value_type, Size_type>::initialize(value_type nuz,
                                                       value_type ravg)
{
895
    /*
frey_m's avatar
frey_m committed
896 897
     * The initialization is based on the analytical solution of
     * using a spherical symmetric beam in the paper:
898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914
     * 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
915
     *
916 917 918 919 920 921 922 923 924 925
     * [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;
926
    //value_type kilo = 1.0e3;
927 928

    // convert mass m_m from MeV/c^2 to eV*s^{2}/m^{2}
929
    value_type m = m_m * mega/(Physics::c * Physics::c);        // [m] = eV*s^{2}/m^{2}, [m_m] = MeV/c^2
930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952

    /* 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)
953 954 955
    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
956

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

960 961 962 963 964 965 966 967 968
    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
969 970
        throw OpalException("SigmaGenerator::initialize()",
                            "Negative alpha value: " + std::to_string(alpha) + " < 0");
971
    }
972 973 974 975

    // 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
976

977 978
    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
979 980


981 982
    // b must be positive, otherwise no real-valued frequency
    if (b < 0)
frey_m's avatar
frey_m committed
983 984
        throw OpalException("SigmaGenerator::initialize()",
                            "Negative value --> No real-valued frequency.");
985 986 987

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

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

993
    if (a < tmp)
frey_m's avatar
frey_m committed
994 995
        throw OpalException("Error in SigmaGenerator::initialize()",
                            "Square root of negative number.");
kraus's avatar
kraus committed
996

adelmann's avatar
adelmann committed
997
    if (h * h * nuz * nuz <= K)
frey_m's avatar
frey_m committed
998 999
        throw OpalException("SigmaGenerator::initialize()",
                            "h^{2} * nu_{z}^{2} <= K (i.e. square root of negative number)");
1000 1001 1002 1003 1004 1005 1006

    value_type Omega = std::sqrt(a + tmp);                // formula (22), [Omega] = 1/m
    value_type omega = std::sqrt(a - tmp);                // formula (22), [omega] = 1/m

    value_type A = h / (Omega * Omega + K);           // formula (26), [A] = m
    value_type B = h / (omega * omega + K);           // formula (26), [B] = m
    value_type invAB = 1.0 / (B - A);                 // [invAB] = 1/m
kraus's avatar
kraus committed
1007

1008 1009 1010
    // construct initial sigma-matrix (formula (29, 30, 31)
    // Remark: We multiply with 10^{6} (= mega) to convert emittances back.
    // 1 m^{2} = 10^{6} mm^{2}
1011