Commit ffd03f26 authored by cortes_c's avatar cortes_c
Browse files

getting rid of the bug which caused a discrepancy between OPAL Single mode...

getting rid of the bug which caused a discrepancy between OPAL Single mode tracking and Closed Orbit Finder
parent 619e4aef
......@@ -8,6 +8,9 @@
*
* @author Matthias Frey
* @version 1.0
* @version 1.1
* Added rotate(), for the case where theta != 0
* Fixed bugg in getOrbit()
*/
#ifndef CLOSEDORBITFINDER_H
......@@ -74,13 +77,13 @@ class ClosedOrbitFinder
bool domain = true);
/// Returns the inverse bending radius (size of container N+1)
container_type& getInverseBendingRadius();
container_type& getInverseBendingRadius(value_type angle = 0);
/// Returns the step lengths of the path (size of container N+1)
container_type& getPathLength();
container_type& getPathLength(value_type angle = 0);
/// Returns the field index (size of container N+1)
container_type& getFieldIndex();
container_type& getFieldIndex(value_type angle = 0);
/// Returns the radial and vertical tunes (in that order)
std::pair<value_type,value_type> getTunes();
......@@ -138,6 +141,9 @@ class ClosedOrbitFinder
/// This function computes nzcross_ which is used to compute the tune in z-direction and the frequency error
void computeVerticalOscillations();
/// This function rotates the calculated closed orbit finder properties to the initial angle
container_type rotate(value_type angle, 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]
/// Stores current momenta in horizontal direction for the solutions of the ODE with different initial values
......@@ -268,7 +274,6 @@ ClosedOrbitFinder<Value_type,
lastOrbitVal_m(0.0), lastMomentumVal_m(0.0),
vertOscDone_m(false), domain_m(domain), stepper_m(), rguess_m(rguess), bField_m(fmapfn, nSector)
{
if ( Emin_m > Emax_m )
throw OpalException("ClosedOrbitFinder::ClosedOrbitFinder()",
"Incorrect cyclotron energy (MeV) bounds: Maximum cyclotron energy smaller than minimum cyclotron energy.");
......@@ -285,9 +290,7 @@ ClosedOrbitFinder<Value_type,
throw OpalException("ClosedOrbitFinder::ClosedOrbitFinder()", "Relativistic factor equal zero.");
// if domain_m = true --> integrate over a single sector
if (domain_m) {
N_m /= nSector_m;
}
if (domain_m) N_m /= nSector_m;
// reserve storage for the orbit and momentum (--> size = 0, capacity = N_m+1)
/*
......@@ -314,22 +317,25 @@ ClosedOrbitFinder<Value_type,
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>::getInverseBendingRadius()
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getInverseBendingRadius(value_type angle)
{
if (angle != 0.0) h_m = rotate(angle, h_m);
return h_m;
}
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>::getPathLength()
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getPathLength(value_type angle)
{
if (angle != 0.0) ds_m = rotate(angle, ds_m);
return ds_m;
}
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>::getFieldIndex()
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getFieldIndex(value_type angle)
{
if (angle != 0.0) fidx_m = rotate(angle, fidx_m);
return fidx_m;
}
......@@ -352,23 +358,8 @@ 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>::getOrbit(value_type angle)
{
container_type r = r_m;
if (angle != 0.0) {
// compute the number of steps per degree
value_type deg_step = N_m / 360.0;
// compute starting point
size_type start = deg_step * angle;
// copy end to start
std::copy(r_m.begin() + start, r_m.end(), r.begin());
// copy start to end
std::copy_n(r_m.begin(), start, r.end() - start);
}
return r;
if (angle != 0.0) r_m = rotate(angle, r_m);
return r_m;
}
template<typename Value_type, typename Size_type, class Stepper>
......@@ -377,18 +368,7 @@ inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_typ
{
container_type pr = pr_m;
if (angle != 0.0) {
// compute the number of steps per degree
value_type deg_step = N_m / 360.0;
// compute starting point
size_type start = deg_step * angle;
// copy end to start
std::copy(pr_m.begin() + start, pr_m.end(), pr.begin());
// copy start to end
std::copy_n(pr_m.begin(), start, pr.end() - start);
}
if (angle != 0.0) pr = rotate(angle, pr_m);
// change units from meters to \beta * \gamma
/* in Gordon paper:
......@@ -407,10 +387,9 @@ inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_typ
*
* The momentum in \beta * \gamma is obtained by dividing by "a"
*/
value_type factor = 1.0 / acon_m(wo_m);
std::for_each(pr.begin(), pr.end(), [factor](value_type p) { return p * factor; });
value_type factor = 1.0 / acon_m(wo_m);
std::for_each(pr.begin(), pr.end(), [factor](value_type& p) { p *= factor; });
return pr;
}
......@@ -433,7 +412,7 @@ template<typename Value_type, typename Size_type, class Stepper>
typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::value_type
ClosedOrbitFinder<Value_type, Size_type, Stepper>::getFrequencyError()
{
// if the vertical oscillations aren't computed, we have to, since there we also compuote the frequency error.
// if the vertical oscillations aren't computed, we have to, since there we also compute the frequency error.
if(!vertOscDone_m)
computeVerticalOscillations();
......@@ -495,7 +474,6 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
const double theta)
{
pr2 = y[1] * y[1];
if (p2 < pr2)
throw OpalException("ClosedOrbitFinder::findOrbit()", "p_{r}^2 > p^{2} (defined in Gordon paper) --> Square root of negative number.");
......@@ -504,7 +482,7 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
invptheta = 1.0 / ptheta;
// interpolate values of magnetic field
bField_m.interpolate(bint, brint, btint, y[0], theta * 180.0 / Physics::pi);
bField_m.interpolate(bint, brint, btint, y[0], theta /* * 180.0 / Physics::pi*/);
bint *= invbcon;
brint *= invbcon;
......@@ -514,7 +492,7 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
// Gordon, formula (5b)
dydt[1] = ptheta - y[0] * bint;
// Gordon, formulas (9a) and (9b)
for (size_type i = 2; i < 5; i += 2) {
for (size_type i = 2; i < 6; i += 2) {
dydt[i] = (y[1] * y[i] + y[0] * p2 * y[i+1] * invptheta * invptheta) * invptheta;
dydt[i+1] = - y[1] * y[i+1] * invptheta - (bint + y[0] * brint) * y[i];
}
......@@ -522,7 +500,6 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
// define initial state container for integration: y = {r, pr, x1, px1, x2, px2}
state_type y(6);
// difference of last and first value of r (1. element) and pr (2. element)
container_type err(2);
// correction term for initial values: r = r + dr, pr = pr + dpr; Gordon, formula (17)
......@@ -595,13 +572,13 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
px_m[1] = 1.0; // px2; Gordon, formula (10)
nxcross_m = 0; // counts the number of crossings of x-axis (excluding first step)
idx = 0; // index for looping over r and pr arrays
//real_theta = 0.0;
// fill container with initial states
y = {init[0],init[1], x_m[0], px_m[0], x_m[1], px_m[1]};
// integrate from 0 to 2*pi (one has to get back to the "origin")
boost::numeric::odeint::integrate_n_steps(stepper_m,orbit_integration,y,0.0,dtheta_m,N_m,store);
// write new state
x_m[0] = y[2];
px_m[0] = y[3];
......@@ -624,6 +601,7 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
// compute amplitude of the error
error = std::sqrt(delta[0] * delta[0] + delta[1] * delta[1] * invgamma4) / r_m[0];
} while (error > accuracy && niterations++ < maxit);
// reset iteration counter
......@@ -662,7 +640,27 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
// returns true if converged, otherwise false
return error < accuracy;
}
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){
container_type orbitPropertyCopy = orbitProperty;
// compute the number of steps per degree
value_type deg_step = N_m / 360.0;
// compute starting point
size_type start = deg_step * angle;
// copy end to start
std::copy(orbitProperty.begin() + start, orbitProperty.end(), orbitPropertyCopy.begin());
// copy start to end
std::copy_n(orbitProperty.begin(), start, orbitPropertyCopy.end() - start);
return orbitPropertyCopy;
}
template<typename Value_type, typename Size_type, class Stepper>
Value_type ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeTune(const std::array<value_type,2>& y,
value_type py2, size_type ncross)
......@@ -746,7 +744,8 @@ void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeOrbitProperties()
for (size_type i = 0; i < N_m; ++i) {
// interpolate magnetic field
bField_m.interpolate(bint, brint, btint, r_m[i], theta * 180.0 / Physics::pi);
bField_m.interpolate(bint, brint, btint, r_m[i], theta);
bint *= invbcon;
brint *= invbcon;
btint *= invbcon;
......@@ -764,7 +763,6 @@ void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeOrbitProperties()
// increase angle
theta += dtheta_m;
}
// compute average radius
ravg_m = std::accumulate(r_m.begin(),r_m.end(),0.0) / value_type(r_m.size());
}
......@@ -804,7 +802,7 @@ void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeVerticalOscillati
invptheta = 1.0 / ptheta;
// intepolate values of magnetic field
bField_m.interpolate(bint, brint, btint, y[0], theta * 180.0 / Physics::pi);
bField_m.interpolate(bint, brint, btint, y[0], theta);
bint *= invbcon;
brint *= invbcon;
......
......@@ -1301,7 +1301,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
<< " PHIINIT= " << CyclotronElement->getPHIinit() << endl;
*gmsg << "* ----------------------------------------------------" << endl;
const double wo = CyclotronElement->getRfFrequ()*1E6/CyclotronElement->getCyclHarm()*2.0*Physics::pi;
const double wo = CyclotronElement->getRfFrequ()*1E6*2.0*Physics::pi/CyclotronElement->getCyclHarm();
const double fmLowE = CyclotronElement->getFMLowE();
const double fmHighE = CyclotronElement->getFMHighE();
......@@ -1312,7 +1312,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
"'CYCLOTRON' definition.");
}
int Nint = 1440;
int Nint = 1440*2;
double scaleFactor = 1.0;
bool writeMap = true;
......
......@@ -50,7 +50,7 @@ void MagneticField::interpolate(double& bint,
// x horizontal
// y longitudinal
// z is vertical
const double xir = (rad - BP.rmin) / (BP.delr);
const double xir = (rad - BP.rmin) / BP.delr;
// ir : the number of path whose radius is less than the 4 points of cell which surround the particle.
const int ir = (int)xir;
......@@ -60,15 +60,13 @@ void MagneticField::interpolate(double& bint,
// wr2 : the relative distance to the outer path radius
const double wr2 = 1.0 - wr1;
double tet_rad = theta;
// the actual angle of particle
// tet_rad = theta / Physics::pi * 180.0;
double tet = theta / Physics::pi * 180.0;
// the corresponding angle on the field map
// Note: this does not work if the start point of field map does not equal zero.
double tet_map = fmod(tet_rad, 360.0 / this->getSymmetry());
double tet_map = fmod(tet, 360.0 / this->getSymmetry());
double xit = tet_map / BP.dtet;
......@@ -82,30 +80,30 @@ void MagneticField::interpolate(double& bint,
it = it + 1;
int r1t1, r2t1, r1t2, r2t2;
// int ntetS = Bfield.ntet + 1;
int ntetS = Bfield.ntet + 1;
// r1t1 : the index of the "min angle, min radius" point in the 2D field array.
// considering the array start with index of zero, minus 1.
// if(myBFieldType_m != FFAGBF) {
if(1) {
/*
For FFAG this does not work
*/
// r1t1 = it + ntetS * ir - 1;
// r1t2 = r1t1 + 1;
// r2t1 = r1t1 + ntetS;
// r2t2 = r2t1 + 1 ;
// } else {
// /*
// With this we have B-field AND this is far more
// intuitive for me ....
// */
r1t1 = it + ntetS * ir - 1;
r1t2 = r1t1 + 1;
r2t1 = r1t1 + ntetS;
r2t2 = r2t1 + 1 ;
} else {
/*
With this we have B-field AND this is far more
intuitive for me ....
*/
r1t1 = idx(ir, it);
r2t1 = idx(ir + 1, it);
r1t2 = idx(ir, it + 1);
r2t2 = idx(ir + 1, it + 1);
// }
}
bint = 0.0;
brint = 0.0;
......@@ -120,16 +118,17 @@ void MagneticField::interpolate(double& bint,
Bfield.dbr[r2t2] * wr1 * wt1);
// dB_{z}/dtheta
btint = Bfield.dbt[r1t1] * wr2 * wt2 +
btint = (Bfield.dbt[r1t1] * wr2 * wt2 +
Bfield.dbt[r2t1] * wr1 * wt2 +
Bfield.dbt[r1t2] * wr2 * wt1 +
Bfield.dbt[r2t2] * wr1 * wt1;
Bfield.dbt[r2t2] * wr1 * wt1);
// B_{z}
bint = Bfield.bfld[r1t1] * wr2 * wt2 +
bint = (Bfield.bfld[r1t1] * wr2 * wt2 +
Bfield.bfld[r2t1] * wr1 * wt2 +
Bfield.bfld[r1t2] * wr2 * wt1 +
Bfield.bfld[r2t2] * wr1 * wt1;
Bfield.bfld[r2t2] * wr1 * wt1);
//bint *= -1.0;
}
}
......
......@@ -57,6 +57,7 @@
#include <gsl/gsl_eigen.h>
#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
#if BOOST_VERSION >= 106000
#include <boost/numeric/odeint/integrate/check_adapter.hpp>
#endif
......@@ -392,7 +393,7 @@ SigmaGenerator<Value_type, Size_type>::SigmaGenerator(value_type I, value_type e
gamma_m(E/m+1.0), gamma2_m(gamma_m*gamma_m),
nh_m(nh), beta_m(std::sqrt(1.0-1.0/gamma2_m)), m_m(m), niterations_m(0), converged_m(false),
Emin_m(Emin), Emax_m(Emax), nSector_m(nSector), N_m(N), nStepsPerSector_m(N/nSector),
error_m(std::numeric_limits<value_type>::max()),
error_m(1),
fieldmap_m(fieldmap), truncOrder_m(truncOrder), write_m(write),
scaleFactor_m(scaleFactor), sigmas_m(0),permutations_m{{0,1,2,3},{1,2,3,0},{0,1,3,2},{1,0,2,3},{0,2,3,1},{0,3,2,1},
{1,0,3,2},{1,3,2,0},{2,1,0,3},{2,0,1,3},{2,3,0,1},{2,3,1,0},{3,1,0,2},{3,2,0,1}, {3,2,1,0},{3,0,1,2}}
......@@ -502,16 +503,16 @@ 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, wo_m, N_m, accuracy,
boost::numeric::odeint::runge_kutta_cash_karp54 <container_type> > cof(E_m, m_m, wo_m, N_m, accuracy,
maxitOrbit, Emin_m, Emax_m,
nSector_m, fieldmap_m, rguess,
type, scaleFactor_m, false);
// properties of one turn
container_type h_turn = cof.getInverseBendingRadius();
container_type h_turn = cof.getInverseBendingRadius(angle);
container_type r_turn = cof.getOrbit(angle);
container_type ds_turn = cof.getPathLength();
container_type fidx_turn = cof.getFieldIndex();
container_type ds_turn = cof.getPathLength(angle);
container_type fidx_turn = cof.getFieldIndex(angle);
tunes = cof.getTunes();
ravg = cof.getAverageRadius(); // average radius
......@@ -529,21 +530,12 @@ template<typename Value_type, typename Size_type>
<< "* vertical tune: " << tunes.second << endl
<< "* ----------------------------" << endl << endl;
// compute the number of steps per degree
value_type deg_step = N_m / 360.0;
// compute starting point of computation
size_type start = deg_step * angle;
// copy properties of the length of one sector (--> nSteps_m)
size_type step = 0;
for (size_type i = 0; i < nSteps_m; ++i) {
step = (start + i) % r_turn.size();
r[i] = r_turn[step];
h[i] = h_turn[step];
fidx[i] = fidx_turn[step];
ds[i] = ds_turn[step];
r[i] = r_turn[i];
h[i] = h_turn[i];
fidx[i] = fidx_turn[i];
ds[i] = ds_turn[i];
}
if(write_m){
......@@ -1116,6 +1108,10 @@ template<typename Value_type, typename Size_type> void SigmaGenerator<Value_typ
}
oneTurn<<std::endl;
oneTurn.close();
std::ofstream error("data/error"+energy+"MeV.dat",std::ios::app);
error << error_m << std::endl;
error.close();
}
template<typename Value_type, typename Size_type> void SigmaGenerator<Value_type, Size_type>::writeMatched(const matrix_type& match, const size_type permutation,const std::vector<matrix_type>& Mcycs, const std::vector<matrix_type>& Mscs, const matrix_type& Mcyc){
......
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