Commit 73120351 authored by frey_m's avatar frey_m

Merge branch '593-closedorbitfinder-sigmagenerator-charge-not-applied-correctly' into 'master'

Resolve "ClosedOrbitFinder/SigmaGenerator: Charge not applied correctly"

See merge request !420
parents 37da44f0 773f6668
......@@ -70,6 +70,7 @@ class ClosedOrbitFinder
/// Sets the initial values for the integration and calls findOrbit().
/*!
* @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 cycl is the cyclotron element
* @param domain is a boolean (default: true). If "true" the closed orbit is computed over a single sector,
......@@ -77,7 +78,7 @@ class ClosedOrbitFinder
* @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,
ClosedOrbitFinder(value_type E0, value_type q, size_type N,
Cyclotron* cycl, bool domain = true, int Nsectors = 1);
/// Returns the inverse bending radius (size of container N+1)
......@@ -154,7 +155,7 @@ class ClosedOrbitFinder
// void computeVerticalOscillations();
/// 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
std::array<value_type,2> x_m; // x_m = [x1, x2]
......@@ -188,6 +189,9 @@ class ClosedOrbitFinder
/// Is the rest mass [MeV / c**2]
value_type E0_m;
// Is the particle charge [e]
value_type q_m;
/// Is the nominal orbital frequency
/* (see paper of Dr. C. Baumgarten: "Transverse-Longitudinal
* Coupling by Space Charge in Cyclotrons" (2012), formula (1))
......@@ -243,8 +247,8 @@ class ClosedOrbitFinder
/*!
* 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) {
return e0 * 1.0e7 / (/* physics::q0 */ 1.0 * Physics::c * Physics::c / wo);
std::function<double(double, double)> bcon_m = [this](double e0, double wo) {
return e0 * 1.0e7 / (q_m * Physics::c * Physics::c / wo);
};
Cyclotron* cycl_m;
......@@ -258,11 +262,13 @@ template<typename Value_type, typename Size_type, class Stepper>
ClosedOrbitFinder<Value_type,
Size_type,
Stepper>::ClosedOrbitFinder(value_type E0,
value_type q,
size_type N, Cyclotron* cycl,
bool domain, int nSectors)
: nxcross_m(0)
, nzcross_m(0)
, E0_m(E0)
, q_m(q)
, wo_m(cycl->getRfFrequ(0)*1E6/cycl->getCyclHarm()*2.0*Physics::pi)
, N_m(N)
, dtheta_m(Physics::two_pi/value_type(N))
......@@ -894,15 +900,21 @@ void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeOrbitProperties(c
template<typename Value_type, typename Size_type, class Stepper>
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;
// compute the number of steps per degree
value_type deg_step = N_m / 360.0;
if (angle < 0.0) {
angle = 360.0 + angle;
}
// compute starting point
size_type start = deg_step * angle;
unsigned int start = deg_step * angle;
start %= orbitProperty.size();
// copy end to start
std::copy(orbitProperty.begin() + start, orbitProperty.end(), orbitPropertyCopy.begin());
......
......@@ -248,7 +248,7 @@ void Distribution::execute() {
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
......@@ -284,7 +284,7 @@ void Distribution::create(size_t &numberOfParticles, double massIneV) {
switch (distrTypeT_m) {
case DistrTypeT::MATCHEDGAUSS:
createMatchedGaussDistribution(numberOfLocalParticles, massIneV);
createMatchedGaussDistribution(numberOfLocalParticles, massIneV, charge);
break;
case DistrTypeT::FROMFILE:
createDistributionFromFile(numberOfParticles, massIneV);
......@@ -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:
......@@ -1250,7 +1253,7 @@ 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, full, Nsectors);
cof_t cof(massIneV*1E-6, charge, Nint, CyclotronElement, full, Nsectors);
cof.findOrbit(accuracy, maxitCOF, E_m*1E-6, denergy, rguess, true);
throw EarlyLeaveException("Distribution::CreateMatchedGaussDistribution()",
......@@ -1266,6 +1269,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
Attributes::getReal(itsAttr[Attrib::Distribution::ET])*1E6,
E_m*1E-6,
massIneV*1E-6,
charge,
CyclotronElement,
Nint,
Nsectors,
......@@ -1373,7 +1377,7 @@ void Distribution::createOpalCycl(PartBunchBase<double, 3> *beam,
emitting_m = false;
// Create distribution.
create(numberOfPartToCreate, beam->getM());
create(numberOfPartToCreate, beam->getM(), beam->getQ());
// this variable is needed to be compatible with OPAL-T
particlesPerDist_m.push_back(tOrZDist_m.size());
......@@ -1453,11 +1457,11 @@ void Distribution::createOpalT(PartBunchBase<double, 3> *beam,
setupEmissionModel(beam);
// Create distributions.
create(particlesPerDist_m.at(0), beam->getM());
create(particlesPerDist_m.at(0), beam->getM(), beam->getQ());
size_t distCount = 1;
for (Distribution* addedDist : addedDistributions_m) {
addedDist->create(particlesPerDist_m.at(distCount), beam->getM());
addedDist->create(particlesPerDist_m.at(distCount), beam->getM(), beam->getQ());
distCount++;
}
......
......@@ -286,7 +286,13 @@ private:
void applyEmissModelAstra(double &px, double &py, double &pz, std::vector<double> &additionalRNs);
void applyEmissModelNone(double &pz);
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 checkEmissionParameters();
void checkIfEmitted();
......@@ -326,7 +332,8 @@ private:
void createDistributionMultiGauss(size_t numberOfParticles, double massIneV);
void createDistributionFromFile(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 fillEBinHistogram();
void fillParticleBins();
......
......@@ -69,6 +69,7 @@ SigmaGenerator::SigmaGenerator(double I,
double ez,
double E,
double m,
double q,
const Cyclotron* cycl,
unsigned int N,
unsigned int Nsectors,
......@@ -82,6 +83,7 @@ SigmaGenerator::SigmaGenerator(double I,
, nh_m(cycl->getCyclHarm())
, beta_m(std::sqrt(1.0-1.0/gamma2_m))
, m_m(m)
, q_m(q)
, niterations_m(0)
, converged_m(false)
, Emin_m(cycl->getFMLowE())
......@@ -131,7 +133,7 @@ SigmaGenerator::SigmaGenerator(double I,
// formula (57)
double lam = Physics::two_pi*Physics::c / (wo_m * nh_m); // wavelength, [lam] = 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
* Physics::c * beta_m * beta_m * gamma_m * gamma2_m);
......@@ -193,7 +195,7 @@ bool SigmaGenerator::match(double accuracy,
std::vector<matrix_t> Mcycs(nSteps_m), Mscs(nSteps_m);
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) ) {
throw OpalException("SigmaGenerator::match()",
......@@ -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
// 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
* 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
// 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
double sig0 = std::sqrt(2.0 * rcyc * e) / gamma_m; // [sig0] = m*sqrt(rad) --> [sig0] = m
......
......@@ -71,6 +71,7 @@ public:
* @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 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 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)
......@@ -79,7 +80,7 @@ public:
* @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,
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);
/// Searches for a matched distribution.
......@@ -160,6 +161,8 @@ private:
double beta_m;
/// Is the mass of the particles, \f$ \left[m\right] = \frac{MeV}{c^{2}} \f$
double m_m;
/// Is the particle charge [e]
double q_m;
/// Is the number of iterations needed for convergence
unsigned int niterations_m;
/// 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