Commit db1cd250 authored by frey_m's avatar frey_m

Issue #273: Use cyclotron element as input argument

modified:   ClosedOrbitFinder.h
modified:   Distribution.cpp
modified:   SigmaGenerator.h
parent ac94b1f6
......@@ -52,26 +52,13 @@ class ClosedOrbitFinder
/*!
* @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 cycl is the cyclotron element
* @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);
ClosedOrbitFinder(value_type E, value_type E0, size_type N,
const Cyclotron* cycl, bool domain = true);
/// Returns the inverse bending radius (size of container N+1)
container_type getInverseBendingRadius(const value_type& angle = 0);
......@@ -116,17 +103,18 @@ class ClosedOrbitFinder
/// 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
* @param rguess initial radius guess in [mm]
*/
bool findOrbit(value_type, size_type);
bool findOrbit(value_type, size_type, value_type = -1.0);
/// Fills in the values of h_m, ds_m, fidx_m. It gets called by in by constructor.
/// Fills in the values of h_m, ds_m, fidx_m.
void computeOrbitProperties();
private:
/// 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)
......@@ -173,6 +161,9 @@ class ClosedOrbitFinder
value_type E0_m;
/// Is the nominal orbital frequency
/* (see paper of Dr. C. Baumgarten: "Transverse-Longitudinal
* Coupling by Space Charge in Cyclotrons" (2012), formula (1))
*/
value_type wo_m;
/// Number of integration steps
size_type N_m;
......@@ -188,11 +179,6 @@ class ClosedOrbitFinder
/// 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;
......@@ -230,9 +216,6 @@ class ClosedOrbitFinder
/// 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"
......@@ -260,32 +243,26 @@ 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)
size_type N, const Cyclotron* cycl,
bool domain)
: nxcross_m(0)
, nzcross_m(0)
, E_m(E)
, E0_m(E0)
, wo_m(wo)
, wo_m(cycl->getRfFrequ()*1E6/cycl->getCyclHarm()*2.0*Physics::pi)
, 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)
, Emin_m(cycl->getFMLowE())
, Emax_m(cycl->getFMHighE())
, nSector_m(cycl->getSymmetry())
, lastOrbitVal_m(0.0)
, lastMomentumVal_m(0.0)
, vertOscDone_m(false)
, domain_m(domain)
, stepper_m()
, rguess_m(rguess)
{
if ( Emin_m > Emax_m )
......@@ -322,16 +299,10 @@ ClosedOrbitFinder<Value_type,
fidx_m.reserve(N_m);
// read in magnetic fieldmap
bField_m.setFieldMapFN(fmapfn);
bField_m.setFieldMapFN(cycl->getFieldMapFN());
bField_m.setSymmetry(nSector_m);
int fieldflag = bField_m.getFieldFlag(type);
bField_m.read(fieldflag, scaleFactor);
// compute closed orbit
converged_m = findOrbit(accuracy, maxit);
// compute h, ds, fidx, rav (average radius)
computeOrbitProperties();
int fieldflag = bField_m.getFieldFlag(cycl->getCyclotronType());
bField_m.read(fieldflag, cycl->getBScale());
}
template<typename Value_type, typename Size_type, class Stepper>
......@@ -445,17 +416,15 @@ typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::value_type
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) {
bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type accuracy,
size_type maxit,
value_type rguess)
{
/* REMARK TO GORDON
* q' = 1/b = 1/bcon
* a' = a = acon
......@@ -576,10 +545,10 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
// radial momentum; Gordon, formula (20)
container_type init;
if (rguess_m < 0)
if (rguess < 0)
init = {beta * acon, 0.0};
else
init = {rguess_m/1000.0, 0.0};
init = {rguess * 0.001, 0.0};
// store initial values for updating values for higher energies
container_type previous_init = {0.0, 0.0};
......
......@@ -1311,12 +1311,9 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
<< " PHIINIT= " << CyclotronElement->getPHIinit() << endl
<< "* ----------------------------------------------------" << endl;
const double wo = CyclotronElement->getRfFrequ()*1E6/CyclotronElement->getCyclHarm()*2.0*Physics::pi;
const double fmLowE = CyclotronElement->getFMLowE();
const double fmHighE = CyclotronElement->getFMHighE();
if ( fmLowE < 0 || fmHighE < 0 ) {
if ( CyclotronElement->getFMLowE() < 0 ||
CyclotronElement->getFMHighE() < 0 )
{
throw OpalException("Distribution::CreateMatchedGaussDistribution()",
"Missing attributes 'FMLOWE' and/or 'FMHIHGE' in "
"'CYCLOTRON' definition.");
......@@ -1329,8 +1326,6 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
double rguess =
Attributes::getReal(itsAttr[Attrib::Distribution::RGUESS]);
int nSector = (int)CyclotronElement->getSymmetry();
double accuracy =
Attributes::getReal(itsAttr[Attrib::Distribution::RESIDUUM]);
......@@ -1340,11 +1335,14 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
typedef boost::numeric::odeint::runge_kutta4<container_t> rk4_t;
typedef ClosedOrbitFinder<double,unsigned int, rk4_t> cof_t;
cof_t cof(E_m*1E-6, massIneV*1E-6, wo, Nint, accuracy,
maxitCOF, fmLowE, fmHighE, nSector,
CyclotronElement->getFieldMapFN(), rguess,
CyclotronElement->getCyclotronType(),
CyclotronElement->getBScale(), false);
cof_t cof(E_m*1E-6, massIneV*1E-6, Nint, CyclotronElement, false);
if ( !cof.findOrbit(accuracy, maxitCOF, rguess) ) {
throw OpalException("Distribution::CreateMatchedGaussDistribution()",
"Closed orbit finder didn't converge.");
}
cof.computeOrbitProperties();
std::pair<double, double> tunes = cof.getTunes();
double ravg = cof.getAverageRadius(); // average radius
......@@ -1374,25 +1372,18 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
Attributes::getReal(itsAttr[Attrib::Distribution::EX])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::EY])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::ET])*1E6,
wo,
E_m*1E-6,
CyclotronElement->getCyclHarm(),
massIneV*1E-6,
fmLowE,
fmHighE,
nSector,
CyclotronElement,
Nint,
CyclotronElement->getFieldMapFN(),
Attributes::getReal(itsAttr[Attrib::Distribution::ORDERMAPS]),
CyclotronElement->getBScale(),
writeMap);
if (siggen->match(accuracy,
Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSSI]),
maxitCOF,
CyclotronElement->getPHIinit(),
CyclotronElement,
rguess,
CyclotronElement->getCyclotronType(),
false, full)) {
std::array<double,3> Emit = siggen->getEmittances();
......
......@@ -101,25 +101,17 @@ public:
* @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$
* @param nSector is the number of sectors (symmetry assumption)
* @param cycl is the cyclotron element
* @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)
* @param fieldmap is the location of the file that specifies the magnetic field
* @param truncOrder is the truncation order for power series of the Hamiltonian
* @param scaleFactor for the magnetic field (default: 1.0)
* @param write is a boolean (default: true). If true all maps of all iterations are stored, otherwise not.
*/
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);
value_type E, value_type m, const Cyclotron* cycl,
size_type N, size_type truncOrder, bool write = true);
/// Searches for a matched distribution.
/*!
......@@ -129,13 +121,13 @@ public:
* 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°)
* @param guess value of radius for closed orbit finder
* @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, value_type angle,
value_type guess, const std::string& type, bool harmonic, bool full);
bool match(value_type accuracy, size_type maxit, size_type maxitOrbit,
const Cyclotron* cycl, value_type rguess, bool harmonic, bool full);
/*!
* Eigenvalue / eigenvector solver
......@@ -336,39 +328,33 @@ SigmaGenerator<Value_type, Size_type>::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,
const Cyclotron* cycl,
size_type N,
const std::string& fieldmap,
size_type truncOrder,
value_type scaleFactor,
bool write)
: I_m(I)
, wo_m(wo)
, wo_m(cycl->getRfFrequ()*1E6/cycl->getCyclHarm()*2.0*Physics::pi)
, E_m(E)
, gamma_m(E/m+1.0)
, gamma2_m(gamma_m*gamma_m)
, nh_m(nh)
, nh_m(cycl->getCyclHarm())
, beta_m(std::sqrt(1.0-1.0/gamma2_m))
, m_m(m)
, niterations_m(0)
, converged_m(false)
, Emin_m(Emin)
, Emax_m(Emax)
, nSector_m(nSector)
, Emin_m(cycl->getFMLowE())
, Emax_m(cycl->getFMHighE())
, nSector_m(cycl->getSymmetry())
, N_m(N)
, nStepsPerSector_m(N/nSector)
, nStepsPerSector_m(N/cycl->getSymmetry())
, nSteps_m(N)
, error_m(std::numeric_limits<value_type>::max())
, fieldmap_m(fieldmap)
, fieldmap_m(cycl->getFieldMapFN())
, truncOrder_m(truncOrder)
, write_m(write)
, scaleFactor_m(scaleFactor)
, scaleFactor_m(cycl->getBScale())
, sigmas_m(nStepsPerSector_m)
, rinit_m(0.0)
, prinit_m(0.0)
......@@ -443,9 +429,8 @@ template<typename Value_type, typename Size_type>
bool SigmaGenerator<Value_type, Size_type>::match(value_type accuracy,
size_type maxit,
size_type maxitOrbit,
value_type angle,
const Cyclotron* cycl,
value_type rguess,
const std::string& type,
bool harmonic, bool full)
{
/* compute the equilibrium orbit for energy E_
......@@ -476,15 +461,20 @@ template<typename Value_type, typename Size_type>
if (!harmonic) {
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,
type, scaleFactor_m, false);
boost::numeric::odeint::runge_kutta4<container_type> > cof(E_m, m_m, N_m, cycl, false);
if ( !cof.findOrbit(accuracy, maxitOrbit, rguess) ) {
throw OpalException("SigmaGenerator::match()",
"Closed orbit finder didn't converge.");
}
cof.computeOrbitProperties();
// 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);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment