27#ifndef CLOSEDORBITFINDER_H
28#define CLOSEDORBITFINDER_H
50#include <boost/numeric/odeint/integrate/integrate_n_steps.hpp>
51#include <boost/filesystem.hpp>
55template<
typename Value_type,
typename Size_type,
class Stepper>
82 Cyclotron* cycl,
bool domain =
true,
int Nsectors = 1);
94 std::pair<value_type,value_type>
getTunes();
136 bool isTuneMode =
false);
161 std::array<value_type,2>
x_m;
165 std::array<value_type,2>
z_m;
261template<
typename Value_type,
typename Size_type,
class Stepper>
267 bool domain,
int nSectors)
280 , nSectors_m(nSectors)
286 throw OpalException(
"ClosedOrbitFinder::ClosedOrbitFinder()",
287 "Incorrect cyclotron energy (MeV) bounds: Maximum cyclotron energy smaller than minimum cyclotron energy.");
312template<
typename Value_type,
typename Size_type,
class Stepper>
317 return rotate(angle, h_m);
322template<
typename Value_type,
typename Size_type,
class Stepper>
327 return rotate(angle, ds_m);
332template<
typename Value_type,
typename Size_type,
class Stepper>
337 return rotate(angle, fidx_m);
341template<
typename Value_type,
typename Size_type,
class Stepper>
344 value_type nur = computeTune(x_m,px_m[1],nxcross_m);
347 value_type nuz = computeTune(z_m,pz_m[1],nzcross_m);
349 return std::make_pair(nur,nuz);
352template<
typename Value_type,
typename Size_type,
class Stepper>
357 return rotate(angle, r_m);
362template<
typename Value_type,
typename Size_type,
class Stepper>
369 pr = rotate(angle, pr);
394template<
typename Value_type,
typename Size_type,
class Stepper>
401template<
typename Value_type,
typename Size_type,
class Stepper>
414template<
typename Value_type,
typename Size_type,
class Stepper>
456 E_fin = cycl_m->getFMHighE();
459 namespace fs = boost::filesystem;
462 std::string tunefile = (dir /
"tunes.dat").
string();
465 std::ofstream out(tunefile, std::ios::out);
468 << std::setw(15) <<
"# energy[MeV]"
469 << std::setw(15) <<
"radius_ini[m]"
470 << std::setw(15) <<
"momentum_ini[Beta Gamma]"
471 << std::setw(15) <<
"radius_avg[m]"
472 << std::setw(15) <<
"nu_r"
473 << std::setw(15) <<
"nu_z"
478 enum Guess {
NONE, FIRST, SECOND};
485 *
gmsg <<
level3 <<
"Start iteration to find closed orbit of energy " << E_fin <<
" MeV "
486 <<
"in steps of " << dE <<
" MeV." <<
endl;
488 for (; E <= E_fin + eps; E += dE) {
504 init = {beta * acon, 0.0, 0.0, 1.0};
506 init[0] = rguess * 0.001;
509 }
else if (guess == FIRST) {
511 init[0] = (beta*acon) * rn1;
514 }
else if (guess == SECOND) {
516 init[0] = (beta*acon) * (rn1 + (rn1-rn2));
517 init[1] = p*(pn1 + (pn1-pn2));
520 std::fill( r_m.begin(), r_m.end(), 0);
521 std::fill( pr_m.begin(), pr_m.end(), 0);
522 std::fill( vz_m.begin(), vz_m.end(), 0);
523 std::fill(vpz_m.begin(), vpz_m.end(), 0);
531 *
gmsg <<
level3 <<
" Try to find orbit for " << E <<
" MeV ... ";
533 if ( !this->findOrbitOfEnergy_m(E, init, error, accuracy, maxit) ) {
534 *
gmsg <<
endl <<
"ClosedOrbitFinder didn't converge for energy " + std::to_string(E) +
" MeV." <<
endl;
544 rn1 = r_m[0] / (acon * beta);
549 this->computeOrbitProperties(E);
550 value_type reo = this->getOrbit( cycl_m->getPHIinit())[0];
551 value_type peo = this->getMomentum(cycl_m->getPHIinit())[0];
552 std::pair<value_type , value_type > tunes = this->getTunes();
555 <<
"* ----------------------------" <<
endl
556 <<
"* Closed orbit info (Gordon units):" <<
endl
558 <<
"* kinetic energy: " << std::setw(12) << E <<
" [MeV]" <<
endl
559 <<
"* average radius: " << std::setw(12) << ravg_m <<
" [m]" <<
endl
560 <<
"* initial radius: " << std::setw(12) << reo <<
" [m]" <<
endl
561 <<
"* initial momentum: " << std::setw(12) << peo <<
" [Beta Gamma]" <<
endl
562 <<
"* frequency error: " << phase_m <<
endl
563 <<
"* horizontal tune: " << tunes.first <<
endl
564 <<
"* vertical tune: " << tunes.second <<
endl
565 <<
"* ----------------------------" <<
endl <<
endl;
567 std::ofstream out(tunefile, std::ios::app);
569 << std::setw(15) << E
570 << std::setw(15) << reo
571 << std::setw(15) << peo
572 << std::setw(15) << ravg_m
573 << std::setw(15) << tunes.first
574 << std::setw(15) << tunes.second <<
std::endl;
579 *
gmsg <<
level3 <<
"Finished closed orbit finder successfully." <<
endl;
595 phase_m *= cycl_m->getSymmetry();
598 return error < accuracy;
601template<
typename Value_type,
typename Size_type,
class Stepper>
613 value_type invbcon = 1.0 / bcon_m(E0_m, wo_m);
629 nxcross_m += (idx > 0) * (y[4] * xold < 0);
633 nzcross_m += (idx > 0) * (y[10] * zold < 0);
656 value_type invgamma4 = 1.0 / (gamma2 * gamma2);
671 throw OpalException(
"ClosedOrbitFinder::findOrbitOfEnergy_m()",
672 "p_{r}^2 > p^{2} (defined in Gordon paper) --> Square root of negative number.");
676 invptheta = 1.0 / ptheta;
678 brint=0.0, btint=0.0, bint=0.0;
679 for (
int i = 0; i<nSectors_m; i++) {
681 double tmpbr, tmpbt, tmpb;
683 cycl_m->apply(y[0], y[6], angle, tmpbr, tmpbt, tmpb);
697 dydt[0] = y[0] * y[1] * invptheta;
699 dydt[1] = ptheta - y[0] * bint;
702 dydt[i] = (y[1] * y[i] + y[0] * p2 * y[i+1] * invptheta * invptheta) * invptheta;
703 dydt[i+1] = - y[1] * y[i+1] * invptheta - (bint + y[0] * brint) * y[i];
708 dydt[i] = y[0] * y[i+1] * invptheta;
709 dydt[i+1] = (y[0] * brint - y[1] * invptheta * btint) * y[i];
713 dydt[12] = y[0] * invptheta * gamma - 1;
734 y = {init[0],init[1],
735 x_m[0], px_m[0], x_m[1], px_m[1],
737 z_m[0], pz_m[0], z_m[1], pz_m[1],
743 boost::numeric::odeint::integrate_n_steps(stepper_m, orbit_integration,y,0.0,dtheta_m,N_m,store);
763 err[0] = r_m[N_m] - r_m[0];
764 err[1] = pr_m[N_m] - pr_m[0];
767 value_type invdenom = 1.0 / (x_m[0] + px_m[1] - 2.0);
768 delta[0] = ((px_m[1] - 1.0) * err[0] - x_m[1] * err[1]) * invdenom;
769 delta[1] = (( x_m[0] - 1.0) * err[1] - px_m[0] * err[0]) * invdenom;
776 error =
std::sqrt(delta[0] * delta[0] + delta[1] * delta[1] * invgamma4) / r_m[0];
780 }
while ((error > accuracy) && (niterations++ < maxit));
782 if (error > accuracy)
783 *
gmsg <<
"findOrbit not converged after " << niterations <<
" iterations with error: " << error <<
". Needed accuracy " << accuracy <<
endl;
785 return (error < accuracy);
788template<
typename Value_type,
typename Size_type,
class Stepper>
800 bool negative = std::signbit(y[1]);
802 bool uneven = (ncross % 2);
812 muPrime = -std::acosh(abscos);
828 if (!negative && std::signbit(
std::sin(mu))) {
830 }
else if (negative && !std::signbit(
std::sin(mu))) {
841 mu *= cycl_m->getSymmetry();
846template<
typename Value_type,
typename Size_type,
class Stepper>
858 value_type invbcon = 1.0 / bcon_m(E0_m, wo_m);
872 cycl_m->apply(r_m[i], vz_m[i], theta, brint, btint, bint);
882 throw OpalException(
"ClosedOrbitFinder::computeOrbitProperties()",
883 "p_{r}^2 > p^{2} " + std::to_string(p) +
" " + std::to_string(pr_m[i]) +
" (defined in Gordon paper) --> Square root of negative number.");
885 ptheta =
std::sqrt(p2 - pr_m[i] * pr_m[i]);
887 fidx_m[i] = (brint * ptheta - btint * pr_m[i] / r_m[i]) / p2;
890 ds_m[i] = std::hypot(r_m[i] * pr_m[i] / ptheta,r_m[i]) * dtheta_m;
897 ravg_m = std::accumulate(r_m.begin(),r_m.end(),0.0) /
value_type(r_m.size());
900template<
typename Value_type,
typename Size_type,
class Stepper>
910 angle = 360.0 + angle;
914 unsigned int start = deg_step * angle;
916 start %= orbitProperty.size();
919 std::copy(orbitProperty.begin() + start, orbitProperty.end(), orbitPropertyCopy.begin());
922 std::copy_n(orbitProperty.begin(), start, orbitPropertyCopy.end() - start);
924 return orbitPropertyCopy;
Tps< T > cos(const Tps< T > &x)
Cosine.
Tps< T > sin(const Tps< T > &x)
Sine.
Tps< T > sqrt(const Tps< T > &x)
Square root.
T * value_type(const SliceIterator< T > &)
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
bool for_each(const BareFieldIterator< T, D > &p, SameFieldID s, C)
PETE_TUTree< FnArcCos, typename T::PETE_Expr_t > acos(const PETE_Expr< T > &l)
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Inform & endl(Inform &inf)
Inform & level3(Inform &inf)
py::list function(PolynomialPatch *patch, py::list point)
constexpr double u_two_pi
The value of.
constexpr double two_pi
The value of.
constexpr double c
The velocity of light in m/s.
constexpr double pi
The value of.
std::string getInputBasename()
get input file name without extension
static OpalData * getInstance()
std::string getAuxiliaryOutputDirectory() const
get the name of the the additional data directory
void read(const double &scaleFactor)
virtual double getFMHighE() const
virtual double getBScale() const
virtual double getSymmetry() const
virtual double getFMLowE() const
std::vector< value_type > container_type
Type of container for storing quantities (path length, orbit, etc.)
std::pair< value_type, value_type > getTunes()
Returns the radial and vertical tunes (in that order)
std::array< value_type, 2 > px_m
Stores current momenta in horizontal direction for the solutions of the ODE with different initial va...
size_type N_m
Number of integration steps.
ClosedOrbitFinder(value_type E0, value_type q, size_type N, Cyclotron *cycl, bool domain=true, int Nsectors=1)
Sets the initial values for the integration and calls findOrbit().
container_type vz_m
Stores the vertical oribt (size: N_m+1)
size_type nzcross_m
Counts the number of zero-line crossings in vertical direction (used for computing vertical tune)
value_type phase_m
Is the phase.
bool findOrbitOfEnergy_m(const value_type &, container_type &, value_type &, const value_type &, size_type)
std::vector< value_type > state_type
Type for holding state of ODE values.
void computeOrbitProperties(const value_type &E)
Fills in the values of h_m, ds_m, fidx_m.
std::array< value_type, 2 > z_m
Stores current position in vertical direction for the solutions of the ODE with different initial val...
container_type fidx_m
Stores the field index.
value_type wo_m
Is the nominal orbital frequency.
container_type ds_m
Stores the step length.
bool findOrbit(value_type accuracy, size_type maxit, value_type ekin, value_type dE=1.0, value_type rguess=-1.0, bool isTuneMode=false)
Computes the closed orbit.
container_type h_m
Stores the inverse bending radius.
std::function< double(double, double)> bcon_m
Cyclotron unit (Tesla)
container_type r_m
Stores the radial orbit (size: N_m+1)
Size_type size_type
Type for specifying sizes.
container_type vpz_m
Stores the vertical momentum.
value_type getAverageRadius()
Returns the average orbit radius.
std::function< double(double)> acon_m
value_type dtheta_m
Is the angle step size.
std::function< void(const state_type &, state_type &, const double)> function_t
container_type getFieldIndex(const value_type &angle=0)
Returns the field index (size of container N+1)
Value_type value_type
Type of variables.
container_type getMomentum(value_type angle=0)
container_type getOrbit(value_type angle=0)
container_type rotate(value_type angle, const container_type &orbitProperty)
This function computes nzcross_ which is used to compute the tune in z-direction and the frequency er...
bool isConverged()
Returns true if a closed orbit could be found.
container_type getPathLength(const value_type &angle=0)
Returns the step lengths of the path (size of container N+1)
std::array< value_type, 2 > pz_m
Stores current momenta in vertical direction for the solutions of the ODE with different initial valu...
std::array< value_type, 2 > x_m
Stores current position in horizontal direction for the solutions of the ODE with different initial v...
value_type ravg_m
Is the average radius.
container_type pr_m
Stores the radial momentum.
container_type getInverseBendingRadius(const value_type &angle=0)
Returns the inverse bending radius (size of container N+1)
value_type getFrequencyError()
Returns the frequency error.
value_type computeTune(const std::array< value_type, 2 > &, value_type, size_type)
This function is called by the function getTunes().
value_type E0_m
Is the rest mass [MeV / c**2].
size_type nxcross_m
Counts the number of zero-line crossings in horizontal direction (used for computing horizontal tune)
Stepper stepper_m
Defines the stepper for integration of the ODE's.
The base class for all OPAL exceptions.
virtual const std::string & what() const
Return the message string for the exception.
virtual const std::string & where() const
Return the name of the method or function which detected the exception.