Commit 761d98c6 authored by frey_m's avatar frey_m
Browse files

Closed orbit fix + GeV input not MeV

parent 25478b1e
...@@ -442,15 +442,12 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc ...@@ -442,15 +442,12 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
* *
*/ */
value_type E; // starting energy value_type E = cycl_m->getFMLowE() * 1.0e3; // starting energy (GeV --> MeV)
value_type E_fin; // final energy value_type E_fin = ekin; // final energy
const value_type eps = dE * 1.0e-1; // articial constant for stopping criteria
if ( isTuneMode ) { if (isTuneMode) {
E = cycl_m->getFMLowE(); E_fin = cycl_m->getFMHighE() * 1.0e3; // (GeV --> MeV)
E_fin = cycl_m->getFMHighE();
} else {
E = ekin;
E_fin = ekin;
} }
namespace fs = boost::filesystem; namespace fs = boost::filesystem;
...@@ -474,12 +471,16 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc ...@@ -474,12 +471,16 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
container_type init; container_type init;
enum Guess {NONE, FIRST, SECOND}; enum Guess {NONE, FIRST, SECOND};
Guess guess = NONE; Guess guess = NONE;
value_type rn1, pn1; // normalised r, pr value of previous closed orbit value_type rn1 = 0.0, pn1 = 0.0; // normalised r, pr value of previous closed orbit
value_type rn2, pn2; // normalised r, pr value of second to previous closed orbit value_type rn2 = 0.0, pn2 = 0.0; // normalised r, pr value of second to previous closed orbit
// iterate until suggested energy (start with minimum energy) // iterate until suggested energy (start with minimum energy)
// increase energy by dE // increase energy by dE
for (; E <= E_fin ; E+=dE) { *gmsg << level3 << "Start iteration to find closed orbit of energy " << E_fin << " MeV "
<< "in steps of " << dE << " MeV." << endl;
for (; E <= E_fin + eps; E += dE) {
error = std::numeric_limits<value_type>::max(); error = std::numeric_limits<value_type>::max();
// energy dependent values // energy dependent values
...@@ -521,11 +522,16 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc ...@@ -521,11 +522,16 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
vz_m[0] = init[2]; vz_m[0] = init[2];
vpz_m[0] = init[3]; vpz_m[0] = init[3];
*gmsg << level3 << " Try to find orbit for " << E << " MeV ... ";
if ( !this->findOrbitOfEnergy_m(E, init, error, accuracy, maxit) ) { if ( !this->findOrbitOfEnergy_m(E, init, error, accuracy, maxit) ) {
*gmsg << "ClosedOrbitFinder didn't converge for energy " + std::to_string(E) + " MeV." << endl; *gmsg << endl << "ClosedOrbitFinder didn't converge for energy " + std::to_string(E) + " MeV." << endl;
guess = NONE; guess = NONE;
continue; continue;
} }
*gmsg << level3 << "Successfully found." << endl;
// store for next initial guess // store for next initial guess
rn2 = rn1; rn2 = rn1;
pn2 = pn1; pn2 = pn1;
...@@ -564,6 +570,8 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc ...@@ -564,6 +570,8 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
} }
} }
*gmsg << level3 << "Finished closed orbit finder successfully." << endl;
/* store last entry, since it is needed in computeVerticalOscillations(), because we have to do the same /* store last entry, since it is needed in computeVerticalOscillations(), because we have to do the same
* number of integrations steps there. * number of integrations steps there.
*/ */
......
...@@ -354,8 +354,8 @@ SigmaGenerator<Value_type, Size_type>::SigmaGenerator(value_type I, ...@@ -354,8 +354,8 @@ SigmaGenerator<Value_type, Size_type>::SigmaGenerator(value_type I,
, m_m(m) , m_m(m)
, niterations_m(0) , niterations_m(0)
, converged_m(false) , converged_m(false)
, Emin_m(cycl->getFMLowE()) , Emin_m(cycl->getFMLowE() * 1.0e3)
, Emax_m(cycl->getFMHighE()) , Emax_m(cycl->getFMHighE() * 1.0e3)
, N_m(N) , N_m(N)
, nSectors_m(Nsectors) , nSectors_m(Nsectors)
, nStepsPerSector_m(N/cycl->getSymmetry()) , nStepsPerSector_m(N/cycl->getSymmetry())
......
...@@ -103,10 +103,10 @@ OpalCyclotron::OpalCyclotron(): ...@@ -103,10 +103,10 @@ OpalCyclotron::OpalCyclotron():
("GEOMETRY", "Boundary Geometry for the Cyclotron"); ("GEOMETRY", "Boundary Geometry for the Cyclotron");
itsAttr[FMLOWE] = Attributes::makeReal itsAttr[FMLOWE] = Attributes::makeReal
("FMLOWE", "Minimal Energy [MeV] the fieldmap can accept. Used in GAUSSMATCHED", -1.0); ("FMLOWE", "Minimal Energy [GeV] the fieldmap can accept. Used in GAUSSMATCHED", -1.0);
itsAttr[FMHIGHE] = Attributes::makeReal itsAttr[FMHIGHE] = Attributes::makeReal
("FMHIGHE","Maximal Energy [MeV] the fieldmap can accept. Used in GAUSSMATCHED", -1.0); ("FMHIGHE","Maximal Energy [GeV] the fieldmap can accept. Used in GAUSSMATCHED", -1.0);
itsAttr[SPIRAL] = Attributes::makeBool itsAttr[SPIRAL] = Attributes::makeBool
("SPIRAL","Flag whether or not this is a spiral inflector simulation", false); ("SPIRAL","Flag whether or not this is a spiral inflector simulation", false);
......
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