SigmaGenerator.h 12.9 KB
Newer Older
1 2
//
// Class: SigmaGenerator.h
3
// The SigmaGenerator class uses the class <b>ClosedOrbitFinder</b> to get the parameters(inverse bending radius, path length
4
// field index and tunes) to initialize the sigma matrix.
5 6
// The main function of this class is <b>match(double, unsigned int)</b>, where it iteratively tries to find a matched
// distribution for given
7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
// 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/>.
//

32 33 34 35
#ifndef SIGMAGENERATOR_H
#define SIGMAGENERATOR_H

#include <array>
adelmann's avatar
adelmann committed
36
#include <string>
37
#include <utility>
38 39
#include <vector>

40 41
#include "AbsBeamline/Cyclotron.h"
#include "FixedAlgebra/FTps.h"
42
#include "Physics/Physics.h"
43

44 45 46 47
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/vector.hpp>

frey_m's avatar
frey_m committed
48

49
/// @brief This class computes the matched distribution
50 51
class SigmaGenerator
{
frey_m's avatar
frey_m committed
52 53
public:
    /// Type for storing maps
54
    typedef boost::numeric::ublas::matrix<double> matrix_type;
kraus's avatar
kraus committed
55

56
    typedef std::complex<double> complex_t;
kraus's avatar
kraus committed
57

58 59
    /// Type for storing complex matrices
    typedef boost::numeric::ublas::matrix<complex_t> cmatrix_type;
kraus's avatar
kraus committed
60

frey_m's avatar
frey_m committed
61
    /// Type for storing the sparse maps
62
    typedef boost::numeric::ublas::compressed_matrix<complex_t,
frey_m's avatar
frey_m committed
63 64 65
                                                     boost::numeric::ublas::row_major
                                                     > sparse_matrix_type;
    /// Type for storing vectors
66
    typedef boost::numeric::ublas::vector<double> vector_type;
frey_m's avatar
frey_m committed
67
    /// Container for storing the properties for each angle
68
    typedef std::vector<double> container_type;
frey_m's avatar
frey_m committed
69
    /// Type of the truncated powere series
70
    typedef FTps<double,2*3> Series;
frey_m's avatar
frey_m committed
71
    /// Type of a map
72
    typedef FVps<double,2*3> Map;
frey_m's avatar
frey_m committed
73
    /// Type of the Hamiltonian for the cyclotron
74
    typedef std::function<Series(double,double,double)> Hamiltonian;
frey_m's avatar
frey_m committed
75
    /// Type of the Hamiltonian for the space charge
76
    typedef std::function<Series(double,double,double)> SpaceCharge;
frey_m's avatar
frey_m committed
77 78 79 80 81 82 83 84 85

    /// 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$
86
     * @param cycl is the cyclotron element
frey_m's avatar
frey_m committed
87 88
     * @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)
89
     * @param Nsectors is the number of sectors that the field map is averaged over (closed orbit computation)
frey_m's avatar
frey_m committed
90 91 92
     * @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.
     */
93 94 95
    SigmaGenerator(double I, double ex, double ey, double ez,
                   double E, double m, const Cyclotron* cycl,
                   unsigned int N, unsigned int Nsectors, unsigned int truncOrder, bool write = true);
frey_m's avatar
frey_m committed
96 97 98 99 100 101 102 103 104

    /// 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
105
     * @param denergy energy step size for closed orbit finder [MeV]
106
     * @param rguess value of radius for closed orbit finder
frey_m's avatar
frey_m committed
107 108 109
     * @param type specifies the magnetic field format (e.g. CARBONCYCL)
     * @param full match over full turn not just single sector
     */
110 111
    bool match(double accuracy, unsigned int maxit, unsigned int maxitOrbit,
               Cyclotron* cycl, double denergy, double rguess, bool full);
kraus's avatar
kraus committed
112

113 114 115 116
    /*!
     * 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
117
     */
118
    void eigsolve_m(const matrix_type& Mturn, sparse_matrix_type& R);
kraus's avatar
kraus committed
119

frey_m's avatar
frey_m committed
120
    /*!
121 122 123 124 125 126 127 128 129
     * @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
130
     *
frey_m's avatar
frey_m committed
131
     * @param Mturn is a 6x6 dimensional one turn transfer matrix
132 133
     * @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
134
     */
135
    void decouple(const matrix_type& Mturn, sparse_matrix_type& R, sparse_matrix_type& invR);
frey_m's avatar
frey_m committed
136 137 138 139 140 141 142

    /// 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
     */
143
    double isEigenEllipse(const matrix_type& Mturn, const matrix_type& sigma);
frey_m's avatar
frey_m committed
144 145 146 147 148

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

    /// Returns the number of iterations needed for convergence (if not converged, it returns zero)
149
    unsigned int getIterations() const;
frey_m's avatar
frey_m committed
150 151

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

    /// 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)
155
    std::array<double,3> getEmittances() const;
kraus's avatar
kraus committed
156

frey_m's avatar
frey_m committed
157 158 159
    const double& getInjectionRadius() const {
        return rinit_m;
    }
kraus's avatar
kraus committed
160

frey_m's avatar
frey_m committed
161 162 163
    const double& getInjectionMomentum() const {
        return prinit_m;
    }
164

165
private:
frey_m's avatar
frey_m committed
166
    /// Stores the value of the current, \f$ \left[I\right] = A \f$
167
    double I_m;
frey_m's avatar
frey_m committed
168
    /// Stores the desired emittances, \f$ \left[\varepsilon_{x}\right] = \left[\varepsilon_{y}\right] = \left[\varepsilon_{z}\right] = m \ rad \f$
169
    std::array<double,3> emittance_m;
frey_m's avatar
frey_m committed
170
    /// Is the orbital frequency, \f$ \left[\omega_{o}\right] = \frac{1}{s} \f$
171
    double wo_m;
frey_m's avatar
frey_m committed
172
    /// Stores the user-define energy, \f$ \left[E\right] = MeV \f$
173
    double E_m;
frey_m's avatar
frey_m committed
174
    /// Relativistic factor (which can be computed out ot the kinetic energy and rest mass (potential energy)), \f$ \left[\gamma\right] = 1 \f$
175
    double gamma_m;
frey_m's avatar
frey_m committed
176
    /// Relativistic factor squared
177
    double gamma2_m;
frey_m's avatar
frey_m committed
178
    /// Harmonic number, \f$ \left[N_{h}\right] = 1 \f$
179
    double nh_m;
frey_m's avatar
frey_m committed
180
    /// Velocity (c/v), \f$ \left[\beta\right] = 1 \f$
181
    double beta_m;
frey_m's avatar
frey_m committed
182
    /// Is the mass of the particles, \f$ \left[m\right] = \frac{MeV}{c^{2}} \f$
183
    double m_m;
frey_m's avatar
frey_m committed
184
    /// Is the number of iterations needed for convergence
185
    unsigned int niterations_m;
frey_m's avatar
frey_m committed
186 187 188
    /// Is true if converged, false otherwise
    bool converged_m;
    /// Minimum energy needed in cyclotron, \f$ \left[E_{min}\right] = MeV \f$
189
    double Emin_m;
frey_m's avatar
frey_m committed
190
    /// Maximum energy reached in cyclotron, \f$ \left[E_{max}\right] = MeV \f$
191
    double Emax_m;
frey_m's avatar
frey_m committed
192
    /// Number of integration steps for closed orbit computation
193
    unsigned int N_m;
194
    /// Number of (symmetric) sectors the field is averaged over
195
    unsigned int nSectors_m;
frey_m's avatar
frey_m committed
196
    /// Number of integration steps per sector (--> also: number of maps per sector)
197
    unsigned int nStepsPerSector_m;
kraus's avatar
kraus committed
198

frey_m's avatar
frey_m committed
199
    /// Number of integration steps in total
200
    unsigned int nSteps_m;
kraus's avatar
kraus committed
201

frey_m's avatar
frey_m committed
202
    /// Error of computation
203
    double error_m;
kraus's avatar
kraus committed
204

frey_m's avatar
frey_m committed
205
    /// Truncation order of the power series for the Hamiltonian (cyclotron and space charge)
206
    unsigned int truncOrder_m;
kraus's avatar
kraus committed
207

frey_m's avatar
frey_m committed
208 209
    /// Decides for writing output (default: true)
    bool write_m;
kraus's avatar
kraus committed
210

frey_m's avatar
frey_m committed
211 212
    /// Stores the stationary distribution (sigma matrix)
    matrix_type sigma_m;
213

frey_m's avatar
frey_m committed
214 215
    /// Vector storing the second moment matrix for each angle
    std::vector<matrix_type> sigmas_m;
216

frey_m's avatar
frey_m committed
217 218
    /// Stores the Hamiltonian of the cyclotron
    Hamiltonian H_m;
kraus's avatar
kraus committed
219

frey_m's avatar
frey_m committed
220 221
    /// Stores the Hamiltonian for the space charge
    SpaceCharge Hsc_m;
kraus's avatar
kraus committed
222

frey_m's avatar
frey_m committed
223 224
    /// 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
225

frey_m's avatar
frey_m committed
226 227 228 229 230 231
    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
232
     *
frey_m's avatar
frey_m committed
233 234 235
     * @param nuz is the vertical tune
     * @param ravg is the average radius of the closed orbit
     */
236
    void initialize(double, double);
kraus's avatar
kraus committed
237

frey_m's avatar
frey_m committed
238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
    /// 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
     */
261
    double L2ErrorNorm(const matrix_type&, const matrix_type&);
kraus's avatar
kraus committed
262 263


frey_m's avatar
frey_m committed
264 265 266 267 268
    /// 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
     */
269
    double L1ErrorNorm(const matrix_type&, const matrix_type&);
frey_m's avatar
frey_m committed
270

271
    double LInfErrorNorm(const matrix_type&, const matrix_type&);
272

frey_m's avatar
frey_m committed
273 274 275 276
    /// Transforms a floating point value to a string
    /*!
     * @param val is the floating point value which is transformed to a string
     */
277
    std::string float2string(double val);
kraus's avatar
kraus committed
278

frey_m's avatar
frey_m committed
279 280 281 282 283 284 285 286 287 288
    /// 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
     */
289 290 291
    void writeOrbitOutput_m(const std::pair<double,double>& tunes,
                            const double& ravg,
                            const double& freqError,
frey_m's avatar
frey_m committed
292 293 294 295 296
                            const container_type& r_turn,
                            const container_type& peo,
                            const container_type& h_turn,
                            const container_type&  fidx_turn,
                            const container_type&  ds_turn);
297 298

    void writeMatrix(std::ofstream&, const matrix_type&);
299 300
};

301

302 303 304
inline
typename SigmaGenerator::matrix_type&
SigmaGenerator::getSigma()
frey_m's avatar
frey_m committed
305
{
306
    return sigma_m;
307
}
308

309

310 311
inline
unsigned int SigmaGenerator::getIterations() const
frey_m's avatar
frey_m committed
312
{
313
    return (converged_m) ? niterations_m : 0;
314 315 316
}


317 318
inline
double SigmaGenerator::getError() const
frey_m's avatar
frey_m committed
319
{
adelmann's avatar
adelmann committed
320 321 322
    return error_m;
}

323 324 325

inline
std::array<double,3> SigmaGenerator::getEmittances() const
frey_m's avatar
frey_m committed
326
{
327 328
    double bgam = gamma_m*beta_m;
    return std::array<double,3>{{
frey_m's avatar
frey_m committed
329 330 331
        emittance_m[0] / Physics::pi / bgam * 1.0e6,
        emittance_m[1] / Physics::pi / bgam * 1.0e6,
        emittance_m[2] / Physics::pi / bgam * 1.0e6
frey_m's avatar
frey_m committed
332
    }};
adelmann's avatar
adelmann committed
333 334
}

kraus's avatar
kraus committed
335
#endif