Commit 4f45cb63 authored by snuverink_j's avatar snuverink_j

add attribute to average field for closed orbit finder

parent a9204e3d
......@@ -60,9 +60,11 @@ class ClosedOrbitFinder
* @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.
* @param Nsectors is an int (default: 1). Number of sectors that the field map is averaged over
* in order to avoid first harmonic. Only valid if domain is false
*/
ClosedOrbitFinder(value_type E0, size_type N,
Cyclotron* cycl, bool domain = true);
Cyclotron* cycl, bool domain = true, int Nsectors = 1);
/// Returns the inverse bending radius (size of container N+1)
container_type getInverseBendingRadius(const value_type& angle = 0);
......@@ -207,6 +209,11 @@ class ClosedOrbitFinder
* over 2*pi
*/
bool domain_m;
/**
* Number of sectors to average the field map over
* in order to avoid first harmonic. Only valid if domain is false
*/
int nSectors_m;
/// Defines the stepper for integration of the ODE's
Stepper stepper_m;
......@@ -238,7 +245,7 @@ ClosedOrbitFinder<Value_type,
Size_type,
Stepper>::ClosedOrbitFinder(value_type E0,
size_type N, Cyclotron* cycl,
bool domain)
bool domain, int nSectors)
: nxcross_m(0)
, nzcross_m(0)
, E0_m(E0)
......@@ -250,6 +257,7 @@ ClosedOrbitFinder<Value_type,
/* , lastOrbitVal_m(0.0) */
/* , lastMomentumVal_m(0.0) */
, domain_m(domain)
, nSectors_m(nSectors)
, stepper_m()
, cycl_m(cycl)
{
......@@ -605,15 +613,27 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbitOfEnergy_m(
{
pr2 = y[1] * y[1];
if (p2 < pr2)
throw OpalException("ClosedOrbitFinder::findOrbit()",
throw OpalException("ClosedOrbitFinder::findOrbitOfEnergy_m()",
"p_{r}^2 > p^{2} (defined in Gordon paper) --> Square root of negative number.");
// Gordon, formula (5c)
ptheta = std::sqrt(p2 - pr2);
ptheta = std::sqrt(p2 - pr2);
invptheta = 1.0 / ptheta;
// interpolate values of magnetic field
cycl_m->apply(y[0], y[6], theta, brint, btint, bint);
// average field over the number of sectors
brint=0.0, btint=0.0, bint=0.0;
for (int i = 0; i<nSectors_m; i++) {
double angle = theta + i * Physics::two_pi / nSectors_m;
double tmpbr, tmpbt, tmpb;
// interpolate values of magnetic field
cycl_m->apply(y[0], y[6], angle, tmpbr, tmpbt, tmpb);
brint += tmpbr;
btint += tmpbt;
bint += tmpb;
}
brint /= nSectors_m;
btint /= nSectors_m;
bint /= nSectors_m;
// multiply by 1 / b
bint *= invbcon;
brint *= invbcon;
btint *= invbcon;
......
......@@ -50,6 +50,7 @@
#include <sys/time.h>
#include <boost/regex.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
extern Inform *gmsg;
......@@ -1250,7 +1251,6 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
/*
ToDo:
- add flag in order to calculate tunes from FMLOWE to FMHIGHE
- eliminate physics and error
*/
......@@ -1281,38 +1281,48 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
*/
Cyclotron* CyclotronElement = const_cast<Cyclotron*>(CyclotronVisitor.front());
bool full = !Attributes::getBool(itsAttr[Attrib::Distribution::SECTOR]);
int Nint = (int)Attributes::getReal(itsAttr[Attrib::Distribution::NSTEPS]);
bool full = !Attributes::getBool(itsAttr[Attrib::Distribution::SECTOR]);
int Nint = (int)Attributes::getReal(itsAttr[Attrib::Distribution::NSTEPS]);
int Nsectors = (int)Attributes::getReal(itsAttr[Attrib::Distribution::NSECTORS]);
if ( Nint < 0 )
throw OpalException("Distribution::CreateMatchedGaussDistribution()",
"Negative number of integration steps");
if ( Nsectors < 0 )
throw OpalException("Distribution::CreateMatchedGaussDistribution()",
"Negative number of sectors");
if ( Nsectors > 1 && full == false )
throw OpalException("Distribution::CreateMatchedGaussDistribution()",
"Averaging over sectors can only be done with SECTOR=FALSE");
*gmsg << "* ----------------------------------------------------" << endl;
*gmsg << "* About to find closed orbit and matched distribution " << endl;
*gmsg << "* I= " << I_m*1E3 << " (mA) E= " << E_m*1E-6 << " (MeV)" << endl;
*gmsg << "* EX= " << Attributes::getReal(itsAttr[Attrib::Distribution::EX])
<< " EY= " << Attributes::getReal(itsAttr[Attrib::Distribution::EY])
<< " ET= " << Attributes::getReal(itsAttr[Attrib::Distribution::ET]) << endl;
if ( full )
*gmsg << "* SECTOR: " << "match using all sectors" << endl;
if ( full ) {
*gmsg << "* SECTOR: " << "match using all sectors, with" << endl
<< "* NSECTORS = " << Nsectors << " to average the field over" << endl;
}
else
*gmsg << "* SECTOR: " << "match using single sector" << endl;
*gmsg << "* NSTEPS = " << Nint << endl
<< "* HN= " << CyclotronElement->getCyclHarm()
<< " PHIINIT= " << CyclotronElement->getPHIinit() << endl
<< "* ----------------------------------------------------" << endl;
if ( CyclotronElement->getFMLowE() < 0 ||
if ( CyclotronElement->getFMLowE() < 0 ||
CyclotronElement->getFMHighE() < 0 )
{
throw OpalException("Distribution::CreateMatchedGaussDistribution()",
"Missing attributes 'FMLOWE' and/or 'FMHIHGE' in "
"Missing attributes 'FMLOWE' and/or 'FMHIGHE' in "
"'CYCLOTRON' definition.");
}
std::size_t maxitCOF =
Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSCO]);
......@@ -1335,13 +1345,9 @@ 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(massIneV*1E-6, Nint, CyclotronElement, false);
cof_t cof(massIneV*1E-6, Nint, CyclotronElement, false, Nsectors);
cof.findOrbit(accuracy, maxitCOF, E_m*1E-6, denergy, rguess, true);
if ( !cof.findOrbit(accuracy, maxitCOF, E_m*1E-6, denergy, rguess, true) ) {
throw OpalException("Distribution::CreateMatchedGaussDistribution()",
"Closed orbit finder didn't converge.");
}
std::exit(0);
}
......@@ -3508,10 +3514,12 @@ 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::NSECTORS]
= Attributes::makeReal("NSECTORS", "Number of sectors to average field, for closed orbit finder ", 1);
itsAttr[Attrib::Distribution::FNAME]
= Attributes::makeString("FNAME", "File for reading in 6D particle "
......
......@@ -187,15 +187,16 @@ namespace Attrib
EY,
ET,
SECTOR, // Matched-Gauss: single sector or full machine
NSECTORS, // Matched-Gauss: number of sectors to average field
NSTEPS, // Matched-Gauss: number of steps for closed orbit finder
RGUESS, // Matched-Gauss: guess for closed orbit finder
DENERGY, // Matched-Gauss: energy step for closed orbit finder
LINE,
RESIDUUM,
MAXSTEPSCO,
MAXSTEPSSI,
ORDERMAPS,
// E2,
RGUESS,
DENERGY,
ID1, // special particle that the user can set
ID2, // special particle that the user can set
SCALABLE,
......
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