SigmaGenerator.h 48.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14
/**
 * @file 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. It further has a private class member of type
 * <b>RDM</b> for decoupling. \n
 * 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.
 *
 * @author Matthias Frey
 * @version 1.0
 */
15 16 17 18 19
#ifndef SIGMAGENERATOR_H
#define SIGMAGENERATOR_H

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

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

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

#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
42
#if BOOST_VERSION >= 106000
gsell's avatar
gsell committed
43
#include <boost/numeric/odeint/integrate/check_adapter.hpp>
44
#endif
45 46 47 48 49 50

#include "rdm.h"

#include "matrix_vector_operation.h"
#include "ClosedOrbitFinder.h"
#include "MapGenerator.h"
adelmann's avatar
adelmann committed
51
#include "Harmonics.h"
52 53 54

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

frey_m's avatar
frey_m committed
55 56
extern Inform *gmsg;

57
/// @brief This class computes the matched distribution
58 59 60
template<typename Value_type, typename Size_type>
class SigmaGenerator
{
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
    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;
        /// Type for storing the sparse maps
        typedef boost::numeric::ublas::compressed_matrix<value_type,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;
adelmann's avatar
adelmann committed
76 77 78 79 80 81 82 83
        /// 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;
84 85 86 87 88 89 90 91 92 93 94 95 96

        /// 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 wo is the orbital frequency, \f$ \left[\omega_{o}\right] = \frac{1}{s} \f$
         * @param E is the energy, \f$ \left[E\right] = MeV \f$
         * @param nh is the harmonic number
         * @param m is the mass of the particles \f$ \left[m\right] = \frac{MeV}{c^{2}} \f$
         * @param Emin is the minimum energy [MeV] needed in cyclotron, \f$ \left[E_{min}\right] = MeV \f$
         * @param Emax is the maximum energy [MeV] reached in cyclotron, \f$ \left[E_{max}\right] = MeV \f$
adelmann's avatar
adelmann committed
97
         * @param nSector is the number of sectors (symmetry assumption)
98 99
         * @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)
100
         * @param fieldmap is the location of the file that specifies the magnetic field
adelmann's avatar
adelmann committed
101
         * @param truncOrder is the truncation order for power series of the Hamiltonian
102
         * @param scaleFactor for the magnetic field (default: 1.0)
adelmann's avatar
adelmann committed
103
         * @param write is a boolean (default: true). If true all maps of all iterations are stored, otherwise not.
104
         */
105 106 107 108 109
        SigmaGenerator(value_type I, value_type ex, value_type ey, value_type ez,
                       value_type wo, value_type E, value_type nh, value_type m,
                       value_type Emin, value_type Emax, size_type nSector,
                       size_type N, const std::string& fieldmap,
                       size_type truncOrder, value_type scaleFactor = 1.0, bool write = true);
110

111
        /// Searches for a matched distribution.
112
        /*!
113
         * Returns "true" if a matched distribution within the accuracy could be found, returns "false" otherwise.
114
         * @param accuracy is used for computing the equilibrium orbit (to a certain accuracy)
115 116
         * @param maxit is the maximum number of iterations (the program stops if within this value no stationary
         *        distribution was found)
117
         * @param maxitOrbit is the maximum number of iterations for finding closed orbit
adelmann's avatar
adelmann committed
118
         * @param angle defines the start of the sector (one can choose any angle between 0° and 360°)
Andreas Adelmann's avatar
Andreas Adelmann committed
119
	 * @param guess value of radius for closed orbit finder
frey_m's avatar
frey_m committed
120
         * @param type specifies the magnetic field format (e.g. CARBONCYCL)
adelmann's avatar
adelmann committed
121
         * @param harmonic is a boolean. If "true" the harmonics are used instead of the closed orbit finder.
122
         */
123 124
        bool match(value_type accuracy, size_type maxit, size_type maxitOrbit, value_type angle,
                   value_type guess, const std::string& type, bool harmonic);
125 126

        /// Block diagonalizes the symplex part of the one turn transfer matrix
127 128
        /** It computes the transformation matrix <b>R</b> and its inverse <b>invR</b>.
         * It returns a vector containing the four eigenvalues
129 130 131 132 133 134 135
         * (alpha,beta,gamma,delta) (see formula (38), paper: Geometrical method of decoupling)
         */
        /*!
         * @param Mturn is a 6x6 dimensional one turn transfer matrix
         * @param R is the 4x4 dimensional transformation matrix (gets computed)
         * @param invR is the 4x4 dimensional inverse transformation (gets computed)
         */
136
        vector_type decouple(const matrix_type& Mturn, sparse_matrix_type& R, sparse_matrix_type& invR);
137

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

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

        /// Returns the number of iterations needed for convergence (if not converged, it returns zero)
adelmann's avatar
adelmann committed
150
        size_type getIterations() const;
kraus's avatar
kraus committed
151

adelmann's avatar
adelmann committed
152 153
        /// Returns the error (if the program didn't converged, one can call this function to check the error)
        value_type getError() const;
kraus's avatar
kraus committed
154

adelmann's avatar
adelmann committed
155 156
        /// 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;
157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185

    private:
        /// 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 (symmetric) sectors
adelmann's avatar
adelmann committed
186 187
        size_type nSector_m;
        /// Number of integration steps for closed orbit computation
188
        size_type N_m;
adelmann's avatar
adelmann committed
189 190 191 192
        /// Number of integration steps per sector (--> also: number of maps per sector)
        size_type nStepsPerSector_m;
        /// Error of computation
        value_type error_m;
kraus's avatar
kraus committed
193

194 195
        /// Location of magnetic fieldmap
        std::string fieldmap_m;
kraus's avatar
kraus committed
196

adelmann's avatar
adelmann committed
197 198
        /// Truncation order of the power series for the Hamiltonian (cyclotron and space charge)
        size_type truncOrder_m;
kraus's avatar
kraus committed
199

200 201
        /// Decides for writing output (default: true)
        bool write_m;
202 203 204
        
        /// Scale the magnetic field values (default: 1.0)
        value_type scaleFactor_m;
205 206 207 208 209 210 211 212 213

        /// Stores the stationary distribution (sigma matrix)
        matrix_type sigma_m;

        /// Vector storing the second moment matrix for each angle
        std::vector<matrix_type> sigmas_m;

        /// <b>RDM</b>-class member used for decoupling
        RDM<value_type, size_type> rdm_m;
kraus's avatar
kraus committed
214

adelmann's avatar
adelmann committed
215 216
        /// Stores the Hamiltonian of the cyclotron
        Hamiltonian H_m;
kraus's avatar
kraus committed
217

adelmann's avatar
adelmann committed
218 219
        /// Stores the Hamiltonian for the space charge
        SpaceCharge Hsc_m;
kraus's avatar
kraus committed
220

adelmann's avatar
adelmann committed
221 222
        /// All variables x, px, y, py, z, delta
        Series x_m, px_m, y_m, py_m, z_m, delta_m;
223 224 225 226 227 228 229

        /// 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.
        /*!
         * @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
230

231 232 233 234 235 236 237 238 239 240 241 242 243 244 245
        /// Reduces the 6x6 matrix to a 4x4 matrix (for more explanations consider the source code comments)
        template<class matrix>
            void reduce(matrix&);
        /// Expands the 4x4 matrix back to dimension 6x6 (for more explanations consider the source code comments)
        template<class matrix>
            void expand(matrix&);

        /// Computes the new initial sigma matrix
        /*!
         * @param M is the 6x6 one turn transfer matrix
         * @param eigen stores the 4 eigenvalues: alpha, beta, gamma, delta (in this order)
         * @param R is the transformation matrix
         * @param invR is the inverse transformation matrix
         */
        matrix_type updateInitialSigma(const matrix_type&, const vector_type&, sparse_matrix_type&, sparse_matrix_type&);
kraus's avatar
kraus committed
246

247 248 249 250 251 252 253 254 255 256 257 258 259
        /// 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&);
frey_m's avatar
frey_m committed
260 261 262 263 264 265 266 267
        
        
        /// 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&);
kraus's avatar
kraus committed
268

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

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

template<typename Value_type, typename Size_type>
kraus's avatar
kraus committed
300
SigmaGenerator<Value_type, Size_type>::SigmaGenerator(value_type I, value_type ex, value_type ey, value_type ez, value_type wo,
301 302
        value_type E, value_type nh, value_type m, value_type Emin, value_type Emax, size_type nSector,
        size_type N, const std::string& fieldmap, size_type truncOrder, value_type scaleFactor, bool write)
303
: I_m(I), wo_m(wo), E_m(E),
304
  gamma_m(E/m+1.0), gamma2_m(gamma_m*gamma_m),
305
  nh_m(nh), beta_m(std::sqrt(1.0-1.0/gamma2_m)), m_m(m), niterations_m(0), converged_m(false),
306 307 308 309
  Emin_m(Emin), Emax_m(Emax), nSector_m(nSector), N_m(N), nStepsPerSector_m(N/nSector),
  error_m(std::numeric_limits<value_type>::max()),
  fieldmap_m(fieldmap), truncOrder_m(truncOrder), write_m(write),
  scaleFactor_m(scaleFactor), sigmas_m(nStepsPerSector_m)
310
{
311 312
    // set emittances (initialization like that due to old compiler version)
    // [ex] = [ey] = [ez] = pi*mm*mrad --> [emittance] = mm mrad
313 314 315
    emittance_m[0] = ex * Physics::pi;
    emittance_m[1] = ey * Physics::pi;
    emittance_m[2] = ez * Physics::pi;
316 317

    // minimum beta*gamma
318
    value_type minGamma = Emin_m / m_m + 1.0;
319 320 321 322 323 324 325
    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
326

adelmann's avatar
adelmann committed
327 328
    // Define the Hamiltonian
    Series::setGlobalTruncOrder(truncOrder_m);
kraus's avatar
kraus committed
329

adelmann's avatar
adelmann committed
330 331 332 333 334 335 336
    // 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
337

adelmann's avatar
adelmann committed
338 339 340 341 342
    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
343

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

adelmann's avatar
adelmann committed
348
        // formula (57)
349 350 351
        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
352

adelmann's avatar
adelmann committed
353
        value_type milli = 1.0e-3;
kraus's avatar
kraus committed
354

adelmann's avatar
adelmann committed
355 356
        // formula (30), (31)
        // [sigma(0,0)] = mm^{2} rad --> [sx] = [sy] = [sz] = mm
kraus's avatar
kraus committed
357
        // multiply with 0.001 to get meter --> [sx] = [sy] = [sz] = m
adelmann's avatar
adelmann committed
358 359 360
        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
361

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

adelmann's avatar
adelmann committed
364 365
        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
366

adelmann's avatar
adelmann committed
367 368 369
        value_type Kx = kxy / sx;
        value_type Ky = kxy / sy;
        value_type Kz = K3 * f / (tmp * sz);
kraus's avatar
kraus committed
370

adelmann's avatar
adelmann committed
371 372 373 374
        return -0.5 * Kx * x_m * x_m
               -0.5 * Ky * y_m * y_m
               -0.5 * Kz * z_m * z_m * gamma2_m;
     };
375 376 377
}

template<typename Value_type, typename Size_type>
frey_m's avatar
frey_m committed
378 379 380 381 382 383 384 385
  bool SigmaGenerator<Value_type, Size_type>::match(value_type accuracy,
                                                    size_type maxit,
                                                    size_type maxitOrbit,
                                                    value_type angle,
                                                    value_type rguess,
                                                    const std::string& type,
                                                    bool harmonic)
{
386 387 388 389 390 391
    /* compute the equilibrium orbit for energy E_
     * and get the the following properties:
     * - inverse bending radius h
     * - step sizes of path ds
     * - tune nuz
     */
kraus's avatar
kraus committed
392

adelmann's avatar
adelmann committed
393
    try {
kraus's avatar
kraus committed
394

adelmann's avatar
adelmann committed
395
        // object for space charge map and cyclotron map
frey_m's avatar
frey_m committed
396 397 398 399 400 401
        MapGenerator<value_type,
                     size_type,
                     Series,
                     Map,
                     Hamiltonian,
                     SpaceCharge> mapgen(nStepsPerSector_m);
kraus's avatar
kraus committed
402

adelmann's avatar
adelmann committed
403 404
        // compute cyclotron map and space charge map for each angle and store them into a vector
        std::vector<matrix_type> Mcycs(nStepsPerSector_m), Mscs(nStepsPerSector_m);
kraus's avatar
kraus committed
405

adelmann's avatar
adelmann committed
406 407 408
        container_type h(nStepsPerSector_m), r(nStepsPerSector_m), ds(nStepsPerSector_m), fidx(nStepsPerSector_m);
        value_type ravg = 0.0, const_ds = 0.0;
        std::pair<value_type,value_type> tunes;
kraus's avatar
kraus committed
409

adelmann's avatar
adelmann committed
410
        if (!harmonic) {
411 412 413 414
            ClosedOrbitFinder<value_type, size_type,
                boost::numeric::odeint::runge_kutta4<container_type> > cof(E_m, m_m, wo_m, N_m, accuracy,
                                                                           maxitOrbit, Emin_m, Emax_m,
                                                                           nSector_m, fieldmap_m, rguess,
frey_m's avatar
frey_m committed
415
                                                                           type, scaleFactor_m, false);
kraus's avatar
kraus committed
416

adelmann's avatar
adelmann committed
417 418
            // properties of one turn
            container_type h_turn = cof.getInverseBendingRadius();
419
            container_type r_turn = cof.getOrbit(angle);
adelmann's avatar
adelmann committed
420 421
            container_type ds_turn = cof.getPathLength();
            container_type fidx_turn = cof.getFieldIndex();
frey_m's avatar
frey_m committed
422
            tunes = cof.getTunes();
adelmann's avatar
adelmann committed
423
            ravg = cof.getAverageRadius();                   // average radius
kraus's avatar
kraus committed
424

425
	    container_type peo = cof.getMomentum(angle);
frey_m's avatar
frey_m committed
426 427 428 429 430
            
            // write properties to file
            if (write_m)
                writeOrbitOutput_m(tunes, ravg, cof.getFrequencyError(),
                                   r_turn, peo, h_turn, fidx_turn, ds_turn);
431

frey_m's avatar
frey_m committed
432 433 434 435 436 437
            // write to terminal
            *gmsg << "* ----------------------------" << endl
                  << "* Closed orbit info (Gordon units):" << endl
                  << "*" << endl
                  << "* average radius: " << ravg << " [m]" << endl
                  << "* initial radius: " << r_turn[0] << " [m]" << endl
438
                  << "* initial momentum: " << peo[0] << " [Beta Gamma]" << endl
frey_m's avatar
frey_m committed
439 440 441 442
                  << "* frequency error: " << cof.getFrequencyError() << endl
                  << "* horizontal tune: " << tunes.first << endl
                  << "* vertical tune: " << tunes.second << endl
                  << "* ----------------------------" << endl << endl;
443

frey_m's avatar
frey_m committed
444 445 446 447 448 449 450 451 452 453 454
            // compute the number of steps per degree
            value_type deg_step = N_m / 360.0;
            // compute starting point of computation
            size_type start = deg_step * angle;

            // copy properties of the length of one sector (--> nStepsPerSector_m)
            std::copy_n(r_turn.begin()+start,nStepsPerSector_m, r.begin());
            std::copy_n(h_turn.begin()+start,nStepsPerSector_m, h.begin());
            std::copy_n(fidx_turn.begin()+start,nStepsPerSector_m, fidx.begin());
            std::copy_n(ds_turn.begin()+start,nStepsPerSector_m, ds.begin());
            
adelmann's avatar
adelmann committed
455
        } else {
frey_m's avatar
frey_m committed
456 457
            *gmsg << "Not yet supported." << endl;
            return false;
458
        }
kraus's avatar
kraus committed
459

Andreas Adelmann's avatar
Andreas Adelmann committed
460

frey_m's avatar
frey_m committed
461 462
        if(Options::cloTuneOnly)
            throw OpalException("Do only CLO and tune calculation","... ");
Andreas Adelmann's avatar
Andreas Adelmann committed
463

adelmann's avatar
adelmann committed
464 465
        // initialize sigma matrices (for each angle one) (first guess)
        initialize(tunes.second,ravg);
kraus's avatar
kraus committed
466

adelmann's avatar
adelmann committed
467 468
        // for writing
        std::ofstream writeMturn, writeMcyc, writeMsc;
kraus's avatar
kraus committed
469

470
        if (write_m) {
kraus's avatar
kraus committed
471

adelmann's avatar
adelmann committed
472
            std::string energy = float2string(E_m);
kraus's avatar
kraus committed
473

frey_m's avatar
frey_m committed
474 475 476
            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);
kraus's avatar
kraus committed
477

adelmann's avatar
adelmann committed
478 479 480
            writeMturn << "--------------------------------" << std::endl;
            writeMturn << "Iteration: 0 " << std::endl;
            writeMturn << "--------------------------------" << std::endl;
kraus's avatar
kraus committed
481

482
            writeMsc << "--------------------------------" << std::endl;
adelmann's avatar
adelmann committed
483
            writeMsc << "Iteration: 0 " << std::endl;
484 485
            writeMsc << "--------------------------------" << std::endl;
        }
kraus's avatar
kraus committed
486

adelmann's avatar
adelmann committed
487 488 489 490 491 492 493 494
        // calculate only for a single sector (a nSector_-th) of the whole cyclotron
        for (size_type i = 0; i < nStepsPerSector_m; ++i) {
            if (!harmonic) {
                Mcycs[i] = mapgen.generateMap(H_m(h[i],h[i]*h[i]+fidx[i],-fidx[i]),ds[i],truncOrder_m);
                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);
            } else {
                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);
            }
kraus's avatar
kraus committed
495

adelmann's avatar
adelmann committed
496 497
            if (write_m) {
                writeMcyc << Mcycs[i] << std::endl;
498
                writeMsc << Mscs[i] << std::endl;
adelmann's avatar
adelmann committed
499
            }
500
        }
kraus's avatar
kraus committed
501

adelmann's avatar
adelmann committed
502
        // one turn matrix
503
        mapgen.combine(Mscs,Mcycs);
adelmann's avatar
adelmann committed
504
        matrix_type Mturn = mapgen.getMap();
kraus's avatar
kraus committed
505

adelmann's avatar
adelmann committed
506
        if (write_m)
507
            writeMturn << Mturn << std::endl;
kraus's avatar
kraus committed
508

adelmann's avatar
adelmann committed
509 510
        // (inverse) transformation matrix
        sparse_matrix_type R, invR;
kraus's avatar
kraus committed
511

adelmann's avatar
adelmann committed
512 513
        // eigenvalues
        vector_type eigen(4);
kraus's avatar
kraus committed
514

adelmann's avatar
adelmann committed
515 516
        // new initial sigma matrix
        matrix_type newSigma(6,6);
kraus's avatar
kraus committed
517

adelmann's avatar
adelmann committed
518 519
        // for exiting loop
        bool stop = false;
kraus's avatar
kraus committed
520

adelmann's avatar
adelmann committed
521
        value_type weight = 0.05;
kraus's avatar
kraus committed
522

adelmann's avatar
adelmann committed
523 524 525
        while (error_m > accuracy && !stop) {
            // decouple transfer matrix and compute (inverse) tranformation matrix
            eigen = decouple(Mturn,R,invR);
kraus's avatar
kraus committed
526

adelmann's avatar
adelmann committed
527 528
            // construct new initial sigma-matrix
            newSigma = updateInitialSigma(Mturn,eigen,R,invR);
kraus's avatar
kraus committed
529

adelmann's avatar
adelmann committed
530 531
            // compute new sigma matrices for all angles (except for initial sigma)
            updateSigma(Mscs,Mcycs);
kraus's avatar
kraus committed
532

adelmann's avatar
adelmann committed
533 534
            // compute error
            error_m = L2ErrorNorm(sigmas_m[0],newSigma);
kraus's avatar
kraus committed
535

adelmann's avatar
adelmann committed
536 537
            // write new initial sigma-matrix into vector
            sigmas_m[0] = weight*newSigma + (1.0-weight)*sigmas_m[0];
kraus's avatar
kraus committed
538

adelmann's avatar
adelmann committed
539 540 541 542 543
            if (write_m) {
                writeMsc << "--------------------------------" << std::endl;
                writeMsc << "Iteration: " << niterations_m + 1 << std::endl;
                writeMsc << "--------------------------------" << std::endl;
            }
kraus's avatar
kraus committed
544

adelmann's avatar
adelmann committed
545 546 547 548 549 550 551
            // compute new space charge maps
            for (size_type i = 0; i < nStepsPerSector_m; ++i) {
                if (!harmonic) {
                    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);
                } else {
                    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);
                }
kraus's avatar
kraus committed
552

adelmann's avatar
adelmann committed
553 554 555
                if (write_m)
                    writeMsc << Mscs[i] << std::endl;
            }
kraus's avatar
kraus committed
556

adelmann's avatar
adelmann committed
557 558 559
            // construct new one turn transfer matrix M
            mapgen.combine(Mscs,Mcycs);
            Mturn = mapgen.getMap();
kraus's avatar
kraus committed
560

adelmann's avatar
adelmann committed
561 562 563 564 565 566
            if (write_m) {
                writeMturn << "--------------------------------" << std::endl;
                writeMturn << "Iteration: " << niterations_m + 1 << std::endl;
                writeMturn << "--------------------------------" << std::endl;
                writeMturn << Mturn << std::endl;
            }
kraus's avatar
kraus committed
567

adelmann's avatar
adelmann committed
568 569
            // check if number of iterations has maxit exceeded.
            stop = (niterations_m++ > maxit);
570
        }
kraus's avatar
kraus committed
571

adelmann's avatar
adelmann committed
572 573 574
        // store converged sigma-matrix
        sigma_m.resize(6,6,false);
        sigma_m.swap(newSigma);
kraus's avatar
kraus committed
575

adelmann's avatar
adelmann committed
576 577
        // returns if the sigma matrix has converged
        converged_m = error_m < accuracy;
kraus's avatar
kraus committed
578

adelmann's avatar
adelmann committed
579 580 581 582 583 584
        // Close files
        if (write_m) {
            writeMturn.close();
            writeMsc.close();
            writeMcyc.close();
        }
kraus's avatar
kraus committed
585

adelmann's avatar
adelmann committed
586 587
    } catch(const std::exception& e) {
        std::cerr << e.what() << std::endl;
588
    }
frey_m's avatar
frey_m committed
589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610
    
    if ( converged_m && write_m ) {
        // write tunes
        std::ofstream writeSigmaMatched("data/MatchedDistributions.dat", std::ios::app);
        
        std::array<double,3> emit = this->getEmittances();
        
        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;
        }
        
frey_m's avatar
frey_m committed
611
        writeSigmaMatched << std::endl;
frey_m's avatar
frey_m committed
612 613
        writeSigmaMatched.close();
    }
kraus's avatar
kraus committed
614

615
    return converged_m;
616 617 618 619
}

template<typename Value_type, typename Size_type>
typename SigmaGenerator<Value_type, Size_type>::vector_type SigmaGenerator<Value_type, Size_type>::decouple(const matrix_type& Mturn, sparse_matrix_type& R,
kraus's avatar
kraus committed
620
        sparse_matrix_type& invR) {
621 622 623 624 625 626 627 628 629 630 631 632 633 634 635
    // copy one turn matrix
    matrix_type M(Mturn);

    // reduce 6x6 matrix to 4x4 matrix
    reduce<matrix_type>(M);

    // compute symplex part
    matrix_type Ms = rdm_m.symplex(M);

    // diagonalize and compute transformation matrices
    rdm_m.diagonalize(Ms,R,invR);

    /*
     * formula (38) in paper of Dr. Christian Baumgarten:
     * Geometrical method of decoupling
kraus's avatar
kraus committed
636
     *
637 638 639 640
     * 		[0, 	alpha, 	0, 	0;
     * F_{d} =	-beta, 	0, 	0, 	0;
     * 		0, 	0, 	0, 	gamma;
     * 		0, 	0, 	-delta,	0]
kraus's avatar
kraus committed
641 642
     *
     *
643 644 645 646 647 648 649
     */
    vector_type eigen(4);
    eigen(0) =   Ms(0,1);       // alpha
    eigen(1) = - Ms(1,0);       // beta
    eigen(2) =   Ms(2,3);       // gamma
    eigen(3) = - Ms(3,2);       // delta
    return eigen;
650 651 652 653
}

template<typename Value_type, typename Size_type>
typename SigmaGenerator<Value_type, Size_type>::value_type SigmaGenerator<Value_type, Size_type>::isEigenEllipse(const matrix_type& Mturn, const matrix_type& sigma) {
654 655 656 657 658
    // compute sigma matrix after one turn
    matrix_type newSigma = matt_boost::gemmm<matrix_type>(Mturn,sigma,boost::numeric::ublas::trans(Mturn));

    // return L2 error
    return L2ErrorNorm(sigma,newSigma);
659 660 661 662
}

template<typename Value_type, typename Size_type>
inline typename SigmaGenerator<Value_type, Size_type>::matrix_type& SigmaGenerator<Value_type, Size_type>::getSigma() {
663
    return sigma_m;
664 665 666
}

template<typename Value_type, typename Size_type>
adelmann's avatar
adelmann committed
667
inline typename SigmaGenerator<Value_type, Size_type>::size_type SigmaGenerator<Value_type, Size_type>::getIterations() const {
668
    return (converged_m) ? niterations_m : size_type(0);
669 670
}

adelmann's avatar
adelmann committed
671 672 673 674 675 676 677 678
template<typename Value_type, typename Size_type>
inline typename SigmaGenerator<Value_type, Size_type>::value_type SigmaGenerator<Value_type, Size_type>::getError() const {
    return error_m;
}

template<typename Value_type, typename Size_type>
inline std::array<Value_type,3> SigmaGenerator<Value_type, Size_type>::getEmittances() const {
    value_type bgam = gamma_m*beta_m;
679
    return std::array<value_type,3>({emittance_m[0]/Physics::pi/bgam, emittance_m[1]/Physics::pi/bgam, emittance_m[2]/Physics::pi/bgam});
adelmann's avatar
adelmann committed
680 681
}

682 683 684 685 686
// -----------------------------------------------------------------------------------------------------------------------
// PRIVATE MEMBER FUNCTIONS
// -----------------------------------------------------------------------------------------------------------------------

template<typename Value_type, typename Size_type>
687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706
void SigmaGenerator<Value_type, Size_type>::initialize(value_type nuz, value_type ravg) {
    /*
     * The initialization is based on the analytical solution of using a spherical symmetric beam in the paper:
     * 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
707
     *
708 709 710 711 712 713 714 715 716 717 718 719 720
     * [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;
    value_type kilo = 1.0e3;

    // convert mass m_m from MeV/c^2 to eV*s^{2}/m^{2}
721
    value_type m = m_m * mega/(Physics::c * Physics::c);        // [m] = eV*s^{2}/m^{2}, [m_m] = MeV/c^2
722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744

    /* 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)
745 746 747
    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
748

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

752 753 754 755 756 757 758 759 760
    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 {
761
        throw OpalException("SigmaGenerator::initialize()", "Negative alpha value: " + std::to_string(alpha) + " < 0");
762
    }
763 764 765 766

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

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


772 773
    // b must be positive, otherwise no real-valued frequency
    if (b < 0)
774
        throw OpalException("SigmaGenerator::initialize()", "Negative value --> No real-valued frequency.");
775 776 777

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

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

782
    if (a < tmp)
783
        throw OpalException("Error in SigmaGenerator::initialize()", "Square root of negative number.");
kraus's avatar
kraus committed
784

adelmann's avatar
adelmann committed
785
    if (h * h * nuz * nuz <= K)
786
        throw OpalException("SigmaGenerator::initialize()", "h^{2} * nu_{z}^{2} <= K (i.e. square root of negative number)");
787 788 789 790 791 792 793

    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
794

795 796 797
    // 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}
798
    matrix_type sigma = boost::numeric::ublas::zero_matrix<value_type>(6);
799 800 801 802
    sigma(0,0) = invAB * (B * ex / Omega + A * ez / omega) * mega;                      // formula (30), [sigma(0,0)] = m^2 rad * 10^{6} = mm^{2} rad = mm mrad
    sigma(0,5) = sigma(5,0) = invAB * (ex / Omega + ez / omega) * mega;                 // [sigma(0,5)] = [sigma(5,0)] = m rad * 10^{6} = mm mrad	// 1000: m --> mm and 1000 to promille
    sigma(1,1) = invAB * (B * ex * Omega + A * ez * omega) * mega;                      // [sigma(1,1)] = rad * 10^{6} = mrad (and promille)
    sigma(1,4) = sigma(4,1) = invAB * (ex * Omega+ez * omega) / (K * gamma2_m) * mega;  // [sigma(1,4)] = [sigma(4,1)] = m rad * 10^{6} = mm mrad
803 804
    sigma(2,2) = ey / (std::sqrt(h * h * nuz * nuz - K)) * mega;                        // formula (31), [sigma(2,2)] = m rad * 10^{6} = mm mrad
    sigma(3,3) = (ey * mega) * (ey * mega) / sigma(2,2);
805 806
    sigma(4,4) = invAB * (A * ex * Omega + B * ez * omega) / (K * gamma2_m) * mega;     // [sigma(4,4)] = m^{2} rad * 10^{6} = mm^{2} rad = mm mrad
    sigma(5,5) = invAB * (ex / (B * Omega) + ez / (A * omega)) * mega;                  // formula (30), [sigma(5,5)] = rad * 10^{6} = mrad (and promille)
kraus's avatar
kraus committed
807

808
    // fill in initial guess of the sigma matrix (for each angle the same guess)
adelmann's avatar
adelmann committed
809
    sigmas_m.resize(nStepsPerSector_m);
810 811 812
    for (typename std::vector<matrix_type>::iterator it = sigmas_m.begin(); it != sigmas_m.end(); ++it) {
        *it = sigma;
    }
kraus's avatar
kraus committed
813

814 815 816 817 818
    if (write_m) {
        std::string energy = float2string(E_m);
        std::ofstream writeInit("data/maps/InitialSigmaPerAngleForEnergy"+energy+"MeV.dat",std::ios::app);
        writeInit << sigma << std::endl;
        writeInit.close();
819 820 821 822 823 824
    }
}

template<typename Value_type, typename Size_type>
template<class matrix>
void SigmaGenerator<Value_type, Size_type>::reduce(matrix& M) {
825
    /* The 6x6 matrix gets reduced to a 4x4 matrix in the following way:
kraus's avatar
kraus committed
826
     *
827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853
     * a11 a12 a13 a14 a15 a16
     * a21 a22 a23 a24 a25 a26          a11 a12 a15 a16
     * a31 a32 a33 a34 a35 a36  -->     a21 a22 a25 a26
     * a41 a42 a43 a44 a45 a46          a51 a52 a55 a56
     * a51 a52 a53 a54 a55 a56          a61 a62 a65 a66
     * a61 a62 a63 a64 a65 a66
     */

    // copy x- and z-direction to a 4x4 matrix_type
    matrix_type M4x4(4,4);
    for (size_type i = 0; i < 2; ++i) {
        // upper left 2x2 [a11,a12;a21,a22]
        M4x4(i,0) = M(i,0);
        M4x4(i,1) = M(i,1);
        // lower left 2x2 [a51,a52;a61,a62]
        M4x4(i + 2,0) = M(i + 4,0);
        M4x4(i + 2,1) = M(i + 4,1);
        // upper right 2x2 [a15,a16;a25,a26]
        M4x4(i,2) = M(i,4);
        M4x4(i,3) = M(i,5);
        // lower right 2x2 [a55,a56;a65,a66]
        M4x4(i + 2,2) = M(i + 4,4);
        M4x4(i + 2,3) = M(i + 4,5);
    }

    M.resize(4,4,false);
    M.swap(M4x4);
854 855 856 857
}

template<typename Value_type, typename Size_type>
template<class matrix>
858 859
void SigmaGenerator<Value_type, Size_type>::expand(matrix& M) {
    /* The 4x4 matrix gets expanded to a 6x6 matrix in the following way:
kraus's avatar
kraus committed
860
     *
861 862
     *                          a11 a12 0 0 a13 a14
     * a11 a12 a13 a14          a21 a22 0 0 a23 a24
kraus's avatar
kraus committed
863
     * a21 a22 a23 a24  -->     0   0   1 0 0   0
864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887
     * a31 a32 a33 a34          0   0   0 1 0   0
     * a41 a42 a43 a44          a31 a32 0 0 a33 a34
     *                          a41 a42 0 0 a43 a44
     */

    matrix M6x6 = boost::numeric::ublas::identity_matrix<value_type>(6,6);

    for (size_type i = 0; i < 2; ++i) {
        // upper left 2x2 [a11,a12;a21,a22]
        M6x6(i,0) = M(i,0);
        M6x6(i,1) = M(i,1);
        // lower left 2x2 [a31,a32;a41,a42]
        M6x6(i + 4,0) = M(i + 2,0);
        M6x6(i + 4,1) = M(i + 2,1);
        // upper right 2x2 [a13,a14;a23,a24]
        M6x6(i,4) = M(i,2);
        M6x6(i,5) = M(i,3);
        // lower right 2x2 [a22,a34;a43,a44]
        M6x6(i + 4,4) = M(i + 2,2);
        M6x6(i + 4,5) = M(i + 2,3);
    }

    // exchange
    M.swap(M6x6);
888 889 890 891
}

template<typename Value_type, typename Size_type>
typename SigmaGenerator<Value_type, Size_type>::matrix_type SigmaGenerator<Value_type, Size_type>::updateInitialSigma(const matrix_type& M, const vector_type& eigen,
892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925
        sparse_matrix_type& R, sparse_matrix_type& invR) {

    /*
     * This function is based on the paper of Dr. Christian Baumgarten:
     * Transverse-Longitudinal Coupling by Space Charge in Cyclotrons (2012)
     */

    /*
     * Function input:
     * - M: one turn transfer matrix
     * - eigen = {alpha, beta, gamma, delta}
     * - R: transformation matrix (in paper: E)
     * - invR: inverse transformation matrix (in paper: E^{-1}
     */

    /* formula (18):
     * sigma = -E*D*E^{-1}*S
     * with diagonal matrix D (stores eigenvalues of sigma*S (emittances apart from +- i),
     * skew-symmetric matrix (formula (13)), and tranformation matrices E, E^{-1}
     */

    // normalize emittances
    value_type invbg = 1.0 / (beta_m * gamma_m);
    value_type ex = emittance_m[0] * invbg;
    value_type ey = emittance_m[1] * invbg;
    value_type ez = emittance_m[2] * invbg;


    // alpha^2-beta*gamma = 1

    /* 0        eigen(0) 0        0
     * eigen(1) 0        0        0
     * 0        0        0        eigen(2)
     * 0        0        eigen(3) 0
kraus's avatar
kraus committed
926
     *
927
     * M = cos(mux)*[1, 0; 0, 1] + sin(mux)*[alpha, beta; -gamma, -alpha], Book, p. 242
kraus's avatar
kraus committed
928
     *
929 930 931 932 933
     * -----------------------------------------------------------------------------------
     * X-DIRECTION and Z-DIRECTION
     * -----------------------------------------------------------------------------------
     * --> eigen(0) = sin(mux)*betax
     * --> eigen(1) = -gammax*sin(mux)
kraus's avatar
kraus committed
934
     *
935
     * What is sin(mux)?   --> alphax = 0 --> -alphax^2+betax*gammax = betax*gammax = 1
kraus's avatar
kraus committed
936
     *
937
     * Convention: betax > 0
kraus's avatar
kraus committed
938
     *
939 940 941
     * 1) betax = 1/gammax
     * 2) eig0 = sin(mux)*betax
     * 3) eig1 = -gammax*sin(mux)
kraus's avatar
kraus committed
942
     *
943 944
     * eig0 = sin(mux)/gammax
     * eig1 = -gammax*sin(mux) <--> 1/gammax = -sin(mux)/eig1
kraus's avatar
kraus committed
945
     *
946 947 948 949
     * eig0 = -sin(mux)^2/eig1 --> -sin(mux)^2 = eig0*eig1      --> sin(mux) = sqrt(-eig0*eig1)
     *                                                          --> gammax = -eig1/sin(mux)
     *                                                          --> betax = eig0/sin(mux)
     */
kraus's avatar
kraus committed
950 951


952 953 954 955
    // x-direction
    value_type alphax = 0.0;
    value_type betax  = std::sqrt(std::fabs(eigen(0) / eigen(1)));
    value_type gammax = 1.0 / betax;
kraus's avatar
kraus committed
956

957 958 959 960 961 962 963 964 965
    // z-direction
    value_type alphaz = 0.0;
    value_type betaz  = std::sqrt(std::fabs(eigen(2) / eigen(3)));
    value_type gammaz = 1.0 / betaz;

    /*
     * -----------------------------------------------------------------------------------
     * Y-DIRECTION
     * -----------------------------------------------------------------------------------
kraus's avatar
kraus committed
966
     *
967 968
     * m22 m23
     * m32 m33
kraus's avatar
kraus committed
969
     *
970 971
     * m22 = cos(muy) + alpha*sin(muy)
     * m33 = cos(muy) - alpha*sin(muy)
kraus's avatar
kraus committed
972
     *
973 974
     * --> cos(muy) = 0.5*(m22 + m33)
     *     sin(muy) = sign(m32)*sqrt(1-cos(muy)^2)
kraus's avatar
kraus committed
975
     *
976
     * m22-m33 = 2*alpha*sin(muy) --> alpha = 0.5*(m22-m33)/sin(muy)
kraus's avatar
kraus committed
977
     *
978 979 980 981 982
     * m23 = betay*sin(muy)     --> betay = m23/sin(muy)
     * m32 = -gammay*sin(muy)   --> gammay = -m32/sin(muy)
     */

    value_type cosy = 0.5 * (M(2,2) + M(3,3));
983

frey_m's avatar
frey_m committed
984
    value_type sign = (std::signbit(M(2,3))) ? value_type(-1) : value_type(1);
985

frey_m's avatar
frey_m committed
986
    value_type invsiny = sign / std::sqrt(std::fabs( 1.0 - cosy * cosy));
kraus's avatar
kraus committed
987

988 989