diff --git a/src/Distribution/CMakeLists.txt b/src/Distribution/CMakeLists.txt
index 81cf1cdc64d776363db893c99d65af9b3b1c6c81..15a96fc25cd292ca7a9f1a71e0e7d27b5925850a 100644
--- a/src/Distribution/CMakeLists.txt
+++ b/src/Distribution/CMakeLists.txt
@@ -12,7 +12,6 @@ add_opal_sources(${_SRCS})
set (HDRS
ClosedOrbitFinder.h
Distribution.h
- Harmonics.h
LaserProfile.h
MapGenerator.h
matrix_vector_operation.h
diff --git a/src/Distribution/Distribution.cpp b/src/Distribution/Distribution.cpp
index e4e1648269780fa2a23fed21d75af2c7d709ac4e..2778e4631b70abb6a62463c28ad985011a061f3f 100644
--- a/src/Distribution/Distribution.cpp
+++ b/src/Distribution/Distribution.cpp
@@ -1297,7 +1297,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
CyclotronElement,
denergy,
rguess,
- false, full)) {
+ full)) {
std::array<double,3> Emit = siggen->getEmittances();
diff --git a/src/Distribution/Harmonics.h b/src/Distribution/Harmonics.h
deleted file mode 100644
index dbd7356ec7d8c85918aff78d1daece910b4bff10..0000000000000000000000000000000000000000
--- a/src/Distribution/Harmonics.h
+++ /dev/null
@@ -1,553 +0,0 @@
-//
-// Class Harmonics
-// This class computes the cyclotron map based on harmonics.
-// All functions are copied and translated to C++ from the original program inj2_ana.c of Dr. C. Baumgarten.
-//
-// Copyright (c) 2014 - 2015, Christian Baumgarten, Paul Scherrer Institut, Villigen PSI, Switzerland
-// Matthias Frey, ETH Zürich
-// All rights reserved
-//
-// 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/>.
-//
-#ifndef HARMONICS_H
-#define HARMONICS_H
-
-#include <cmath>
-#include <fstream>
-#include <string>
-#include <utility>
-#include <vector>
-
-#include "Physics/Physics.h"
-
-#include <boost/numeric/ublas/matrix.hpp>
-
-#include "matrix_vector_operation.h"
-
-template<typename Value_type, typename Size_type>
-class Harmonics
-{
-public:
- /// Type of variable
- typedef Value_type value_type;
- /// Type of size
- typedef Size_type size_type;
- /// Type for specifying matrices
- typedef boost::numeric::ublas::matrix<value_type> matrix_type;
-
-
- /// Defines starting point of computation
- enum START { SECTOR_ENTRANCE, SECTOR_CENTER, SECTOR_EXIT, VALLEY_CENTER };
-
- /// Constructor
- /*!
- * @param wo is the nominal orbital frequency in Hz
- * @param Emin is the minimal energy in MeV
- * @param Emax is the maximal energy in MeV
- * @param nth is the number of angle steps
- * @param nr is the number of radial steps
- * @param nSector is the number of sectors
- * @param E is the kinetic energy [MeV]
- * @param E0 is the potential energy [MeV]
- */
- Harmonics(value_type, value_type, value_type, size_type, size_type, size_type, value_type, value_type);
-
- /// Compute all maps of the cyclotron using harmonics
- /*!
- * @param filename1 name of file 1
- * @param filename2 name of file 2
- * @param startmod (in center of valley, start of sector, ...)
- */
- std::vector<matrix_type> computeMap(const std::string, const std::string, const int);
-
- /// Returns the radial and vertical tunes
- std::pair<value_type, value_type> getTunes();
-
- /// Returns the radius
- value_type getRadius();
-
- /// Returns step size
- value_type getPathLength();
-
-private:
-
- /// Stores the nominal orbital frequency in Hz
- value_type wo_m;
- /// Stores the minimum energy in MeV
- value_type Emin_m;
- /// Stores the maximum energy in MeV
- value_type Emax_m;
- /// Stores the number of angle splits
- size_type nth_m;
- /// Stores the number of radial splits
- size_type nr_m;
- /// Stores the number of sectors
- size_type nSector_m;
- /// Stores the tunes (radial, vertical)
- std::pair<value_type, value_type> tunes_m;
- /// Stores the radius
- value_type radius_m;
- /// Stores the path length
- value_type ds_m;
- /// Stores the energy for which we perform the computation
- value_type E_m;
- /// Stores the potential energy [MeV]
- value_type E0_m;
-
- /// Compute some matrix (ask Dr. C. Baumgarten)
- matrix_type __Mix6(value_type, value_type, value_type);
- /// Compute some matrix (ask Dr. C. Baumgarten)
- matrix_type __Md6(value_type, value_type);
- /// Compute some matrix (ask Dr. C. Baumgarten)
- matrix_type __Mb6(value_type, value_type, value_type);
- /// Compute some matrix (ask Dr. C. Baumgarten)
- matrix_type __Mb6k(value_type, value_type, value_type, value_type);
-};
-
-// -----------------------------------------------------------------------------------------------------------------------
-// PUBLIC MEMBER FUNCTIONS
-// -----------------------------------------------------------------------------------------------------------------------
-
-template<typename Value_type, typename Size_type>
-Harmonics<Value_type, Size_type>::Harmonics(value_type wo, value_type Emin, value_type Emax,
- size_type nth, size_type nr, size_type nSector, value_type E, value_type E0)
- : wo_m(wo), Emin_m(Emin), Emax_m(Emax), nth_m(nth), nr_m(nr), nSector_m(nSector), E_m(E), E0_m(E0)
-{}
-
-template<typename Value_type, typename Size_type>
-std::vector<typename Harmonics<Value_type, Size_type>::matrix_type> Harmonics<Value_type, Size_type>::computeMap(const std::string filename1, const std::string filename2, const int startmod) {
- // i --> different energies
- // j --> different angles
-
- // ---------------
-
-
-// unsigned int nth_m = 360; //270;
- unsigned int N = 4;
- value_type two_pi = 2.0 * M_PI;
- value_type tpiN = two_pi / value_type(N);
- value_type piN = M_PI / value_type(N);
-// unsigned int nr_m = 180;
- unsigned int cnt = 0;
- bool set = false;
-
- value_type gap = 0.05;
- value_type K1 = 0.45;
-
- std::vector<matrix_type> Mcyc(nth_m);
-
- value_type beta_m;
-
- std::vector<value_type> gamma(nr_m), E(nr_m), PC(nr_m), nur(nr_m), nuz(nr_m);
-
- std::vector<value_type> R(nr_m), r(nr_m);
- std::vector<value_type> Bmag(nr_m), k(nr_m);
- std::vector<value_type> alpha(nr_m), beta(nr_m);
-
- //----------------------------------------------------------------------
- // read files
- std::ifstream infile2(filename2);
- std::ifstream infile1(filename1);
-
- if (!infile1.is_open() || !infile2.is_open()) {
- std::cerr << "Couldn't open files!" << std::endl;
- std::exit(0);
- }
-
- unsigned int n = 0;
- value_type rN1, cN1, sN1, aN1, phN1, rN2, cN2, sN2, aN2, phN2;
-
- while (infile1 >> rN1 >> cN1 >> sN1 >> aN1 >> phN1 &&
- infile2 >> rN2 >> cN2 >> sN2 >> aN2 >> phN2) {
-
- R[n] = rN1 * 0.01;
- alpha[n] = 2.0 * std::acos(aN2 / aN1) / value_type(N);
- beta[n] = std::atan(sN1 / cN1) / value_type(N);
- value_type f4 = aN1 * M_PI * 0.5 / std::sin(0.5 * alpha[n] * value_type(N));
- value_type f8 = aN2 * M_PI / std::sin(alpha[n] * value_type(N));
-
- Bmag[n] = 0.5 * (f4 + f8);
-
- n++;
- }
-
- infile1.close();
- infile2.close();
- //----------------------------------------------------------------------
-
- value_type s = std::sin(piN);
- value_type c = std::cos(piN);
-
- std::vector<value_type> dalp(nr_m), dbet(nr_m), len2(nr_m), len(nr_m), t1eff(nr_m), t2eff(nr_m);
- std::vector<value_type> t1(nr_m), t2(nr_m), t3(nr_m);
-
- dalp[0] = (alpha[1] - alpha[0]) / (R[1] - R[0]);
- dalp[n-1] = (alpha[n-1] - alpha[n-2]) / (R[n-1] - R[n-2]);
- dbet[0] = (beta[1] - alpha[0]) / (R[1] - R[0]);
- dbet[n-1] = (beta[n-1] - beta[n-2]) / (R[n-1] - R[n-2]);
-
- for (unsigned int i = 1; i < n - 1; ++i) {
- dalp[i] = R[i] * (alpha[i+1] - alpha[i-1]) / (R[i+1] - R[i-1]);
- dbet[i] = R[i] * (beta[i+1] - beta[i-1]) / (R[i+1] - R[i-1]);
- }
-
- for (unsigned int i = 0; i < n; ++i) {
- value_type cot = 1.0 / std::tan(0.5 * alpha[i]);
-
- value_type eps1 = std::atan(dbet[i] - 0.5 * dalp[i]);
- value_type eps2 = std::atan(dbet[i] + 0.5 * dalp[i]);
-
- value_type phi = piN - 0.5 * alpha[i];
-
- // bending radius
- r[i] = R[i] * sin(0.5 * alpha[i]) / s;
- len2[i] = R[i] * std::sin(phi);
- len[i] = two_pi * r[i] + 2.0 * len2[i] * value_type(N);
-
- value_type g1 = phi + eps1;
- value_type g2 = - phi + eps2;
-
- t1[i] = std::tan(g1);
- t2[i] = std::tan(g2);
- t3[i] = std::tan(phi);
-
- value_type fac = gap / r[i] * K1;
- value_type psi = fac * (1.0 + std::sin(g1) * std::sin(g1)) / std::cos(g1);
-
- t1eff[i] = std::tan(g1 - psi);
-
- psi = fac * (1.0 + std::sin(g2) * std::sin(g2)) / std::cos(g2);
- t2eff[i] = std::tan(g2 + psi);
-
- beta_m = wo_m / Physics::c * len[i] / two_pi;
- gamma[i] = 1.0 / std::sqrt(1.0 - beta_m * beta_m);
- E[i] = E0_m * (gamma[i] - 1.0);
- PC[i] = E0_m * gamma[i] * beta_m;
- Bmag[i] = E0_m * 1.0e6 / Physics::c * beta_m * gamma[i] / r[i] * 10.0;
-
- if (!set && E[i] >= E_m) {
- set = true;
- cnt = i;
- }
- }
-
- // (average) field gradient
- for (unsigned int i = 0; i < n; ++i) {
- value_type dBdr;
- if (i == 0)
- dBdr = (Bmag[1] - Bmag[0]) / (R[1] - R[0]);
- else if (i == n - 1)
- dBdr = (Bmag[n-1] - Bmag[n-2]) / (R[n-1] - R[n-2]);
- else
- dBdr = (Bmag[i+1] - Bmag[i-1]) / (R[i+1] - R[i-1]);
-
- value_type bb = std::sin(piN - 0.5 * alpha[i]) / std::sin(0.5 * alpha[i]);
- value_type eps = 2.0 * bb / (1.0 + bb * bb);
- value_type BaV = Bmag[i];
-
- k[i] = r[i] * r[i] / (R[i] * BaV) * dBdr * (1.0 + bb * value_type(N) / M_PI * s);
- }
-
- for (unsigned int i = 0; i < n; ++i) {
- if ((k[1] > -0.999) && (E[i] > Emin_m) && (E[i] < Emax_m + 1.0)) {
- value_type kx = std::sqrt(1.0 + k[i]);
- value_type Cx = std::cos(tpiN * kx);
- value_type Sx = std::sin(tpiN * kx);
- value_type pm, ky, Cy, Sy;
-
- if (k[i] >= 0.0) {
- pm = 1.0;
- ky = std::sqrt(std::abs(k[i]));
- Cy = std::cosh(tpiN * ky);
- Sy = std::sinh(tpiN * ky);
- } else {
- pm = -1.0;
- ky = std::sqrt(std::abs(k[i]));
- Cy = std::cos(tpiN * ky);
- Sy = std::sin(tpiN * ky);
- }
-
- value_type tav = 0.5 * (t1[i] - t2[i]);
- value_type tmul = t1[i] * t2[i];
- value_type lrho = 2.0 * len2[i] / r[i];
-
- nur[i] = std::acos(Cx * (1.0 + lrho * tav) + Sx / kx * (tav - lrho * 0.5 * (kx * kx + tmul))) / tpiN;
-
- tav = 0.5 * (t1eff[i] - t2eff[i]);
- tmul = t1eff[i] * t2eff[i];
-
- nuz[i] = std::acos(Cy * (1.0 - lrho * tav) + Sy / ky * (0.5 * lrho * (pm * ky * ky - tmul) - tav)) / tpiN;
-
-// std::cout << E[i] << " " << R[i] << " " << nur[i] << " " << nuz[i] << std::endl;
- }
- }
- // -------------------------------------------------------------------------
-
- int i = cnt;
-
- tunes_m = { nur[i], nuz[i] };
- radius_m = R[i];
-
-// for (int i = 9; i < 10; ++i) { // n = 1 --> only for one energy
- value_type smax, slim, ss, gam2;
-
- smax = len[i] / value_type(N);
- slim = r[i] * tpiN;
- ds_m = smax / value_type(nth_m);
- ss = 0.0;
- gam2 = gamma[i] * gamma[i];
-
- matrix_type Mi = __Mix6( t1[i], t1eff[i], r[i]);
- matrix_type Mx = __Mix6(- t2[i], - t2eff[i], r[i]);
- matrix_type Md = __Md6(ds_m, gam2);
- matrix_type Mb = __Mb6k(ds_m / r[i], r[i], k[i], gam2);
-
- unsigned int j = 0;
-
- matrix_type M0 = boost::numeric::ublas::zero_matrix<value_type>(6);
- matrix_type M1 = boost::numeric::ublas::zero_matrix<value_type>(6);
-
- switch (startmod) {
- case 0: // SECTOR_ENTRANCE
- ss = slim - ds_m;
- Mcyc[j++] = boost::numeric::ublas::prod(Mb, Mi);
-
- while (ss > ds_m) {
- Mcyc[j++] = Mb;
- ss -= ds_m;
- }
-
- M0 = __Mb6k(ss / r[i], r[i], k[i], gam2);
- M1 = __Md6(ds_m - ss, gam2);
-
- Mcyc[j++] = matt_boost::gemmm<matrix_type>(M1, Mx, M0);
-
- ss = smax - slim - (ds_m - ss);
-
- while ((ss > ds_m) && (j < nth_m)) {
- Mcyc[j++] = Md;
- ss -= ds_m;
- }
-
- if (j < nth_m)
- Mcyc[j++] = __Md6(ss, gam2);
-
- break;
-
- case 1: // SECTOR_CENTER
- j = 0;
- ss = slim * 0.5;
-
- while (ss > ds_m) {
- Mcyc[j++] = Mb;
- ss -= ds_m;
- }
-
- M0 = __Mb6k(ss / r[i], r[i], k[i], gam2);
- M1 = __Md6(ds_m - ss, gam2);
-
- Mcyc[j++] = matt_boost::gemmm<matrix_type>(M1, Mx, M0);
-
- ss = smax - slim - (ds_m - ss);
-
- while (ss > ds_m) {
- Mcyc[j++] = Md;
- ss -= ds_m;
- }
-
- M0 = __Md6(ss, gam2);
- M1 = __Mb6k((ds_m - ss) / r[i], r[i], k[i], gam2);
-
- Mcyc[j++] = matt_boost::gemmm<matrix_type>(M1, Mi, M0);
-
- ss = slim * 0.5 - (ds_m - ss);
-
- while ((ss > ds_m) && (j < nth_m)) {
- Mcyc[j++] = Mb;
- ss -= ds_m;
- }
-
- if (j < nth_m)
- Mcyc[j++] = __Mb6k(ss / r[i], r[i], k[i], gam2);
-
- break;
-
- case 2: // SECTOR_EXIT
- j = 0;
-
- Mcyc[j++] = boost::numeric::ublas::prod(Md,Mx);
-
- ss = smax - slim - ds_m;
-
- while (ss > ds_m) {
- Mcyc[j++] = Md;
- ss -= ds_m;
- }
-
- M0 = __Md6(ss, gam2);
- M1 = __Mb6k((ds_m - ss) / r[i], r[i], k[i], gam2);
-
- Mcyc[j++] = matt_boost::gemmm<matrix_type>(M1, Mi, M0);
-
- ss = slim - (ds_m - ss);
-
- while ((ss > ds_m) && (j < nth_m)) {
- Mcyc[j++] = Mb;
- ss -= ds_m;
- }
-
- if (j < nth_m)
- Mcyc[j++] = __Mb6k(ss / r[i], r[i], k[i], gam2);
-
- break;
-
- case 3: // VALLEY_CENTER
- default:
- j = 0;
- ss = (smax - slim) * 0.5;
-
- while (ss > ds_m) {
- Mcyc[j++] = Md;
- ss -= ds_m;
- }
-
- M0 = __Md6(ss, gam2);
- M1 = __Mb6k((ds_m - ss) / r[i], r[i], k[i], gam2);
-
- Mcyc[j++] = matt_boost::gemmm<matrix_type>(M1, Mi, M0);
-
- ss = slim - (ds_m - ss);
- while (ss > ds_m) {
- Mcyc[j++] = Mb;
- ss -= ds_m;
- }
-
- M0 = __Mb6k(ss / r[i], r[i], k[i], gam2);
- M1 = __Md6(ds_m - ss, gam2);
-
- Mcyc[j++] = matt_boost::gemmm<matrix_type>(M1, Mx, M0);
-
- ss = (smax - slim) * 0.5 - (ds_m - ss);
- while ((ss > ds_m) && (j < nth_m)) {
- Mcyc[j++] = Md;
- ss -= ds_m;
- }
-
- if (j < nth_m)
- Mcyc[j++] = __Md6(ss, gam2);
-
- break;
- } // end of switch
-// } // end of for
-
- return std::vector<matrix_type>(Mcyc);
-}
-
-template<typename Value_type, typename Size_type>
-inline std::pair<Value_type, Value_type> Harmonics<Value_type, Size_type>::getTunes() {
- return tunes_m;
-}
-
-template<typename Value_type, typename Size_type>
-inline typename Harmonics<Value_type, Size_type>::value_type Harmonics<Value_type, Size_type>::getRadius() {
- return radius_m;
-}
-
-template<typename Value_type, typename Size_type>
-inline typename Harmonics<Value_type, Size_type>::value_type Harmonics<Value_type, Size_type>::getPathLength() {
- return ds_m;
-}
-
-
-// -----------------------------------------------------------------------------------------------------------------------
-// PRIVATE MEMBER FUNCTIONS
-// -----------------------------------------------------------------------------------------------------------------------
-
-template<typename Value_type, typename Size_type>
-typename Harmonics<Value_type, Size_type>::matrix_type Harmonics<Value_type, Size_type>::__Mix6(value_type tx, value_type ty, value_type r) {
- matrix_type M = boost::numeric::ublas::identity_matrix<value_type>(6);
-
- M(1,0) = tx / r;
- M(3,2) = - ty / r;
-
- return M;
-}
-
-template<typename Value_type, typename Size_type>
-typename Harmonics<Value_type, Size_type>::matrix_type Harmonics<Value_type, Size_type>::__Md6(value_type l, value_type gam2) {
- matrix_type M = boost::numeric::ublas::identity_matrix<value_type>(6);
-
- M(0,1) = l;
- M(2,3) = l;
- M(4,5) = l / gam2;
-
- return M;
-}
-
-template<typename Value_type, typename Size_type>
-typename Harmonics<Value_type, Size_type>::matrix_type Harmonics<Value_type, Size_type>::__Mb6(value_type phi, value_type r, value_type gam2) {
- matrix_type M = boost::numeric::ublas::identity_matrix<value_type>(6);
-
- value_type C = std::cos(phi);
- value_type S = std::sin(phi);
- value_type l =r * phi;
-
- M(0,0) = M(1,1) = C;
- M(0,1) = S * r;
- M(1,0) = - S / r;
- M(2,3) = l;
- M(4,5) = l / gam2 - l + r * S;
- M(4,0) = - (M(1,5) = S);
- M(4,1) = - (M(0,5) = r * (1.0 - C));
-
- return M;
-}
-
-template<typename Value_type, typename Size_type>
-typename Harmonics<Value_type, Size_type>::matrix_type Harmonics<Value_type, Size_type>::__Mb6k(value_type phi, value_type r, value_type k, value_type gam2) {
- matrix_type M = boost::numeric::ublas::identity_matrix<value_type>(6);
-
- value_type C,S;
-
- if (k == 0.0) return __Mb6(phi,r,gam2);
-
- value_type fx = std::sqrt(1.0 + k);
- value_type fy = std::sqrt(std::abs(k));
- value_type c = std::cos(phi * fx);
- value_type s = std::sin(phi * fx);
- value_type l = r * phi;
-
-
- if (k > 0.0) {
- C = std::cosh(phi * fy);
- S = std::sinh(phi * fy);
- } else {
- C = std::cos(phi * fy);
- S = std::sin(phi * fy);
- }
-
- M(0,0) = M(1,1) = c;
- M(0,1) = s * r / fx;
- M(1,0) = - s * fx / r;
- M(2,2) = M(3,3) = C;
- M(2,3) = S * r / fy;
-
- value_type sign = (std::signbit(k)) ? value_type(-1) : value_type(1);
-
- M(3,2) = sign * S * fy / r;
- M(4,5) = l / gam2 - r / (1.0 + k) * (phi - s / fx);
- M(4,0)= - (M(1,5) = s / fx);
- M(4,1)= - (M(0,5) = r * (1.0 - c) / (1.0 + k));
-
- return M;
-}
-
-#endif
\ No newline at end of file
diff --git a/src/Distribution/SigmaGenerator.h b/src/Distribution/SigmaGenerator.h
index a8d6e148de9358daceff6641d6b7cef2b12bfc6c..c3b43acc1eea6e568653aac157afee1aa668175a 100644
--- a/src/Distribution/SigmaGenerator.h
+++ b/src/Distribution/SigmaGenerator.h
@@ -70,7 +70,6 @@
#include "matrix_vector_operation.h"
#include "ClosedOrbitFinder.h"
#include "MapGenerator.h"
-#include "Harmonics.h"
extern Inform *gmsg;
@@ -140,11 +139,10 @@ public:
* @param denergy energy step size for closed orbit finder [MeV]
* @param rguess value of radius for closed orbit finder
* @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
*/
bool match(value_type accuracy, size_type maxit, size_type maxitOrbit,
- Cyclotron* cycl, value_type denergy, value_type rguess, bool harmonic, bool full);
+ Cyclotron* cycl, value_type denergy, value_type rguess, bool full);
/*!
* Eigenvalue / eigenvector solver
@@ -198,7 +196,7 @@ public:
return prinit_m;
}
- private:
+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$
@@ -304,6 +302,8 @@ public:
*/
value_type L1ErrorNorm(const matrix_type&, const matrix_type&);
+ value_type LInfErrorNorm(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
@@ -328,6 +328,8 @@ public:
const container_type& h_turn,
const container_type& fidx_turn,
const container_type& ds_turn);
+
+ void writeMatrix(std::ofstream&, const matrix_type&);
};
// -----------------------------------------------------------------------------------------------------------------------
@@ -442,7 +444,7 @@ template<typename Value_type, typename Size_type>
Cyclotron* cycl,
value_type denergy,
value_type rguess,
- bool harmonic, bool full)
+ bool full)
{
/* compute the equilibrium orbit for energy E_
* and get the the following properties:
@@ -467,61 +469,61 @@ template<typename Value_type, typename Size_type>
std::vector<matrix_type> Mcycs(nSteps_m), Mscs(nSteps_m);
container_type h(nSteps_m), r(nSteps_m), ds(nSteps_m), fidx(nSteps_m);
- value_type ravg = 0.0, const_ds = 0.0;
- std::pair<value_type,value_type> tunes;
- if (!harmonic) {
- ClosedOrbitFinder<value_type, size_type,
- boost::numeric::odeint::runge_kutta4<container_type> > cof(m_m, N_m, cycl, false, nSectors_m);
+ ClosedOrbitFinder<value_type, size_type,
+ boost::numeric::odeint::runge_kutta4<container_type> > cof(m_m, N_m, cycl, false, nSectors_m);
- if ( !cof.findOrbit(accuracy, maxitOrbit, E_m, denergy, rguess) ) {
- throw OpalException("SigmaGenerator::match()",
- "Closed orbit finder didn't converge.");
- }
+ if ( !cof.findOrbit(accuracy, maxitOrbit, E_m, denergy, rguess) ) {
+ throw OpalException("SigmaGenerator::match()",
+ "Closed orbit finder didn't converge.");
+ }
+
+ cof.computeOrbitProperties(E_m);
- cof.computeOrbitProperties(E_m);
-
- // properties of one turn
- tunes = cof.getTunes();
- ravg = cof.getAverageRadius(); // average radius
-
- value_type angle = cycl->getPHIinit();
- container_type h_turn = cof.getInverseBendingRadius(angle);
- container_type r_turn = cof.getOrbit(angle);
- container_type ds_turn = cof.getPathLength(angle);
- container_type fidx_turn = cof.getFieldIndex(angle);
-
- container_type peo = cof.getMomentum(angle);
-
- // 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
- << "* Closed orbit info:" << endl
- << "*" << endl
- << "* average radius: " << ravg << " [m]" << endl
- << "* initial radius: " << r_turn[0] << " [m]" << endl
- << "* initial momentum: " << peo[0] << " [Beta Gamma]" << endl
- << "* frequency error: " << cof.getFrequencyError()*100 <<" [ % ] "<< endl
- << "* horizontal tune: " << tunes.first << endl
- << "* vertical tune: " << tunes.second << endl
- << "* ----------------------------" << endl << endl;
-
- // 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());
-
- rinit_m = r[0];
- prinit_m = peo[0];
-
- } else {
- *gmsg << "Not yet supported." << endl;
- return false;
+ // properties of one turn
+ std::pair<value_type,value_type> tunes = cof.getTunes();
+ value_type ravg = cof.getAverageRadius(); // average radius
+
+ value_type angle = cycl->getPHIinit();
+ container_type h_turn = cof.getInverseBendingRadius(angle);
+ container_type r_turn = cof.getOrbit(angle);
+ container_type ds_turn = cof.getPathLength(angle);
+ container_type fidx_turn = cof.getFieldIndex(angle);
+
+ container_type peo = cof.getMomentum(angle);
+
+ // 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
+ << "* Closed orbit info:" << endl
+ << "*" << endl
+ << "* average radius: " << ravg << " [m]" << endl
+ << "* initial radius: " << r_turn[0] << " [m]" << endl
+ << "* initial momentum: " << peo[0] << " [Beta Gamma]" << endl
+ << "* frequency error: " << cof.getFrequencyError()*100 <<" [ % ] "<< endl
+ << "* horizontal tune: " << tunes.first << endl
+ << "* vertical tune: " << tunes.second << endl
+ << "* ----------------------------" << endl << endl;
+
+ // 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());
+
+ rinit_m = r[0];
+ prinit_m = peo[0];
+
+ std::string fpath = Util::combineFilePath({
+ OpalData::getInstance()->getAuxiliaryOutputDirectory(),
+ "maps"
+ });
+ if (!boost::filesystem::exists(fpath)) {
+ boost::filesystem::create_directory(fpath);
}
// initialize sigma matrices (for each angle one) (first guess)
@@ -536,65 +538,54 @@ template<typename Value_type, typename Size_type>
std::string fname = Util::combineFilePath({
OpalData::getInstance()->getAuxiliaryOutputDirectory(),
- "OneTurnMapForEnergy" + energy + "MeV.dat"
+ "maps",
+ "OneTurnMapsForEnergy" + energy + "MeV.dat"
});
writeMturn.open(fname, std::ios::app);
fname = Util::combineFilePath({
OpalData::getInstance()->getAuxiliaryOutputDirectory(),
- "SpaceChargeMapPerAngleForEnergy" + energy + "MeV.dat"
+ "maps",
+ "SpaceChargeMapPerAngleForEnergy" + energy + "MeV_iteration_0.dat"
});
writeMsc.open(fname, std::ios::app);
fname = Util::combineFilePath({
OpalData::getInstance()->getAuxiliaryOutputDirectory(),
+ "maps",
"CyclotronMapPerAngleForEnergy" + energy + "MeV.dat"
});
writeMcyc.open(fname, 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
for (size_type i = 0; i < nSteps_m; ++i) {
- if (!harmonic) {
- Mcycs[i] = mapgen.generateMap(H_m(h[i],
- h[i]*h[i]+fidx[i],
- -fidx[i]),
- ds[i],truncOrder_m);
+ 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);
- }
- if (write_m) {
- writeMcyc << Mcycs[i] << std::endl;
- writeMsc << Mscs[i] << std::endl;
- }
+ 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);
+
+ writeMatrix(writeMcyc, Mcycs[i]);
+ writeMatrix(writeMsc, Mscs[i]);
}
+ if (write_m)
+ writeMsc.close();
+
// one turn matrix
mapgen.combine(Mscs,Mcycs);
matrix_type Mturn = mapgen.getMap();
- if (write_m)
- writeMturn << Mturn << std::endl;
+ writeMatrix(writeMturn, Mturn);
// (inverse) transformation matrix
sparse_matrix_type R(6, 6), invR(6, 6);
@@ -623,40 +614,40 @@ template<typename Value_type, typename Size_type>
// 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
for (size_type i = 0; i < nSteps_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);
+ 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);
+
+ if (write_m) {
+
+ std::string energy = float2string(E_m);
+
+ std::string fname = Util::combineFilePath({
+ OpalData::getInstance()->getAuxiliaryOutputDirectory(),
+ "maps",
+ "SpaceChargeMapPerAngleForEnergy" + energy + "MeV_iteration_"
+ + std::to_string(niterations_m + 1)
+ + ".dat"
+ });
+
+ writeMsc.open(fname, std::ios::out);
}
- if (write_m)
- writeMsc << Mscs[i] << std::endl;
+ writeMatrix(writeMsc, Mscs[i]);
+
+ if (write_m) {
+ writeMsc.close();
+ }
}
// 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;
- }
+ writeMatrix(writeMturn, Mturn);
// check if number of iterations has maxit exceeded.
stop = (niterations_m++ > maxit);
@@ -672,7 +663,6 @@ template<typename Value_type, typename Size_type>
// Close files
if (write_m) {
writeMturn.close();
- writeMsc.close();
writeMcyc.close();
}
@@ -1134,25 +1124,6 @@ SigmaGenerator<Value_type, Size_type>::updateInitialSigma(const matrix_type& /*M
sigma(i,i) *= -1.0;
}
- if (write_m) {
- std::string energy = float2string(E_m);
-
- std::string fname = Util::combineFilePath({
- OpalData::getInstance()->getAuxiliaryOutputDirectory(),
- "maps",
- "SigmaPerAngleForEnergy" + energy + "MeV.dat"
- });
-
- std::ofstream writeSigma(fname, std::ios::app);
-
- writeSigma << "--------------------------------" << std::endl;
- writeSigma << "Iteration: " << niterations_m + 1 << std::endl;
- writeSigma << "--------------------------------" << std::endl;
-
- writeSigma << sigma << std::endl;
- writeSigma.close();
- }
-
return sigma;
}
@@ -1170,13 +1141,17 @@ void SigmaGenerator<Value_type, Size_type>::updateSigma(const std::vector<matrix
std::string fname = Util::combineFilePath({
OpalData::getInstance()->getAuxiliaryOutputDirectory(),
"maps",
- "SigmaPerAngleForEnergy" + energy + "MeV.dat"
+ "SigmaPerAngleForEnergy" + energy + "MeV_iteration_"
+ + std::to_string(niterations_m + 1)
+ + ".dat"
});
writeSigma.open(fname,std::ios::app);
}
// initial sigma is already computed
+ writeMatrix(writeSigma, sigmas_m[0]);
+
for (size_type i = 1; i < nSteps_m; ++i) {
// transfer matrix for one angle
M = boost::numeric::ublas::prod(Mscs[i - 1],Mcycs[i - 1]);
@@ -1184,12 +1159,10 @@ void SigmaGenerator<Value_type, Size_type>::updateSigma(const std::vector<matrix
sigmas_m[i] = matt_boost::gemmm<matrix_type>(M,sigmas_m[i - 1],
boost::numeric::ublas::trans(M));
- if (write_m)
- writeSigma << sigmas_m[i] << std::endl;
+ writeMatrix(writeSigma, sigmas_m[i]);
}
if (write_m) {
- writeSigma << std::endl;
writeSigma.close();
}
}
@@ -1217,16 +1190,32 @@ typename SigmaGenerator<Value_type, Size_type>::value_type
// compute difference
matrix_type diff = newS - oldS;
- std::for_each(diff.data().begin(), diff.data().end(),
+ std::transform(diff.data().begin(), diff.data().end(), diff.data().begin(),
[](value_type& val) {
return std::abs(val);
});
- // sum squared error up and take square root
return std::accumulate(diff.data().begin(), diff.data().end(), 0.0);
}
+template<typename Value_type, typename Size_type>
+typename SigmaGenerator<Value_type, Size_type>::value_type
+ SigmaGenerator<Value_type, Size_type>::LInfErrorNorm(const matrix_type& oldS,
+ const matrix_type& newS)
+{
+ // compute difference
+ matrix_type diff = newS - oldS;
+
+ std::transform(diff.data().begin(), diff.data().end(), diff.data().begin(),
+ [](value_type& val) {
+ return std::abs(val);
+ });
+
+ return *std::max_element(diff.data().begin(), diff.data().end());
+}
+
+
template<typename Value_type, typename Size_type>
std::string SigmaGenerator<Value_type, Size_type>::float2string(value_type val) {
std::ostringstream out;
@@ -1325,4 +1314,17 @@ void SigmaGenerator<Value_type, Size_type>::writeOrbitOutput_m(
writeProperties.close();
}
+
+template<typename Value_type, typename Size_type>
+void SigmaGenerator<Value_type, Size_type>::writeMatrix(std::ofstream& os, const matrix_type& m) {
+ if (!write_m)
+ return;
+
+ for (size_type i = 0; i < 6; ++i) {
+ for (size_type j = 0; j < 6; ++j)
+ os << m(i, j) << " ";
+ }
+ os << std::endl;
+}
+
#endif
\ No newline at end of file