Commit 0c6740d7 authored by frey_m's avatar frey_m

COF: add 'denergy' as input argument

parent 39c89e34
......@@ -105,9 +105,10 @@ class ClosedOrbitFinder
/*!
* @param accuracy specifies the accuracy of the closed orbit
* @param maxit is the maximal number of iterations done for finding the closed orbit
* @param energy step increase [MeV]
* @param rguess initial radius guess in [mm]
*/
bool findOrbit(value_type, size_type, value_type = -1.0);
bool findOrbit(value_type, size_type, value_type = 1.0, value_type = -1.0);
/// Fills in the values of h_m, ds_m, fidx_m.
void computeOrbitProperties();
......@@ -395,6 +396,7 @@ typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::value_type
template<typename Value_type, typename Size_type, class Stepper>
bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type accuracy,
size_type maxit,
value_type dE,
value_type rguess)
{
/* REMARK TO GORDON
......@@ -514,20 +516,9 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
*
*/
// step size of energy
value_type dE;
if (cycl_m->getFMLowE() == cycl_m->getFMHighE())
dE = 0.0;
else
dE = (E_m - cycl_m->getFMLowE()) / (cycl_m->getFMHighE() - cycl_m->getFMLowE());
// iterate until suggested energy (start with minimum energy)
value_type E = cycl_m->getFMLowE();
// 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)
......
......@@ -1325,6 +1325,13 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
double rguess =
Attributes::getReal(itsAttr[Attrib::Distribution::RGUESS]);
double denergy = 1000.0 *
Attributes::getReal(itsAttr[Attrib::Distribution::DENERGY]);
if ( denergy < 0.0 )
throw OpalException("Distribution:CreateMatchedGaussDistribution()",
"DENERGY < 0");
double accuracy =
Attributes::getReal(itsAttr[Attrib::Distribution::RESIDUUM]);
......@@ -1337,7 +1344,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
cof_t cof(E_m*1E-6, massIneV*1E-6, Nint, CyclotronElement, false);
if ( !cof.findOrbit(accuracy, maxitCOF, rguess) ) {
if ( !cof.findOrbit(accuracy, maxitCOF, denergy, rguess) ) {
throw OpalException("Distribution::CreateMatchedGaussDistribution()",
"Closed orbit finder didn't converge.");
}
......@@ -1383,6 +1390,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSSI]),
maxitCOF,
CyclotronElement,
denergy,
rguess,
false, full)) {
......@@ -3565,6 +3573,9 @@ void Distribution::setAttributes() {
itsAttr[Attrib::Distribution::RGUESS]
= Attributes::makeReal("RGUESS", "Guess value of radius (m) for closed orbit finder ", -1);
itsAttr[Attrib::Distribution::DENERGY]
= Attributes::makeReal("DENERGY", "Energy step size for closed orbit finder [GeV]", 0.001);
itsAttr[Attrib::Distribution::FNAME]
......
......@@ -194,6 +194,7 @@ namespace Attrib
ORDERMAPS,
// E2,
RGUESS,
DENERGY,
ID1, // special particle that the user can set
ID2, // special particle that the user can set
SCALABLE,
......
......@@ -121,13 +121,14 @@ 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 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 rguess, bool harmonic, bool full);
Cyclotron* cycl, value_type denergy, value_type rguess, bool harmonic, bool full);
/*!
* Eigenvalue / eigenvector solver
......@@ -430,6 +431,7 @@ template<typename Value_type, typename Size_type>
size_type maxit,
size_type maxitOrbit,
Cyclotron* cycl,
value_type denergy,
value_type rguess,
bool harmonic, bool full)
{
......@@ -463,7 +465,7 @@ template<typename Value_type, typename Size_type>
ClosedOrbitFinder<value_type, size_type,
boost::numeric::odeint::runge_kutta4<container_type> > cof(E_m, m_m, N_m, cycl, false);
if ( !cof.findOrbit(accuracy, maxitOrbit, rguess) ) {
if ( !cof.findOrbit(accuracy, maxitOrbit, denergy, rguess) ) {
throw OpalException("SigmaGenerator::match()",
"Closed orbit finder didn't converge.");
}
......
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