Commit 81e56d96 authored by frey_m's avatar frey_m

Adapted code to new MagneticField class.

modified:   ClosedOrbitFinder.h
modified:   Distribution.cpp
deleted:    MagneticField.cpp
modified:   MagneticField.h
modified:   SigmaGenerator.h
parent be287ffc
......@@ -64,12 +64,13 @@ class ClosedOrbitFinder
* @param Emax is the maximum energy [MeV] reached in cyclotron
* @param nSector is the number of sectors (--> symmetry) of cyclotron
* @param fmapfn is the location of the file that specifies the magnetic field
* @param guesss value of radius for closed orbit finder
* @param guesss value of radius for closed orbit finder
* @param scaleFactor for the magnetic field (default: 1.0)
* @param domain is a boolean (default: true). If "true" the closed orbit is computed over a single sector,
* otherwise over 2*pi.
*/
ClosedOrbitFinder(value_type, value_type, value_type, size_type, value_type, size_type, value_type, value_type, size_type,
const std::string&, value_type, bool = true);
const std::string&, value_type, value_type scaleFactor = 1.0, bool = true);
/// Returns the inverse bending radius (size of container N+1)
container_type& getInverseBendingRadius();
......@@ -251,17 +252,19 @@ class ClosedOrbitFinder
// PUBLIC MEMBER FUNCTIONS
// -----------------------------------------------------------------------------------------------------------------------
template<typename Value_type, typename Size_type, class Stepper>
ClosedOrbitFinder<Value_type, Size_type, Stepper>::ClosedOrbitFinder(value_type E, value_type E0, value_type wo, size_type N,
value_type accuracy, size_type maxit,
value_type Emin, value_type Emax, size_type nSector,
const std::string& fmapfn,
value_type rguess,
bool domain)
template<typename Value_type, typename Size_type, class Stepper>
ClosedOrbitFinder<Value_type,
Size_type,
Stepper>::ClosedOrbitFinder(value_type E, value_type E0,
value_type wo, size_type N,
value_type accuracy, size_type maxit,
value_type Emin, value_type Emax,
size_type nSector, const std::string& fmapfn,
value_type rguess, value_type scaleFactor, bool domain)
: nxcross_m(0), nzcross_m(0), E_m(E), E0_m(E0), wo_m(wo), N_m(N), dtheta_m(Physics::two_pi/value_type(N)),
gamma_m(E/E0+1.0), ravg_m(0), phase_m(0), converged_m(false), Emin_m(Emin), Emax_m(Emax), nSector_m(nSector),
lastOrbitVal_m(0.0), lastMomentumVal_m(0.0),
vertOscDone_m(false), domain_m(domain), stepper_m(), rguess_m(rguess)
vertOscDone_m(false), domain_m(domain), stepper_m(), rguess_m(rguess), bField_m(fmapfn)
{
if ( Emin_m > Emax_m )
......@@ -298,7 +301,7 @@ ClosedOrbitFinder<Value_type, Size_type, Stepper>::ClosedOrbitFinder(value_type
fidx_m.reserve(N_m);
// read in magnetic fieldmap
bField_m.read(fmapfn);
bField_m.read(scaleFactor);
// compute closed orbit
converged_m = findOrbit(accuracy, maxit);
......
......@@ -48,7 +48,6 @@
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include "MagneticField.h"
extern Inform *gmsg;
......@@ -1333,11 +1332,6 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
hE = E_m*1E-6;
}
int nr, nth, nsc;
double rmin, dr, dth;
MagneticField::ReadHeader(&nr, &nth, &rmin, &dr, &dth, &nsc, Attributes::getString(itsAttr[AttributesT::FMAPFN]));
int Nint = 1000;
bool writeMap = true;
......@@ -1352,11 +1346,7 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
lE,
hE,
(int)Attributes::getReal(itsAttr[AttributesT::MAGSYM]),
rmin,
Nint,
nth,
nr,
dr,
Attributes::getString(itsAttr[AttributesT::FMAPFN]),
Attributes::getReal(itsAttr[AttributesT::ORDERMAPS]),
writeMap);
......
#include "MagneticField.h"
template<typename Value_type>
void MagneticField<Value_type>::read(const double &scaleFactor) {
FILE *f = NULL;
*gmsg << "* ----------------------------------------------" << endl;
*gmsg << "* READ IN CARBON CYCLOTRON FIELD MAP " << endl;
*gmsg << "* ----------------------------------------------" << endl;
BP.Bfact = scaleFactor;
if((f = fopen(fmapfn_m.c_str(), "r")) == NULL) {
throw GeneralClassicException("Cyclotron::getFieldFromFile_Carbon",
"failed to open file '" + fmapfn_m + "', please check if it exists");
}
CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.rmin));
*gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.delr));
//if the value is negative, the actual value is its reciprocal.
if(BP.delr < 0.0) BP.delr = 1.0 / (-BP.delr);
*gmsg << "* Stepsize in radial direction: " << BP.delr << " [mm]" << endl;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.tetmin));
*gmsg << "* Minimal angle of measured field map: " << BP.tetmin << " [deg]" << endl;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.dtet));
//if the value is negative, the actual value is its reciprocal.
if(BP.dtet < 0.0) BP.dtet = 1.0 / (-BP.dtet);
*gmsg << "* Stepsize in azimuthal direction: " << BP.dtet << " [deg]" << endl;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%d", &Bfield.ntet));
*gmsg << "* Grid points along azimuth (ntet): " << Bfield.ntet << endl;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%d", &Bfield.nrad));
*gmsg << "* Grid points along radius (nrad): " << Bfield.nrad << endl;
//Bfield.ntetS = Bfield.ntet;
Bfield.ntetS = Bfield.ntet + 1;
*gmsg << "* Accordingly, total grid point along azimuth: " << Bfield.ntetS << endl;
//Bfield.ntot = idx(Bfield.nrad - 1, Bfield.ntet) + 1;
Bfield.ntot = Bfield.nrad * Bfield.ntetS;
*gmsg << "* Adding a guard cell along azimuth" << endl;
*gmsg << "* Total stored grid point number ((ntet+1) * nrad) : " << Bfield.ntot << endl;
Bfield.bfld.resize(Bfield.ntot);
Bfield.dbt.resize(Bfield.ntot);
Bfield.dbtt.resize(Bfield.ntot);
Bfield.dbttt.resize(Bfield.ntot);
*gmsg << "* rescaling of the fields with factor: " << BP.Bfact << endl;
std::vector<Value_type> r(Bfield.nrad), theta(Bfield.ntet), bval(Bfield.ntot);
for(int i = 0; i < Bfield.nrad; i++) {
for(int k = 0; k < Bfield.ntet; k++) {
CHECK_CYC_FSCANF_EOF(fscanf(f, "%16lE", &(Bfield.bfld[idx(i, k)])));
Bfield.bfld[idx(i, k)] *= BP.Bfact;
}
}
std::fstream fp1;
fp1.open("gnu.out", std::ios::out);
for (int i = 0; i < Bfield.nrad; ++i) {
for (int k = 0; k < Bfield.ntet; ++k) {
fp1 << BP.rmin + (i * BP.delr) << " \t " << k * (BP.tetmin + BP.dtet) << " \t " << Bfield.bfld[idx(i, k)] << endl;
}
}
fp1.close();
const gsl_interp2d_type *T = gsl_interp2d_bilinear;
spline = gsl_spline2d_alloc(T, Bfield.nrad, Bfield.ntet);
xacc = gsl_interp_accel_alloc();
yacc = gsl_interp_accel_alloc();
for (int i = 0; i < Bfield.ntet; ++i)
theta[i] = i * (BP.tetmin + BP.dtet);
int cnt = 0;
for (int i = 0; i < Bfield.nrad; ++i) {
r[i] = BP.rmin * 0.001 + (i * BP.delr * 0.001);
for (int k = 0; k < Bfield.ntet; ++k) {
bval[Bfield.nrad * k + i] = Bfield.bfld[idx(i, k)];
}
}
gsl_spline2d_init(spline, &r[0], &theta[0], &bval[0], r.size(), theta.size());
fclose(f);
*gmsg << "* Field Maps read successfully!" << endl << endl;
}
template<typename Value_type>
void MagneticField<Value_type>::interpolate(Value_type& bint,
Value_type& brint,
Value_type& btint,
Value_type r,
Value_type theta) {
bint = 0;
brint = 0;
btint = 0;
// return if we are out of the radius range
if ( r < BP.rmin * 0.001 || r > (BP.rmin + (Bfield.nrad - 1)* BP.delr) * 0.001 )
return;
// scale the angle to the correct range
if ( theta < 0 )
theta = 360.0 + theta;
if ( theta > BP.tetmin + (Bfield.ntet - 1) * BP.dtet ) {
Value_type max = BP.tetmin + (Bfield.ntet - 1) * BP.dtet;
theta -= int(theta / max) * max;
}
bint = gsl_spline2d_eval(spline, r, theta, xacc, yacc);
// derivative w.r.t. the radius, i.e. dB/dr
brint = gsl_spline2d_eval_deriv_x(spline, r, theta, xacc, yacc);
// derivative w.r.t. the azimuth angle, i.e. dB/dtheta
btint = gsl_spline2d_eval_deriv_y(spline, r, theta, xacc, yacc);
return 0;
}
\ No newline at end of file
......@@ -25,19 +25,32 @@
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
// #include "Utilities/LogicalError.h"
#include "Utilities/GeneralClassicException.h"
// #include <Ippl.h>
/// Memory allocator
#define MALLOC(n,a) ((a*)malloc((n)*sizeof(a)))
extern Inform *gmsg;
#define CHECK_CYC_FSCANF_EOF(arg) if(arg == EOF)\
throw GeneralClassicException("Cyclotron::getFieldFromFile",\
"fscanf returned EOF at " #arg);
struct BfieldData {
template<typename Value_type>
class MagneticField
{
public:
struct BPositions {
// this 4 parameters are need to be read from field file.
double rmin, delr;
double tetmin, dtet;
// Radii and step width of initial Grid
std::vector<double> rarr;
// int ThetaPeriodicity; // Periodicity of Magnetic field
double Bfact; // MULTIPLICATION FACTOR FOR MAGNETIC FIELD
};
struct BfieldData {
std::string filename;
// known from file: field and three theta derivatives
//~ double *bfld; //Bz
......@@ -91,22 +104,8 @@ struct BfieldData {
// Mean and Maximas
double bacc, dbtmx, dbttmx, dbtttmx;
};
template<typename Value_type>
class MagneticField
{
public:
struct BPositions {
// this 4 parameters are need to be read from field file.
double rmin, delr;
double tetmin, dtet;
// Radii and step width of initial Grid
std::vector<double> rarr;
// int ThetaPeriodicity; // Periodicity of Magnetic field
double Bfact; // MULTIPLICATION FACTOR FOR MAGNETIC FIELD
};
MagneticField(const std::string fmapfn);
/*!
* Read in a CARBONCYCL fieldmap.
......@@ -121,7 +120,7 @@ public:
* @param theta is the azimuth angle where to interpolate
*/
void interpolate(Value_type& bint, Value_type& brint, Value_type& btint,
const Value_type& r, const Value_type& theta);
const Value_type& r, Value_type theta);
~MagneticField() {
......@@ -129,6 +128,10 @@ public:
gsl_interp_accel_free(xacc);
gsl_interp_accel_free(yacc);
}
inline int getNradSteps();
inline int getNtetSteps();
private:
......@@ -149,6 +152,145 @@ private:
};
template<typename Value_type>
MagneticField<Value_type>::MagneticField(const std::string fmapfn) :
fmapfn_m(fmapfn)
{ }
template<typename Value_type>
void MagneticField<Value_type>::read(const Value_type& scaleFactor) {
FILE *f = NULL;
*gmsg << "* ----------------------------------------------" << endl;
*gmsg << "* READ IN CARBON CYCLOTRON FIELD MAP " << endl;
*gmsg << "* ----------------------------------------------" << endl;
BP.Bfact = scaleFactor;
if((f = fopen(fmapfn_m.c_str(), "r")) == NULL) {
throw GeneralClassicException("Cyclotron::getFieldFromFile_Carbon",
"failed to open file '" + fmapfn_m + "', please check if it exists");
}
CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.rmin));
*gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.delr));
//if the value is negative, the actual value is its reciprocal.
if(BP.delr < 0.0) BP.delr = 1.0 / (-BP.delr);
*gmsg << "* Stepsize in radial direction: " << BP.delr << " [mm]" << endl;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.tetmin));
*gmsg << "* Minimal angle of measured field map: " << BP.tetmin << " [deg]" << endl;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.dtet));
//if the value is negative, the actual value is its reciprocal.
if(BP.dtet < 0.0) BP.dtet = 1.0 / (-BP.dtet);
*gmsg << "* Stepsize in azimuthal direction: " << BP.dtet << " [deg]" << endl;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%d", &Bfield.ntet));
*gmsg << "* Grid points along azimuth (ntet): " << Bfield.ntet << endl;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%d", &Bfield.nrad));
*gmsg << "* Grid points along radius (nrad): " << Bfield.nrad << endl;
//Bfield.ntetS = Bfield.ntet;
Bfield.ntetS = Bfield.ntet + 1;
*gmsg << "* Accordingly, total grid point along azimuth: " << Bfield.ntetS << endl;
//Bfield.ntot = idx(Bfield.nrad - 1, Bfield.ntet) + 1;
Bfield.ntot = Bfield.nrad * Bfield.ntetS;
*gmsg << "* Adding a guard cell along azimuth" << endl;
*gmsg << "* Total stored grid point number ((ntet+1) * nrad) : " << Bfield.ntot << endl;
Bfield.bfld.resize(Bfield.ntot);
Bfield.dbt.resize(Bfield.ntot);
Bfield.dbtt.resize(Bfield.ntot);
Bfield.dbttt.resize(Bfield.ntot);
*gmsg << "* rescaling of the fields with factor: " << BP.Bfact << endl;
std::vector<Value_type> r(Bfield.nrad), theta(Bfield.ntet), bval(Bfield.ntot);
for(int i = 0; i < Bfield.nrad; i++) {
for(int k = 0; k < Bfield.ntet; k++) {
CHECK_CYC_FSCANF_EOF(fscanf(f, "%16lE", &(Bfield.bfld[idx(i, k)])));
Bfield.bfld[idx(i, k)] *= BP.Bfact;
}
}
std::fstream fp1;
fp1.open("gnu.out", std::ios::out);
for (int i = 0; i < Bfield.nrad; ++i) {
for (int k = 0; k < Bfield.ntet; ++k) {
fp1 << BP.rmin + (i * BP.delr) << " \t " << k * (BP.tetmin + BP.dtet)
<< " \t " << Bfield.bfld[idx(i, k)] << std::endl;
}
}
fp1.close();
const gsl_interp2d_type *T = gsl_interp2d_bilinear;
spline = gsl_spline2d_alloc(T, Bfield.nrad, Bfield.ntet);
xacc = gsl_interp_accel_alloc();
yacc = gsl_interp_accel_alloc();
for (int i = 0; i < Bfield.ntet; ++i)
theta[i] = i * (BP.tetmin + BP.dtet);
for (int i = 0; i < Bfield.nrad; ++i) {
r[i] = BP.rmin * 0.001 + (i * BP.delr * 0.001);
for (int k = 0; k < Bfield.ntet; ++k) {
bval[Bfield.nrad * k + i] = Bfield.bfld[idx(i, k)];
}
}
gsl_spline2d_init(spline, &r[0], &theta[0], &bval[0], r.size(), theta.size());
fclose(f);
*gmsg << "* Field Maps read successfully!" << endl << endl;
}
template<typename Value_type>
void MagneticField<Value_type>::interpolate(Value_type& bint,
Value_type& brint,
Value_type& btint,
const Value_type &r,
Value_type theta) {
bint = 0;
brint = 0;
btint = 0;
// return if we are out of the radius range
if ( r < BP.rmin * 0.001 || r > (BP.rmin + (Bfield.nrad - 1)* BP.delr) * 0.001 )
return;
// scale the angle to the correct range
if ( theta < 0 )
theta = 360.0 + theta;
if ( theta > BP.tetmin + (Bfield.ntet - 1) * BP.dtet ) {
Value_type max = BP.tetmin + (Bfield.ntet - 1) * BP.dtet;
theta -= int(theta / max) * max;
}
bint = gsl_spline2d_eval(spline, r, theta, xacc, yacc);
// derivative w.r.t. the radius, i.e. dB/dr
brint = gsl_spline2d_eval_deriv_x(spline, r, theta, xacc, yacc);
// derivative w.r.t. the azimuth angle, i.e. dB/dtheta
btint = gsl_spline2d_eval_deriv_y(spline, r, theta, xacc, yacc);
}
template<typename Value_type>
int MagneticField<Value_type>::getNradSteps() {
return Bfield.nrad;
}
template<typename Value_type>
int MagneticField<Value_type>::getNtetSteps() {
return Bfield.ntet;
}
#endif
......@@ -95,21 +95,25 @@ class SigmaGenerator
* @param Emin is the minimum energy [MeV] needed in cyclotron, \f$ \left[E_{min}\right] = MeV \f$
* @param Emax is the maximum energy [MeV] reached in cyclotron, \f$ \left[E_{max}\right] = MeV \f$
* @param nSector is the number of sectors (symmetry assumption)
* @param rmin is the minimal radius of the cyclotron, \f$ \left[r_{min}\right] = m \f$
* @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)
* @param ntheta is the number of angle splits (fieldmap variable)
* @param nradial is the number of radial splits (fieldmap variable)
* @param dr is the radial step size, \f$ \left[\Delta r}\right] = m \f$
* @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)
* @param fieldmap is the location of the file that specifies the magnetic field
* @param truncOrder is the truncation order for power series of the Hamiltonian
* @param scaleFactor for the magnetic field (default: 1.0)
* @param write is a boolean (default: true). If true all maps of all iterations are stored, otherwise not.
*/
SigmaGenerator(value_type, value_type, value_type, value_type, value_type, value_type, value_type, value_type, value_type, value_type, size_type, value_type, size_type, size_type, size_type, value_type, const std::string&, size_type, bool = true);
SigmaGenerator(value_type I, value_type ex, value_type ey, value_type ez,
value_type wo, value_type E, value_type nh, value_type m,
value_type Emin, value_type Emax, size_type nSector,
size_type N, const std::string& fieldmap,
size_type truncOrder, value_type scaleFactor = 1.0, bool write = true);
/// Searches for a matched distribution. Returns "true" if a matched distribution within the accuracy could be found, returns "false" otherwise
/// Searches for a matched distribution.
/*!
* Returns "true" if a matched distribution within the accuracy could be found, returns "false" otherwise.
* @param accuracy is used for computing the equilibrium orbit (to a certain accuracy)
* @param maxit is the maximum number of iterations (the program stops if within this value no stationary distribution was found)
* @param maxit is the maximum number of iterations (the program stops if within this value no stationary
* distribution was found)
* @param maxitOrbit is the maximum number of iterations for finding closed orbit
* @param angle defines the start of the sector (one can choose any angle between 0° and 360°)
* @param guess value of radius for closed orbit finder
......@@ -118,7 +122,8 @@ class SigmaGenerator
bool match(value_type, size_type, size_type, value_type, value_type, bool);
/// Block diagonalizes the symplex part of the one turn transfer matrix
/** It computes the transformation matrix <b>R</b> and its inverse <b>invR</b>. It returns a vector containing the four eigenvalues
/** It computes the transformation matrix <b>R</b> and its inverse <b>invR</b>.
* It returns a vector containing the four eigenvalues
* (alpha,beta,gamma,delta) (see formula (38), paper: Geometrical method of decoupling)
*/
/*!
......@@ -128,8 +133,9 @@ class SigmaGenerator
*/
vector_type decouple(const matrix_type&, sparse_matrix_type&, sparse_matrix_type&);
/// Checks if the sigma-matrix is an eigenellipse and returns the L2 error (The idea of this function is taken from Dr. Christian Baumgarten's program).
/// Checks if the sigma-matrix is an eigenellipse and returns the L2 error.
/*!
* The idea of this function is taken from Dr. Christian Baumgarten's program.
* @param Mturn is the one turn transfer matrix
* @param sigma is the sigma matrix to be tested
*/
......@@ -176,18 +182,10 @@ class SigmaGenerator
value_type Emax_m;
/// Number of (symmetric) sectors
size_type nSector_m;
/// Minimal radius of cyclotron, \f$ \left[r_{min}\right] = m \f$
value_type rmin_m;
/// Number of integration steps for closed orbit computation
size_type N_m;
/// Number of integration steps per sector (--> also: number of maps per sector)
size_type nStepsPerSector_m;
/// Number of angle splits (fieldmap)
size_type ntheta_m;
/// Number of radial splits (fieldmap)
size_type nradial_m;
/// Radial step size, \f$ \left[\Delta r\right] = m \f$
value_type dr_m;
/// Error of computation
value_type error_m;
......@@ -199,6 +197,9 @@ class SigmaGenerator
/// Decides for writing output (default: true)
bool write_m;
/// Scale the magnetic field values (default: 1.0)
value_type scaleFactor_m;
/// Stores the stationary distribution (sigma matrix)
matrix_type sigma_m;
......@@ -268,14 +269,15 @@ class SigmaGenerator
template<typename Value_type, typename Size_type>
SigmaGenerator<Value_type, Size_type>::SigmaGenerator(value_type I, value_type ex, value_type ey, value_type ez, value_type wo,
value_type E, value_type nh, value_type m, value_type Emin, value_type Emax, size_type nSector, value_type rmin,
size_type N, size_type ntheta, size_type nradial, value_type dr, const std::string& fieldmap, size_type truncOrder, bool write)
value_type E, value_type nh, value_type m, value_type Emin, value_type Emax, size_type nSector,
size_type N, const std::string& fieldmap, size_type truncOrder, value_type scaleFactor, bool write)
: I_m(I), wo_m(wo), E_m(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), rmin_m(rmin), N_m(N), nStepsPerSector_m(N/nSector), ntheta_m(ntheta),
nradial_m(nradial), dr_m(dr), error_m(std::numeric_limits<value_type>::max()),
fieldmap_m(fieldmap), truncOrder_m(truncOrder), write_m(write), sigmas_m(nStepsPerSector_m)
Emin_m(Emin), Emax_m(Emax), nSector_m(nSector), N_m(N), nStepsPerSector_m(N/nSector),
error_m(std::numeric_limits<value_type>::max()),
fieldmap_m(fieldmap), truncOrder_m(truncOrder), write_m(write),
scaleFactor_m(scaleFactor), sigmas_m(nStepsPerSector_m)
{
// set emittances (initialization like that due to old compiler version)
// [ex] = [ey] = [ez] = pi*mm*mrad --> [emittance] = mm mrad
......@@ -365,11 +367,11 @@ template<typename Value_type, typename Size_type>
std::pair<value_type,value_type> tunes;
if (!harmonic) {
ClosedOrbitFinder<value_type, size_type, boost::numeric::odeint::runge_kutta4<container_type> > cof(E_m, m_m, wo_m, N_m, accuracy,
maxitOrbit, Emin_m, Emax_m,
nSector_m, rmin_m, ntheta_m,
nradial_m, dr_m, fieldmap_m, rguess,
false);
ClosedOrbitFinder<value_type, size_type,
boost::numeric::odeint::runge_kutta4<container_type> > cof(E_m, m_m, wo_m, N_m, accuracy,
maxitOrbit, Emin_m, Emax_m,
nSector_m, fieldmap_m, rguess,
scaleFactor_m, false);
// properties of one turn
container_type h_turn = cof.getInverseBendingRadius();
......@@ -457,7 +459,13 @@ template<typename Value_type, typename Size_type>
std::copy_n(ds_turn.begin()+start,nStepsPerSector_m, ds.begin());
} else {
Harmonics<value_type, size_type> H(wo_m,Emin_m, Emax_m, ntheta_m, nradial_m, nSector_m, E_m, m_m);
MagneticField<double> bField(fieldmap_m);
bField.read(scaleFactor_m);
int nradial = bField.getNradSteps();
int ntheta = bField.getNtetSteps();
Harmonics<value_type, size_type> H(wo_m,Emin_m, Emax_m, ntheta, nradial, nSector_m, E_m, m_m);
Mcycs = H.computeMap("data/inj2sym_mainharms.4","data/inj2sym_mainharms.8",4);
ravg = H.getRadius();
tunes = H.getTunes();
......
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