Commit 6172e02c authored by snuverink_j's avatar snuverink_j
Browse files

Merge branch 'cof' into 'master'

Cof

Closes #282

See merge request !60
parents 34805a45 99ee2b64
This diff is collapsed.
......@@ -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,32 +1345,8 @@ 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, Nint, CyclotronElement, false);
if ( !cof.findOrbit(accuracy, maxitCOF, denergy, 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
container_t reo = cof.getOrbit(CyclotronElement->getPHIinit());
container_t peo = cof.getMomentum(CyclotronElement->getPHIinit());
*gmsg << "* ----------------------------" << endl
<< "* Closed orbit info (Gordon units):" << endl
<< "*" << endl
<< "* average radius: " << ravg << " [m]" << endl
<< "* initial radius: " << reo[0] << " [m]" << endl
<< "* initial momentum: " << peo[0] << " [Beta Gamma]" << endl
<< "* frequency error: " << cof.getFrequencyError() << endl
<< "* horizontal tune: " << tunes.first << endl
<< "* vertical tune: " << tunes.second << endl
<< "* ----------------------------" << endl << endl;
cof_t cof(massIneV*1E-6, Nint, CyclotronElement, false, Nsectors);
cof.findOrbit(accuracy, maxitCOF, E_m*1E-6, denergy, rguess, true);
std::exit(0);
}
......@@ -3528,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,
......
......@@ -463,14 +463,14 @@ 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, N_m, cycl, false);
boost::numeric::odeint::runge_kutta4<container_type> > cof(m_m, N_m, cycl, false);
if ( !cof.findOrbit(accuracy, maxitOrbit, denergy, rguess) ) {
if ( !cof.findOrbit(accuracy, maxitOrbit, E_m, denergy, rguess) ) {
throw OpalException("SigmaGenerator::match()",
"Closed orbit finder didn't converge.");
}
cof.computeOrbitProperties();
cof.computeOrbitProperties(E_m);
// properties of one turn
tunes = cof.getTunes();
......@@ -544,8 +544,8 @@ template<typename Value_type, typename Size_type>
Mcycs[i] = mapgen.generateMap(H_m(h[i],
h[i]*h[i]+fidx[i],
-fidx[i]),
ds[i],truncOrder_m);
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)),
......
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