Commit 078763a6 authored by cortes_c's avatar cortes_c
Browse files

Fix Cyclotron transfer matrices

parent 52af4fa7
......@@ -2734,7 +2734,7 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
globalToLocal(itsBunch_m->R, phi, psi, meanR);
//globalToLocal(itsBunch_m->R, phi, psi, meanR * Vector_t(0.001));
globalToLocal(itsBunch_m->P, phi, psi); // P should only be rotated
globalToLocal(extB_m, phi, psi);
globalToLocal(extE_m, phi, psi);
}
......@@ -3546,4 +3546,4 @@ void ParallelCyclotronTracker::injectBunch_m(bool& flagTransition) {
// After this, numBunch_m is wrong but not needed anymore...
numBunch_m--;
}
}
\ No newline at end of file
}
......@@ -90,7 +90,7 @@ public:
virtual ~ParallelCyclotronTracker();
/// Apply the algorithm to an Ring
/// Apply the algorithm to a Ring
virtual void visitRing(const Ring &ring);
/// Apply the algorithm to a Cyclotorn
......@@ -267,7 +267,7 @@ private:
double sinRefTheta_m;
double cosRefTheta_m;
/// The number of bunches specified in TURNS of RUN commond
/// The number of bunches specified in TURNS of RUN command
int numBunch_m;
// 0 for single bunch (default),
......
/**
* @file ClosedOrbitFinder.h
* The algorithm is based on the paper of M. M. Gordon: "Computation of closed orbits and basic focusing properties for
* sector-focused cyclotrons and the design of 'cyclops'" (1983)
* As template arguments one chooses the type of the variables and the integrator for the ODEs. The supported steppers can
* be found on
* http://www.boost.org/ where it is part of the library Odeint.
*
* @author Matthias Frey
* @version 1.0
*/
#ifndef CLOSEDORBITFINDER_H
#define CLOSEDORBITFINDER_H
#include <algorithm>
#include <array>
#include <cmath>
#include <functional>
#include <limits>
#include <numeric>
#include <string>
#include <utility>
#include <vector>
#include "Utilities/Options.h"
#include "Utilities/Options.h"
#include "Utilities/OpalException.h"
// #include "physics.h"
#include "MagneticField.h"
// include headers for integration
#include <boost/numeric/odeint/integrate/integrate_n_steps.hpp>
/// Finds a closed orbit of a cyclotron for a given energy
template<typename Value_type, typename Size_type, class Stepper>
class ClosedOrbitFinder
{
public:
/// Type of variables
typedef Value_type value_type;
/// Type for specifying sizes
typedef Size_type size_type;
/// Type of container for storing quantities (path length, orbit, etc.)
typedef std::vector<value_type> container_type;
/// Type for holding state of ODE values
typedef std::vector<value_type> state_type;
/// Sets the initial values for the integration and calls findOrbit().
/*!
* @param E is the energy [MeV] to which the closed orbit should be found
* @param E0 is the potential energy (particle energy at rest) [MeV].
* @param wo is the nominal orbital frequency (see paper of Dr. C. Baumgarten: "Transverse-Longitudinal
* Coupling by Space Charge in Cyclotrons" (2012), formula (1))
* @param N specifies the number of splits (2pi/N), i.e number of integration steps
* @param accuracy specifies the accuracy of the closed orbit
* @param maxit is the maximal number of iterations done. Program stops if closed orbit not found within this time.
* @param Emin is the minimum energy [MeV] needed in cyclotron
* @param Emax is the maximum energy [MeV] reached in cyclotron
* @param nSector is the number of sectors (--> symmetry) of cyclotron
* @param fmapfn is the location of the file that specifies the magnetic field
* @param guess value of radius for closed orbit finder
* @param type specifies the field format (e.g. CARBONCYCL)
* @param scaleFactor for the magnetic field (default: 1.0)
* @param domain is a boolean (default: true). If "true" the closed orbit is computed over a single sector,
* otherwise over 2*pi.
*/
ClosedOrbitFinder(value_type E, value_type E0, value_type wo, size_type N,
value_type accuracy, size_type maxit, value_type Emin, value_type Emax,
size_type nSector, const std::string& fmapfn, value_type guess,
const std::string& type, value_type scaleFactor = 1.0,
bool domain = true);
/// Returns the inverse bending radius (size of container N+1)
container_type& getInverseBendingRadius();
/// Returns the step lengths of the path (size of container N+1)
container_type& getPathLength();
/// Returns the field index (size of container N+1)
container_type& getFieldIndex();
/// Returns the radial and vertical tunes (in that order)
std::pair<value_type,value_type> getTunes();
/// Returns the closed orbit (size of container N+1) starting at specific angle (only makes sense when computing
/// the closed orbit for a whole turn) (default value: 0°).
/// Attention: It computes the starting index of the array. If it's not an integer it just cuts the floating point
/// part, i.e. it takes the next starting index below. There's no interpolation of the radius.
/*!
* @param angle is the start angle for the output. Has to be within [0°,360°[ (default: 0°).
*/
container_type getOrbit(value_type angle = 0);
/// Returns the momentum of the orbit (size of container N+1)starting at specific angle (only makes sense when
/// computing the closed orbit for a whole turn) (default value: 0°), \f$ \left[ p_{r} \right] = \si{m}\f$.
/// Attention: It computes the starting index of the array. If it's not an integer it just cuts the floating point
/// part, i.e. it takes the next starting index below. There's no interpolation of the momentum.
/*!
* @param angle is the start angle for the output. Has to be within [0°,360°[ (default: 0°).
* @returns the momentum in \f$ \beta * \gamma \f$ units
*/
container_type getMomentum(value_type angle = 0);
/// Returns the relativistic factor gamma
value_type getGamma();
/// Returns the average orbit radius
value_type getAverageRadius();
/// Returns the frequency error
value_type getFrequencyError();
/// Returns true if a closed orbit could be found
bool isConverged();
private:
/// Computes the closed orbit
/*!
* @param accuracy specifies the accuracy of the closed orbit
* @param maxit is the maximal number of iterations done for finding the closed orbit
*/
bool findOrbit(value_type, size_type);
/// Fills in the values of h_m, ds_m, fidx_m. It gets called by in by constructor.
void computeOrbitProperties();
/// This function is called by the function getTunes().
/*! Transfer matrix Y = [y11, y12; y21, y22] (see Gordon paper for more details).
* @param y are the positions (elements y11 and y12 of Y)
* @param py2 is the momentum of the second solution (element y22 of Y)
* @param ncross is the number of sign changes (\#crossings of zero-line)
*/
value_type computeTune(const std::array<value_type,2>&, value_type, size_type);
/// This function computes nzcross_ which is used to compute the tune in z-direction and the frequency error
void computeVerticalOscillations();
/// Stores current position in horizontal direction for the solutions of the ODE with different initial values
std::array<value_type,2> x_m; // x_m = [x1, x2]
/// Stores current momenta in horizontal direction for the solutions of the ODE with different initial values
std::array<value_type,2> px_m; // px_m = [px1, px2]
/// Stores current position in vertical direction for the solutions of the ODE with different initial values
std::array<value_type,2> z_m; // z_m = [z1, z2]
/// Stores current momenta in vertical direction for the solutions of the ODE with different initial values
std::array<value_type,2> pz_m; // pz_m = [pz1, pz2]
/// Stores the inverse bending radius
container_type h_m;
/// Stores the step length
container_type ds_m;
/// Stores the radial orbit (size: N_m+1)
container_type r_m;
/// Stores the radial momentum
container_type pr_m;
/// Stores the field index
container_type fidx_m;
/// Counts the number of zero-line crossings in horizontal direction (used for computing horizontal tune)
size_type nxcross_m;
/// Counts the number of zero-line crossings in vertical direction (used for computing vertical tune)
size_type nzcross_m; //#crossings of zero-line in x- and z-direction
/// Is the energy for which the closed orbit should be found
value_type E_m;
/// Is the potential energy [MeV]
value_type E0_m;
/// Is the nominal orbital frequency
value_type wo_m;
/// Number of integration steps
size_type N_m;
/// Is the angle step size
value_type dtheta_m;
/// Is the relativistic factor
value_type gamma_m;
/// Is the average radius
value_type ravg_m;
/// Is the phase
value_type phase_m;
/**
* Boolean which tells if a closed orbit for this configuration could be found (get set by the function findOrbit)
*/
bool converged_m;
/// Minimum energy needed in cyclotron
value_type Emin_m;
/// Maximum energy reached in cyclotron
value_type Emax_m;
/// Number of sectors (symmetry)
size_type nSector_m;
/**
* Stores the last orbit value (since we have to return to the beginning to check the convergence in the
* findOrbit() function. This last value is then deleted from the array but is stored in lastOrbitVal_m to
* compute the vertical oscillations)
*/
value_type lastOrbitVal_m;
/**
* Stores the last momentum value (since we have to return to the beginning to check the convergence in the
* findOrbit() function. This last value is then deleted from the array but is stored in lastMomentumVal_m to
* compute the vertical oscillations)
*/
value_type lastMomentumVal_m;
/**
* Boolean which is true if computeVerticalOscillations() executed, otherwise false. This is used for checking in
* getTunes() and getFrequencyError().
*/
bool vertOscDone_m;
/**
* Boolean which is true by default. "true": orbit integration over one sector only, "false": integration
* over 2*pi
*/
bool domain_m;
/// Defines the stepper for integration of the ODE's
Stepper stepper_m;
/// a guesss for the clo finder
value_type rguess_m;
/*!
* This quantity is defined in the paper "Transverse-Longitudinal Coupling by Space Charge in Cyclotrons"
* of Dr. Christian Baumgarten (2012)
* The lambda function takes the orbital frequency \f$ \omega_{o} \f$ (also defined in paper) as input argument.
*/
std::function<double(double)> acon_m = [](double wo) { return Physics::c / wo; };
/// Cyclotron unit \f$ \left[T\right] \f$ (Tesla)
/*!
* The lambda function takes the orbital frequency \f$ \omega_{o} \f$ as input argument.
*/
std::function<double(double, double)> bcon_m = [](double e0, double wo) {
return e0 * 1.0e7 / (/* physics::q0 */ 1.0 * Physics::c * Physics::c / wo);
};
MagneticField bField_m;
};
// -----------------------------------------------------------------------------------------------------------------------
// PUBLIC MEMBER FUNCTIONS
// -----------------------------------------------------------------------------------------------------------------------
template<typename Value_type, typename Size_type, class Stepper>
ClosedOrbitFinder<Value_type,
Size_type,
Stepper>::ClosedOrbitFinder(value_type E, value_type E0,
value_type wo, size_type N,
value_type accuracy, size_type maxit,
value_type Emin, value_type Emax,
size_type nSector, const std::string& fmapfn,
value_type rguess, const std::string& type,
value_type scaleFactor, bool domain)
: nxcross_m(0), nzcross_m(0), E_m(E), E0_m(E0), wo_m(wo), N_m(N), dtheta_m(Physics::two_pi/value_type(N)),
gamma_m(E/E0+1.0), ravg_m(0), phase_m(0), converged_m(false), Emin_m(Emin), Emax_m(Emax), nSector_m(nSector),
lastOrbitVal_m(0.0), lastMomentumVal_m(0.0),
vertOscDone_m(false), domain_m(domain), stepper_m(), rguess_m(rguess), bField_m(fmapfn, nSector)
{
if ( Emin_m > Emax_m )
throw OpalException("ClosedOrbitFinder::ClosedOrbitFinder()",
"Incorrect cyclotron energy (MeV) bounds: Maximum cyclotron energy smaller than minimum cyclotron energy.");
// // Even if the numbers are equal --> if statement is true.
// if ( E_m < Emin_m )
// throw OpalException("ClosedOrbitFinder::ClosedOrbitFinder()", "Kinetic energy smaller than minimum cyclotron energy");
if ( E_m > Emax_m )
throw OpalException("ClosedOrbitFinder::ClosedOrbitFinder()", "Kinetic energy exceeds cyclotron energy");
// velocity: beta = v/c = sqrt(1-1/(gamma*gamma))
if (gamma_m == 0)
throw OpalException("ClosedOrbitFinder::ClosedOrbitFinder()", "Relativistic factor equal zero.");
// if domain_m = true --> integrate over a single sector
if (domain_m) {
N_m /= nSector_m;
}
// reserve storage for the orbit and momentum (--> size = 0, capacity = N_m+1)
/*
* we need N+1 storage, since dtheta = 2pi/N (and not 2pi/(N-1)) that's why we need N+1 integration steps
* to return to the origin (but the return size is N_m)
*/
r_m.reserve(N_m + 1);
pr_m.reserve(N_m + 1);
// reserve memory of N_m for the properties (--> size = 0, capacity = N_m)
h_m.reserve(N_m);
ds_m.reserve(N_m);
fidx_m.reserve(N_m);
// read in magnetic fieldmap
bField_m.read(type, scaleFactor);
// compute closed orbit
converged_m = findOrbit(accuracy, maxit);
// compute h, ds, fidx, rav (average radius)
computeOrbitProperties();
}
template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type&
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getInverseBendingRadius()
{
return h_m;
}
template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type&
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getPathLength()
{
return ds_m;
}
template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type&
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getFieldIndex()
{
return fidx_m;
}
template<typename Value_type, typename Size_type, class Stepper>
std::pair<Value_type,Value_type> ClosedOrbitFinder<Value_type, Size_type, Stepper>::getTunes() {
// compute radial tune
value_type nur = computeTune(x_m,px_m[1],nxcross_m);
// compute nzcross_m
if (!vertOscDone_m)
computeVerticalOscillations();
// compute vertical tune
value_type nuz = computeTune(z_m,pz_m[1],nzcross_m);
return std::make_pair(nur,nuz);
}
template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getOrbit(value_type angle)
{
container_type r = r_m;
if (angle != 0.0) {
// compute the number of steps per degree
value_type deg_step = N_m / 360.0;
// compute starting point
size_type start = deg_step * angle;
// copy end to start
std::copy(r_m.begin() + start, r_m.end(), r.begin());
// copy start to end
std::copy_n(r_m.begin(), start, r.end() - start);
}
return r;
}
template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getMomentum(value_type angle)
{
container_type pr = pr_m;
if (angle != 0.0) {
// compute the number of steps per degree
value_type deg_step = N_m / 360.0;
// compute starting point
size_type start = deg_step * angle;
// copy end to start
std::copy(pr_m.begin() + start, pr_m.end(), pr.begin());
// copy start to end
std::copy_n(pr_m.begin(), start, pr.end() - start);
}
// change units from meters to \beta * \gamma
/* in Gordon paper:
*
* p = \gamma * \beta * a
*
* where a = c / \omega_{0} with \omega_{0} = 2 * \pi * \nu_{0} = 2 * \pi * \nu_{rf} / h
*
* c: speed of light
* h: harmonic number
* v_{rf}: nomial rf frequency
*
* Units:
*
* [a] = m --> [p] = m
*
* The momentum in \beta * \gamma is obtained by dividing by "a"
*/
value_type factor = 1.0 / acon_m(wo_m);
std::for_each(pr.begin(), pr.end(), [factor](value_type p) { return p * factor; });
return pr;
}
template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::value_type
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getGamma()
{
return gamma_m;
}
template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::value_type
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getAverageRadius()
{
return ravg_m;
}
template<typename Value_type, typename Size_type, class Stepper>
typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::value_type
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getFrequencyError()
{
// if the vertical oscillations aren't computed, we have to, since there we also compuote the frequency error.
if(!vertOscDone_m)
computeVerticalOscillations();
return phase_m;
}
template<typename Value_type, typename Size_type, class Stepper>
inline bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::isConverged() {
return converged_m;
}
// -----------------------------------------------------------------------------------------------------------------------
// PRIVATE MEMBER FUNCTIONS
// -----------------------------------------------------------------------------------------------------------------------
template<typename Value_type, typename Size_type, class Stepper>
bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type accuracy, size_type maxit) {
/* REMARK TO GORDON
* q' = 1/b = 1/bcon
* a' = a = acon
*/
// READ IN MAGNETIC FIELD: ONLY FOR STAND-ALONE PROGRAM
value_type bint, brint, btint;
// resize vectors (--> size = N_m+1, capacity = N_m+1), note: we do N_m+1 integration steps
r_m.resize(N_m+1);
pr_m.resize(N_m+1);
// store acon and bcon locally
value_type acon = acon_m(wo_m); // [acon] = m
value_type invbcon = 1.0 / bcon_m(E0_m, wo_m); // [bcon] = MeV*s/(C*m^2) = 10^6 T = 10^7 kG (kilo Gauss)
// helper constants
value_type p2; // p^2 = p*p
value_type pr2; // squared radial momentum (pr^2 = pr*pr)
value_type ptheta, invptheta; // Gordon, formula (5c)
value_type invdenom; // denominator for computing dr,dpr
value_type xold = 0.0; // for counting nxcross
// index for reaching next element of the arrays r and pr (no nicer way found yet)
size_type idx = 0;
// observer for storing the current value after each ODE step (e.g. Runge-Kutta step) into the containers of r and pr
auto store = [&](state_type& y, const value_type t)
{
r_m[idx] = y[0];
pr_m[idx] = y[1];
// count number of crossings (excluding starting point --> idx>0)
nxcross_m += (idx > 0) * (y[4] * xold < 0);
xold = y[4];
++idx;
};
// define the six ODEs (using lambda function)
std::function<void(const state_type&, state_type&, const double)> orbit_integration = [&](const state_type &y,
state_type &dydt,
const double theta)
{
pr2 = y[1] * y[1];
if (p2 < pr2)
throw OpalException("ClosedOrbitFinder::findOrbit()", "p_{r}^2 > p^{2} (defined in Gordon paper) --> Square root of negative number.");
// Gordon, formula (5c)
ptheta = std::sqrt(p2 - pr2);
invptheta = 1.0 / ptheta;
// interpolate values of magnetic field
bField_m.interpolate(bint, brint, btint, y[0], theta * 180.0 / Physics::pi);
bint *= invbcon;
brint *= invbcon;
// Gordon, formula (5a)
dydt[0] = y[0] * y[1] * invptheta;
// Gordon, formula (5b)
dydt[1] = ptheta - y[0] * bint;
// Gordon, formulas (9a) and (9b)
for (size_type i = 2; i < 5; i += 2) {
dydt[i] = (y[1] * y[i] + y[0] * p2 * y[i+1] * invptheta * invptheta) * invptheta;
dydt[i+1] = - y[1] * y[i+1] * invptheta - (bint + y[0] * brint) * y[i];
}
};
// define initial state container for integration: y = {r, pr, x1, px1, x2, px2}
state_type y(6);
// difference of last and first value of r (1. element) and pr (2. element)
container_type err(2);
// correction term for initial values: r = r + dr, pr = pr + dpr; Gordon, formula (17)
container_type delta = {0.0, 0.0};
// amplitude of error; Gordon, formula (18) (a = a')
value_type error = std::numeric_limits<value_type>::max();
// if niterations > maxit --> stop iteration
size_type niterations = 0;
/*
* Christian:
* N = 1440 ---> N = 720 ---> dtheta = 2PI/720 --> nsteps = 721
*
* 0, 2, 4, ... ---> jeden zweiten berechnene: 1, 3, 5, ... interpolieren --> 1440 Werte
*
* Matthias:
* N = 1440 --> dtheta = 2PI/1440 --> nsteps = 1441
*
* 0, 1, 2, 3, 4, 5, ... --> 1440 Werte
*
*/
// step size of energy
value_type dE;
if (Emin_m == Emax_m)
dE = 0.0;
else
dE = (E_m - Emin_m) / (Emax_m - Emin_m);
// iterate until suggested energy (start with minimum energy)
value_type E = Emin_m;
// energy increase not more than 0.25
dE = (dE > 0.25) ? 0.25 : dE;
// energy dependent values
value_type en = E / E0_m; // en = E/E0 = E/(mc^2) (E0 is potential energy)
value_type p = acon * std::sqrt(en * (2.0 + en)); // momentum [p] = m; Gordon, formula (3)
value_type gamma2 = (1.0 + en) * (1.0 + en); // = gamma^2
value_type beta = std::sqrt(1.0 - 1.0 / gamma2);
p2 = p * p; // p^2 = p*p
value_type invgamma4 = 1.0 / (gamma2 * gamma2); // = 1/gamma^4
// set initial values for radius and radial momentum for lowest energy Emin
// orbit, [r] = m; Gordon, formula (20)
// radial momentum; Gordon, formula (20)
container_type init;
if (rguess_m < 0)
init = {beta * acon, 0.0};