SigmaGenerator.h 12.9 KB
 snuverink_j committed Mar 30, 2020 1 2 // // Class: SigmaGenerator.h  frey_m committed Jul 08, 2020 3 // The SigmaGenerator class uses the class ClosedOrbitFinder to get the parameters(inverse bending radius, path length  snuverink_j committed Mar 30, 2020 4 // field index and tunes) to initialize the sigma matrix.  frey_m committed Jul 08, 2020 5 6 // The main function of this class is match(double, unsigned int), where it iteratively tries to find a matched // distribution for given  snuverink_j committed Mar 30, 2020 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 . //  adelmann committed Nov 27, 2014 32 33 34 35 #ifndef SIGMAGENERATOR_H #define SIGMAGENERATOR_H #include  adelmann committed Dec 03, 2014 36 #include  adelmann committed Jan 06, 2015 37 #include  adelmann committed Nov 27, 2014 38 39 #include  frey_m committed Jul 08, 2020 40 41 #include "AbsBeamline/Cyclotron.h" #include "FixedAlgebra/FTps.h"  frey_m committed Apr 04, 2016 42 #include "Physics/Physics.h"  snuverink_j committed Mar 30, 2020 43   adelmann committed Nov 27, 2014 44 45 46 47 #include #include #include  frey_m committed Jan 19, 2016 48   adelmann committed Jan 06, 2015 49 /// @brief This class computes the matched distribution  adelmann committed Nov 27, 2014 50 51 class SigmaGenerator {  frey_m committed Aug 23, 2018 52 53 public: /// Type for storing maps  frey_m committed Jul 08, 2020 54  typedef boost::numeric::ublas::matrix matrix_type;  kraus committed May 02, 2019 55   frey_m committed Jul 08, 2020 56  typedef std::complex complex_t;  kraus committed May 02, 2019 57   frey_m committed Aug 23, 2018 58 59  /// Type for storing complex matrices typedef boost::numeric::ublas::matrix cmatrix_type;  kraus committed May 02, 2019 60   frey_m committed Aug 23, 2018 61  /// Type for storing the sparse maps  frey_m committed Aug 23, 2018 62  typedef boost::numeric::ublas::compressed_matrix sparse_matrix_type; /// Type for storing vectors  frey_m committed Jul 08, 2020 66  typedef boost::numeric::ublas::vector vector_type;  frey_m committed Aug 23, 2018 67  /// Container for storing the properties for each angle  frey_m committed Jul 08, 2020 68  typedef std::vector container_type;  frey_m committed Aug 23, 2018 69  /// Type of the truncated powere series  frey_m committed Jul 08, 2020 70  typedef FTps Series;  frey_m committed Aug 23, 2018 71  /// Type of a map  frey_m committed Jul 08, 2020 72  typedef FVps Map;  frey_m committed Aug 23, 2018 73  /// Type of the Hamiltonian for the cyclotron  frey_m committed Jul 08, 2020 74  typedef std::function Hamiltonian;  frey_m committed Aug 23, 2018 75  /// Type of the Hamiltonian for the space charge  frey_m committed Jul 08, 2020 76  typedef std::function SpaceCharge;  frey_m committed Aug 23, 2018 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$  frey_m committed Dec 21, 2018 86  * @param cycl is the cyclotron element  frey_m committed Aug 23, 2018 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)  snuverink_j committed Mar 25, 2019 89  * @param Nsectors is the number of sectors that the field map is averaged over (closed orbit computation)  frey_m committed Aug 23, 2018 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. */  frey_m committed Jul 08, 2020 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 committed Aug 23, 2018 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 committed Jan 21, 2019 105  * @param denergy energy step size for closed orbit finder [MeV]  frey_m committed Dec 21, 2018 106  * @param rguess value of radius for closed orbit finder  frey_m committed Aug 23, 2018 107 108 109  * @param type specifies the magnetic field format (e.g. CARBONCYCL) * @param full match over full turn not just single sector */  frey_m committed Jul 08, 2020 110 111  bool match(double accuracy, unsigned int maxit, unsigned int maxitOrbit, Cyclotron* cycl, double denergy, double rguess, bool full);  kraus committed May 02, 2019 112   frey_m committed Aug 23, 2018 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 committed Aug 23, 2018 117  */  frey_m committed Aug 23, 2018 118  void eigsolve_m(const matrix_type& Mturn, sparse_matrix_type& R);  kraus committed May 02, 2019 119   frey_m committed Aug 23, 2018 120  /*!  frey_m committed Aug 23, 2018 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 R and its inverse invR.  kraus committed May 02, 2019 130  *  frey_m committed Aug 23, 2018 131  * @param Mturn is a 6x6 dimensional one turn transfer matrix  frey_m committed Aug 23, 2018 132 133  * @param R is the 6x6 dimensional transformation matrix (gets computed) * @param invR is the 6x6 dimensional inverse transformation (gets computed)  frey_m committed Aug 23, 2018 134  */  frey_m committed Aug 23, 2018 135  void decouple(const matrix_type& Mturn, sparse_matrix_type& R, sparse_matrix_type& invR);  frey_m committed Aug 23, 2018 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 */  frey_m committed Jul 08, 2020 143  double isEigenEllipse(const matrix_type& Mturn, const matrix_type& sigma);  frey_m committed Aug 23, 2018 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)  frey_m committed Jul 08, 2020 149  unsigned int getIterations() const;  frey_m committed Aug 23, 2018 150 151  /// Returns the error (if the program didn't converged, one can call this function to check the error)  frey_m committed Jul 08, 2020 152  double getError() const;  frey_m committed Aug 23, 2018 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)  frey_m committed Jul 08, 2020 155  std::array getEmittances() const;  kraus committed May 02, 2019 156   frey_m committed Aug 23, 2018 157 158 159  const double& getInjectionRadius() const { return rinit_m; }  kraus committed May 02, 2019 160   frey_m committed Aug 23, 2018 161 162 163  const double& getInjectionMomentum() const { return prinit_m; }  adelmann committed Jan 06, 2015 164   frey_m committed Jul 08, 2020 165 private:  frey_m committed Aug 23, 2018 166  /// Stores the value of the current, \f$\left[I\right] = A \f$  frey_m committed Jul 08, 2020 167  double I_m;  frey_m committed Jul 08, 2020 168  /// Stores the desired emittances, \f$\left[\varepsilon_{x}\right] = \left[\varepsilon_{y}\right] = \left[\varepsilon_{z}\right] = m \ rad \f$  frey_m committed Jul 08, 2020 169  std::array emittance_m;  frey_m committed Aug 23, 2018 170  /// Is the orbital frequency, \f$\left[\omega_{o}\right] = \frac{1}{s} \f$  frey_m committed Jul 08, 2020 171  double wo_m;  frey_m committed Aug 23, 2018 172  /// Stores the user-define energy, \f$\left[E\right] = MeV \f$  frey_m committed Jul 08, 2020 173  double E_m;  frey_m committed Aug 23, 2018 174  /// Relativistic factor (which can be computed out ot the kinetic energy and rest mass (potential energy)), \f$\left[\gamma\right] = 1 \f$  frey_m committed Jul 08, 2020 175  double gamma_m;  frey_m committed Aug 23, 2018 176  /// Relativistic factor squared  frey_m committed Jul 08, 2020 177  double gamma2_m;  frey_m committed Aug 23, 2018 178  /// Harmonic number, \f$\left[N_{h}\right] = 1 \f$  frey_m committed Jul 08, 2020 179  double nh_m;  frey_m committed Aug 23, 2018 180  /// Velocity (c/v), \f$\left[\beta\right] = 1 \f$  frey_m committed Jul 08, 2020 181  double beta_m;  frey_m committed Aug 23, 2018 182  /// Is the mass of the particles, \f$\left[m\right] = \frac{MeV}{c^{2}} \f$  frey_m committed Jul 08, 2020 183  double m_m;  frey_m committed Aug 23, 2018 184  /// Is the number of iterations needed for convergence  frey_m committed Jul 08, 2020 185  unsigned int niterations_m;  frey_m committed Aug 23, 2018 186 187 188  /// Is true if converged, false otherwise bool converged_m; /// Minimum energy needed in cyclotron, \f$\left[E_{min}\right] = MeV \f$  frey_m committed Jul 08, 2020 189  double Emin_m;  frey_m committed Aug 23, 2018 190  /// Maximum energy reached in cyclotron, \f$\left[E_{max}\right] = MeV \f$  frey_m committed Jul 08, 2020 191  double Emax_m;  frey_m committed Aug 23, 2018 192  /// Number of integration steps for closed orbit computation  frey_m committed Jul 08, 2020 193  unsigned int N_m;  snuverink_j committed Mar 25, 2019 194  /// Number of (symmetric) sectors the field is averaged over  frey_m committed Jul 08, 2020 195  unsigned int nSectors_m;  frey_m committed Aug 23, 2018 196  /// Number of integration steps per sector (--> also: number of maps per sector)  frey_m committed Jul 08, 2020 197  unsigned int nStepsPerSector_m;  kraus committed May 02, 2019 198   frey_m committed Aug 23, 2018 199  /// Number of integration steps in total  frey_m committed Jul 08, 2020 200  unsigned int nSteps_m;  kraus committed May 02, 2019 201   frey_m committed Aug 23, 2018 202  /// Error of computation  frey_m committed Jul 08, 2020 203  double error_m;  kraus committed Apr 01, 2015 204   frey_m committed Aug 23, 2018 205  /// Truncation order of the power series for the Hamiltonian (cyclotron and space charge)  frey_m committed Jul 08, 2020 206  unsigned int truncOrder_m;  kraus committed Apr 01, 2015 207   frey_m committed Aug 23, 2018 208 209  /// Decides for writing output (default: true) bool write_m;  kraus committed May 02, 2019 210   frey_m committed Aug 23, 2018 211 212  /// Stores the stationary distribution (sigma matrix) matrix_type sigma_m;  frey_m committed Jul 18, 2018 213   frey_m committed Aug 23, 2018 214 215  /// Vector storing the second moment matrix for each angle std::vector sigmas_m;  adelmann committed Jan 06, 2015 216   frey_m committed Aug 23, 2018 217 218  /// Stores the Hamiltonian of the cyclotron Hamiltonian H_m;  kraus committed Apr 01, 2015 219   frey_m committed Aug 23, 2018 220 221  /// Stores the Hamiltonian for the space charge SpaceCharge Hsc_m;  kraus committed Apr 01, 2015 222   frey_m committed Aug 23, 2018 223 224  /// All variables x, px, y, py, z, delta Series x_m, px_m, y_m, py_m, z_m, delta_m;  kraus committed May 02, 2019 225   frey_m committed Aug 23, 2018 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 committed May 02, 2019 232  *  frey_m committed Aug 23, 2018 233 234 235  * @param nuz is the vertical tune * @param ravg is the average radius of the closed orbit */  frey_m committed Jul 08, 2020 236  void initialize(double, double);  kraus committed May 02, 2019 237   frey_m committed Aug 23, 2018 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&, const std::vector&); /// 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 */  frey_m committed Jul 08, 2020 261  double L2ErrorNorm(const matrix_type&, const matrix_type&);  kraus committed May 02, 2019 262 263   frey_m committed Aug 23, 2018 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 */  frey_m committed Jul 08, 2020 269  double L1ErrorNorm(const matrix_type&, const matrix_type&);  frey_m committed Aug 23, 2018 270   frey_m committed Jul 08, 2020 271  double LInfErrorNorm(const matrix_type&, const matrix_type&);  frey_m committed Jul 08, 2020 272   frey_m committed Aug 23, 2018 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 */  frey_m committed Jul 08, 2020 277  std::string float2string(double val);  kraus committed May 02, 2019 278   frey_m committed Aug 23, 2018 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 */  frey_m committed Jul 08, 2020 289 290 291  void writeOrbitOutput_m(const std::pair& tunes, const double& ravg, const double& freqError,  frey_m committed Aug 23, 2018 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);  frey_m committed Jul 08, 2020 297 298  void writeMatrix(std::ofstream&, const matrix_type&);  adelmann committed Nov 27, 2014 299 300 };  frey_m committed Jul 18, 2018 301   frey_m committed Jul 08, 2020 302 303 304 inline typename SigmaGenerator::matrix_type& SigmaGenerator::getSigma()  frey_m committed Aug 23, 2018 305 {  frey_m committed Jul 08, 2020 306  return sigma_m;  cortes_c committed Jun 03, 2018 307 }  frey_m committed Jul 18, 2018 308   adelmann committed Nov 27, 2014 309   frey_m committed Jul 08, 2020 310 311 inline unsigned int SigmaGenerator::getIterations() const  frey_m committed Aug 23, 2018 312 {  frey_m committed Jul 08, 2020 313  return (converged_m) ? niterations_m : 0;  adelmann committed Nov 27, 2014 314 315 316 }  frey_m committed Jul 08, 2020 317 318 inline double SigmaGenerator::getError() const  frey_m committed Aug 23, 2018 319 {  adelmann committed Mar 29, 2015 320 321 322  return error_m; }  frey_m committed Jul 08, 2020 323 324 325  inline std::array SigmaGenerator::getEmittances() const  frey_m committed Aug 23, 2018 326 {  frey_m committed Jul 08, 2020 327 328  double bgam = gamma_m*beta_m; return std::array{{  frey_m committed Jul 08, 2020 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 committed Aug 23, 2018 332  }};  adelmann committed Mar 29, 2015 333 334 }  kraus committed May 02, 2019 335 #endif