Commit 773f6668 authored by frey_m's avatar frey_m

Resolve "ClosedOrbitFinder/SigmaGenerator: Charge not applied correctly"

parent 37da44f0
...@@ -70,6 +70,7 @@ class ClosedOrbitFinder ...@@ -70,6 +70,7 @@ class ClosedOrbitFinder
/// Sets the initial values for the integration and calls findOrbit(). /// Sets the initial values for the integration and calls findOrbit().
/*! /*!
* @param E0 is the potential energy (particle energy at rest) [MeV]. * @param E0 is the potential energy (particle energy at rest) [MeV].
* @param q is the particle charge [e]
* @param N specifies the number of splits (2pi/N), i.e number of integration steps * @param N specifies the number of splits (2pi/N), i.e number of integration steps
* @param cycl is the cyclotron element * @param cycl is the cyclotron element
* @param domain is a boolean (default: true). If "true" the closed orbit is computed over a single sector, * @param domain is a boolean (default: true). If "true" the closed orbit is computed over a single sector,
...@@ -77,7 +78,7 @@ class ClosedOrbitFinder ...@@ -77,7 +78,7 @@ class ClosedOrbitFinder
* @param Nsectors is an int (default: 1). Number of sectors that the field map is averaged over * @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 * in order to avoid first harmonic. Only valid if domain is false
*/ */
ClosedOrbitFinder(value_type E0, size_type N, ClosedOrbitFinder(value_type E0, value_type q, size_type N,
Cyclotron* cycl, bool domain = true, int Nsectors = 1); Cyclotron* cycl, bool domain = true, int Nsectors = 1);
/// Returns the inverse bending radius (size of container N+1) /// Returns the inverse bending radius (size of container N+1)
...@@ -154,7 +155,7 @@ class ClosedOrbitFinder ...@@ -154,7 +155,7 @@ class ClosedOrbitFinder
// void computeVerticalOscillations(); // void computeVerticalOscillations();
/// This function rotates the calculated closed orbit finder properties to the initial angle /// This function rotates the calculated closed orbit finder properties to the initial angle
container_type rotate(value_type angle, container_type& orbitProperty); container_type rotate(value_type angle, const container_type& orbitProperty);
/// Stores current position in horizontal direction for the solutions of the ODE with different initial values /// Stores current position in horizontal direction for the solutions of the ODE with different initial values
std::array<value_type,2> x_m; // x_m = [x1, x2] std::array<value_type,2> x_m; // x_m = [x1, x2]
...@@ -188,6 +189,9 @@ class ClosedOrbitFinder ...@@ -188,6 +189,9 @@ class ClosedOrbitFinder
/// Is the rest mass [MeV / c**2] /// Is the rest mass [MeV / c**2]
value_type E0_m; value_type E0_m;
// Is the particle charge [e]
value_type q_m;
/// Is the nominal orbital frequency /// Is the nominal orbital frequency
/* (see paper of Dr. C. Baumgarten: "Transverse-Longitudinal /* (see paper of Dr. C. Baumgarten: "Transverse-Longitudinal
* Coupling by Space Charge in Cyclotrons" (2012), formula (1)) * Coupling by Space Charge in Cyclotrons" (2012), formula (1))
...@@ -243,8 +247,8 @@ class ClosedOrbitFinder ...@@ -243,8 +247,8 @@ class ClosedOrbitFinder
/*! /*!
* The lambda function takes the orbital frequency \f$ \omega_{o} \f$ as input argument. * The lambda function takes the orbital frequency \f$ \omega_{o} \f$ as input argument.
*/ */
std::function<double(double, double)> bcon_m = [](double e0, double wo) { std::function<double(double, double)> bcon_m = [this](double e0, double wo) {
return e0 * 1.0e7 / (/* physics::q0 */ 1.0 * Physics::c * Physics::c / wo); return e0 * 1.0e7 / (q_m * Physics::c * Physics::c / wo);
}; };
Cyclotron* cycl_m; Cyclotron* cycl_m;
...@@ -258,11 +262,13 @@ template<typename Value_type, typename Size_type, class Stepper> ...@@ -258,11 +262,13 @@ template<typename Value_type, typename Size_type, class Stepper>
ClosedOrbitFinder<Value_type, ClosedOrbitFinder<Value_type,
Size_type, Size_type,
Stepper>::ClosedOrbitFinder(value_type E0, Stepper>::ClosedOrbitFinder(value_type E0,
value_type q,
size_type N, Cyclotron* cycl, size_type N, Cyclotron* cycl,
bool domain, int nSectors) bool domain, int nSectors)
: nxcross_m(0) : nxcross_m(0)
, nzcross_m(0) , nzcross_m(0)
, E0_m(E0) , E0_m(E0)
, q_m(q)
, wo_m(cycl->getRfFrequ(0)*1E6/cycl->getCyclHarm()*2.0*Physics::pi) , wo_m(cycl->getRfFrequ(0)*1E6/cycl->getCyclHarm()*2.0*Physics::pi)
, N_m(N) , N_m(N)
, dtheta_m(Physics::two_pi/value_type(N)) , dtheta_m(Physics::two_pi/value_type(N))
...@@ -894,15 +900,21 @@ void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeOrbitProperties(c ...@@ -894,15 +900,21 @@ void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeOrbitProperties(c
template<typename Value_type, typename Size_type, class Stepper> template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type
ClosedOrbitFinder<Value_type, Size_type, Stepper>::rotate(value_type angle, container_type &orbitProperty) { ClosedOrbitFinder<Value_type, Size_type, Stepper>::rotate(value_type angle, const container_type &orbitProperty) {
container_type orbitPropertyCopy = orbitProperty; container_type orbitPropertyCopy = orbitProperty;
// compute the number of steps per degree // compute the number of steps per degree
value_type deg_step = N_m / 360.0; value_type deg_step = N_m / 360.0;
if (angle < 0.0) {
angle = 360.0 + angle;
}
// compute starting point // compute starting point
size_type start = deg_step * angle; unsigned int start = deg_step * angle;
start %= orbitProperty.size();
// copy end to start // copy end to start
std::copy(orbitProperty.begin() + start, orbitProperty.end(), orbitPropertyCopy.begin()); std::copy(orbitProperty.begin() + start, orbitProperty.end(), orbitPropertyCopy.begin());
......
...@@ -248,7 +248,7 @@ void Distribution::execute() { ...@@ -248,7 +248,7 @@ void Distribution::execute() {
void Distribution::update() { void Distribution::update() {
} }
void Distribution::create(size_t &numberOfParticles, double massIneV) { void Distribution::create(size_t &numberOfParticles, double massIneV, double charge) {
/* /*
* If Options::cZero is true than we reflect generated distribution * If Options::cZero is true than we reflect generated distribution
...@@ -284,7 +284,7 @@ void Distribution::create(size_t &numberOfParticles, double massIneV) { ...@@ -284,7 +284,7 @@ void Distribution::create(size_t &numberOfParticles, double massIneV) {
switch (distrTypeT_m) { switch (distrTypeT_m) {
case DistrTypeT::MATCHEDGAUSS: case DistrTypeT::MATCHEDGAUSS:
createMatchedGaussDistribution(numberOfLocalParticles, massIneV); createMatchedGaussDistribution(numberOfLocalParticles, massIneV, charge);
break; break;
case DistrTypeT::FROMFILE: case DistrTypeT::FROMFILE:
createDistributionFromFile(numberOfParticles, massIneV); createDistributionFromFile(numberOfParticles, massIneV);
...@@ -1152,7 +1152,10 @@ void Distribution::createDistributionFromFile(size_t /*numberOfParticles*/, doub ...@@ -1152,7 +1152,10 @@ void Distribution::createDistributionFromFile(size_t /*numberOfParticles*/, doub
} }
void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, double massIneV) { void Distribution::createMatchedGaussDistribution(size_t numberOfParticles,
double massIneV,
double charge)
{
/* /*
ToDo: ToDo:
...@@ -1250,7 +1253,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub ...@@ -1250,7 +1253,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
typedef boost::numeric::odeint::runge_kutta4<container_t> rk4_t; typedef boost::numeric::odeint::runge_kutta4<container_t> rk4_t;
typedef ClosedOrbitFinder<double,unsigned int, rk4_t> cof_t; typedef ClosedOrbitFinder<double,unsigned int, rk4_t> cof_t;
cof_t cof(massIneV*1E-6, Nint, CyclotronElement, full, Nsectors); cof_t cof(massIneV*1E-6, charge, Nint, CyclotronElement, full, Nsectors);
cof.findOrbit(accuracy, maxitCOF, E_m*1E-6, denergy, rguess, true); cof.findOrbit(accuracy, maxitCOF, E_m*1E-6, denergy, rguess, true);
throw EarlyLeaveException("Distribution::CreateMatchedGaussDistribution()", throw EarlyLeaveException("Distribution::CreateMatchedGaussDistribution()",
...@@ -1266,6 +1269,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub ...@@ -1266,6 +1269,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
Attributes::getReal(itsAttr[Attrib::Distribution::ET])*1E6, Attributes::getReal(itsAttr[Attrib::Distribution::ET])*1E6,
E_m*1E-6, E_m*1E-6,
massIneV*1E-6, massIneV*1E-6,
charge,
CyclotronElement, CyclotronElement,
Nint, Nint,
Nsectors, Nsectors,
...@@ -1373,7 +1377,7 @@ void Distribution::createOpalCycl(PartBunchBase<double, 3> *beam, ...@@ -1373,7 +1377,7 @@ void Distribution::createOpalCycl(PartBunchBase<double, 3> *beam,
emitting_m = false; emitting_m = false;
// Create distribution. // Create distribution.
create(numberOfPartToCreate, beam->getM()); create(numberOfPartToCreate, beam->getM(), beam->getQ());
// this variable is needed to be compatible with OPAL-T // this variable is needed to be compatible with OPAL-T
particlesPerDist_m.push_back(tOrZDist_m.size()); particlesPerDist_m.push_back(tOrZDist_m.size());
...@@ -1453,11 +1457,11 @@ void Distribution::createOpalT(PartBunchBase<double, 3> *beam, ...@@ -1453,11 +1457,11 @@ void Distribution::createOpalT(PartBunchBase<double, 3> *beam,
setupEmissionModel(beam); setupEmissionModel(beam);
// Create distributions. // Create distributions.
create(particlesPerDist_m.at(0), beam->getM()); create(particlesPerDist_m.at(0), beam->getM(), beam->getQ());
size_t distCount = 1; size_t distCount = 1;
for (Distribution* addedDist : addedDistributions_m) { for (Distribution* addedDist : addedDistributions_m) {
addedDist->create(particlesPerDist_m.at(distCount), beam->getM()); addedDist->create(particlesPerDist_m.at(distCount), beam->getM(), beam->getQ());
distCount++; distCount++;
} }
......
...@@ -286,7 +286,13 @@ private: ...@@ -286,7 +286,13 @@ private:
void applyEmissModelAstra(double &px, double &py, double &pz, std::vector<double> &additionalRNs); void applyEmissModelAstra(double &px, double &py, double &pz, std::vector<double> &additionalRNs);
void applyEmissModelNone(double &pz); void applyEmissModelNone(double &pz);
void applyEmissModelNonEquil(double eZ, double &px, double &py, double &pz, std::vector<double> &additionalRNs); void applyEmissModelNonEquil(double eZ, double &px, double &py, double &pz, std::vector<double> &additionalRNs);
void create(size_t &numberOfParticles, double massIneV); /*!
* Create the particle distribution.
* @param numberOfParticles to create
* @param massIneV particle charge in eV
* @param charge of the particle type in elementary charge
*/
void create(size_t &numberOfParticles, double massIneV, double charge);
void calcPartPerDist(size_t numberOfParticles); void calcPartPerDist(size_t numberOfParticles);
void checkEmissionParameters(); void checkEmissionParameters();
void checkIfEmitted(); void checkIfEmitted();
...@@ -326,7 +332,8 @@ private: ...@@ -326,7 +332,8 @@ private:
void createDistributionMultiGauss(size_t numberOfParticles, double massIneV); void createDistributionMultiGauss(size_t numberOfParticles, double massIneV);
void createDistributionFromFile(size_t numberOfParticles, double massIneV); void createDistributionFromFile(size_t numberOfParticles, double massIneV);
void createDistributionGauss(size_t numberOfParticles, double massIneV); void createDistributionGauss(size_t numberOfParticles, double massIneV);
void createMatchedGaussDistribution(size_t numberOfParticles, double massIneV); void createMatchedGaussDistribution(size_t numberOfParticles,
double massIneV, double charge);
void sampleUniformDisk(gsl_qrng* quasiRandGen2D, double& x1, double& x2); void sampleUniformDisk(gsl_qrng* quasiRandGen2D, double& x1, double& x2);
void fillEBinHistogram(); void fillEBinHistogram();
void fillParticleBins(); void fillParticleBins();
......
...@@ -69,6 +69,7 @@ SigmaGenerator::SigmaGenerator(double I, ...@@ -69,6 +69,7 @@ SigmaGenerator::SigmaGenerator(double I,
double ez, double ez,
double E, double E,
double m, double m,
double q,
const Cyclotron* cycl, const Cyclotron* cycl,
unsigned int N, unsigned int N,
unsigned int Nsectors, unsigned int Nsectors,
...@@ -82,6 +83,7 @@ SigmaGenerator::SigmaGenerator(double I, ...@@ -82,6 +83,7 @@ SigmaGenerator::SigmaGenerator(double I,
, nh_m(cycl->getCyclHarm()) , nh_m(cycl->getCyclHarm())
, beta_m(std::sqrt(1.0-1.0/gamma2_m)) , beta_m(std::sqrt(1.0-1.0/gamma2_m))
, m_m(m) , m_m(m)
, q_m(q)
, niterations_m(0) , niterations_m(0)
, converged_m(false) , converged_m(false)
, Emin_m(cycl->getFMLowE()) , Emin_m(cycl->getFMLowE())
...@@ -131,7 +133,7 @@ SigmaGenerator::SigmaGenerator(double I, ...@@ -131,7 +133,7 @@ SigmaGenerator::SigmaGenerator(double I,
// formula (57) // formula (57)
double lam = Physics::two_pi*Physics::c / (wo_m * nh_m); // wavelength, [lam] = m double lam = Physics::two_pi*Physics::c / (wo_m * nh_m); // wavelength, [lam] = m
// [K3] = m // [K3] = m
double K3 = 3.0 * I_m * lam double K3 = 3.0 * std::abs(q_m) * I_m * lam
/ (20.0 * std::sqrt(5.0) * Physics::pi * Physics::epsilon_0 * m / (20.0 * std::sqrt(5.0) * Physics::pi * Physics::epsilon_0 * m
* Physics::c * beta_m * beta_m * gamma_m * gamma2_m); * Physics::c * beta_m * beta_m * gamma_m * gamma2_m);
...@@ -193,7 +195,7 @@ bool SigmaGenerator::match(double accuracy, ...@@ -193,7 +195,7 @@ bool SigmaGenerator::match(double accuracy,
std::vector<matrix_t> Mcycs(nSteps_m), Mscs(nSteps_m); std::vector<matrix_t> Mcycs(nSteps_m), Mscs(nSteps_m);
ClosedOrbitFinder<double, unsigned int, ClosedOrbitFinder<double, unsigned int,
boost::numeric::odeint::runge_kutta4<container_t> > cof(m_m, N_m, cycl, false, nSectors_m); boost::numeric::odeint::runge_kutta4<container_t> > cof(m_m, q_m, N_m, cycl, false, nSectors_m);
if ( !cof.findOrbit(accuracy, maxitOrbit, E_m, denergy, rguess) ) { if ( !cof.findOrbit(accuracy, maxitOrbit, E_m, denergy, rguess) ) {
throw OpalException("SigmaGenerator::match()", throw OpalException("SigmaGenerator::match()",
...@@ -519,13 +521,13 @@ void SigmaGenerator::initialize(double nuz, double ravg) ...@@ -519,13 +521,13 @@ void SigmaGenerator::initialize(double nuz, double ravg)
double lam = Physics::two_pi * Physics::c / (wo_m * nh_m); // wavelength, [lam] = m double lam = Physics::two_pi * Physics::c / (wo_m * nh_m); // wavelength, [lam] = m
// m * c^3 --> c^2 in [m] = eV / c^2 cancel --> m * c in denominator // m * c^3 --> c^2 in [m] = eV / c^2 cancel --> m * c in denominator
double K3 = 3.0 * I_m * lam double K3 = 3.0 * std::abs(q_m) * I_m * lam
/ (20.0 * std::sqrt(5.0) * Physics::pi * Physics::epsilon_0 * m / (20.0 * std::sqrt(5.0) * Physics::pi * Physics::epsilon_0 * m
* Physics::c * beta_m * beta_m * gamma2_m * gamma_m); // [K3] = m * Physics::c * beta_m * beta_m * gamma2_m * gamma_m); // [K3] = m
// c in denominator cancels with unit of [m] = eV / c^2 --> we need to multiply // c in denominator cancels with unit of [m] = eV / c^2 --> we need to multiply
// with c in order to get dimensionless quantity // with c in order to get dimensionless quantity
double alpha = Physics::mu_0 * I_m * Physics::c / (5.0 * std::sqrt(10.0) * m double alpha = std::abs(q_m) * Physics::mu_0 * I_m * Physics::c / (5.0 * std::sqrt(10.0) * m
* gamma_m * nh_m) * std::sqrt(rcyc * rcyc * rcyc / (e * e * e)); // [alpha] = 1/rad --> [alpha] = 1 * gamma_m * nh_m) * std::sqrt(rcyc * rcyc * rcyc / (e * e * e)); // [alpha] = 1/rad --> [alpha] = 1
double sig0 = std::sqrt(2.0 * rcyc * e) / gamma_m; // [sig0] = m*sqrt(rad) --> [sig0] = m double sig0 = std::sqrt(2.0 * rcyc * e) / gamma_m; // [sig0] = m*sqrt(rad) --> [sig0] = m
......
...@@ -71,6 +71,7 @@ public: ...@@ -71,6 +71,7 @@ public:
* @param ez is the emittance in z-direction (vertical), \f$ \left[\varepsilon_{z}\right] = \pi\ mm\ mrad \f$ * @param ez is the emittance in z-direction (vertical), \f$ \left[\varepsilon_{z}\right] = \pi\ mm\ mrad \f$
* @param E is the energy, \f$ \left[E\right] = MeV \f$ * @param E is the energy, \f$ \left[E\right] = MeV \f$
* @param m is the mass of the particles \f$ \left[m\right] = \frac{MeV}{c^{2}} \f$ * @param m is the mass of the particles \f$ \left[m\right] = \frac{MeV}{c^{2}} \f$
* @param q is the particle charge [e]
* @param cycl is the cyclotron element * @param cycl is the cyclotron element
* @param N is the number of integration steps (closed orbit computation). That's why its also the number * @param N is the number of integration steps (closed orbit computation). That's why its also the number
* of maps (for each integration step a map) * of maps (for each integration step a map)
...@@ -79,7 +80,7 @@ public: ...@@ -79,7 +80,7 @@ public:
* @param write is a boolean (default: true). If true all maps of all iterations are stored, otherwise not. * @param write is a boolean (default: true). If true all maps of all iterations are stored, otherwise not.
*/ */
SigmaGenerator(double I, double ex, double ey, double ez, SigmaGenerator(double I, double ex, double ey, double ez,
double E, double m, const Cyclotron* cycl, double E, double m, double q, const Cyclotron* cycl,
unsigned int N, unsigned int Nsectors, unsigned int truncOrder, bool write = true); unsigned int N, unsigned int Nsectors, unsigned int truncOrder, bool write = true);
/// Searches for a matched distribution. /// Searches for a matched distribution.
...@@ -160,6 +161,8 @@ private: ...@@ -160,6 +161,8 @@ private:
double beta_m; double beta_m;
/// Is the mass of the particles, \f$ \left[m\right] = \frac{MeV}{c^{2}} \f$ /// Is the mass of the particles, \f$ \left[m\right] = \frac{MeV}{c^{2}} \f$
double m_m; double m_m;
/// Is the particle charge [e]
double q_m;
/// Is the number of iterations needed for convergence /// Is the number of iterations needed for convergence
unsigned int niterations_m; unsigned int niterations_m;
/// Is true if converged, false otherwise /// Is true if converged, false otherwise
......
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