Commit 4e4a8a69 authored by frey_m's avatar frey_m

modified: ClosedOrbitFinder.h

new file:   CycMagneticField.h
modified:   MagneticField.h
modified:   MapGenerator.h
modified:   SigmaGenerator.h
parent 81e56d96
......@@ -31,6 +31,8 @@
#include "MagneticField.h" // ONLY FOR STAND-ALONE PROGRAM
#include "CycMagneticField.h"
#include <fstream>
......@@ -141,9 +143,9 @@ class ClosedOrbitFinder
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
std::array<value_type,2> px_m; // px_m = [px1, px2]
/// Stores current position in longitudinal direction for the solutions of the ODE with different initial values
/// Stores current position in vertical direction for the solutions of the ODE with different initial values
std::array<value_type,2> z_m; // z_m = [z1, z2]
/// Stores current momenta in longitudinal direction for the solutions of the ODE with different initial values
/// Stores current momenta in vertical direction for the solutions of the ODE with different initial values
std::array<value_type,2> pz_m; // pz_m = [pz1, pz2]
/// Stores the inverse bending radius
......@@ -798,7 +800,13 @@ void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeVerticalOscillati
invptheta = 1.0 / ptheta;
// intepolate values of magnetic field
std::cout << "vertical: " << y[4] << std::endl; std::cin.get();
bField_m.interpolate(bint, brint, btint, y[0], theta * 180.0 / Physics::pi);
std::cout << bint << " " << brint << " " << btint << std::endl;
bint *= invbcon;
brint *= invbcon;
btint *= invbcon;
......
#ifndef CYCMAGNETICFIELD_H
#define CYCMAGNETICFIELD_H
#include <cmath>
#include "AbsBeamline/Cyclotron.h"
class CycMagneticField : public Cyclotron {
public:
void interpolate(double& bint,
double& brint,
double& btint,
const double& r,
const double& theta);
private:
bool interpolate_m(const double& r,
const double& theta,
const double& z,
double& bint
);
};
void CycMagneticField::interpolate(double& bint,
double& brint,
double& btint,
const double& r,
const double& theta,
const double& z
)
{
interpolate_m(r, theta, z, bint);
// derivative w.r.t. the radius
double dr = 0.0001
double b1, b2;
double tmp = r + dr;
interpolate_m(tmp, theta, z, b1);
tmp = r - dr;
interpolate_m(tmp, theta, z, b2);
brint = (b2 - b1) / (2.0 * dr);
tmp = theta + dr;
interpolate_m(r, tmp, z, b1);
tmp = theta - dr;
interpolate_m(r, tmp, z, b2);
btint = (b2 - b1) / (2.0 * dr);
}
bool CycMagneticField::interpolate_m(const double& rad,
const double& theta,
const double& z,
double& bint
)
{
// x horizontal
// y longitudinal
// z is vertical
const double xir = (rad - BP.rmin) / BP.delr;
// ir : the mumber of path whoes radius is less then the 4 points of cell which surrond particle.
const int ir = (int)xir;
// wr1 : the relative distance to the inner path radius
const double wr1 = xir - (double)ir;
// wr2 : the relative distance to the outer path radius
const double wr2 = 1.0 - wr1;
// double tet_rad = theta;
// the actual angle of particle
theta = theta / 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.
tet_map = fmod(theta, 360.0 / symmetry_m);
xit = tet_map / BP.dtet;
int it = (int) xit;
// *gmsg << R << " tet_map= " << tet_map << " ir= " << ir << " it= " << it << " bf= " << Bfield.bfld[idx(ir,it)] << endl;
const double wt1 = xit - (double)it;
const double wt2 = 1.0 - wt1;
// it : the number of point on the inner path whoes angle is less then the particle' corresponding angle.
// include zero degree point
it = it + 1;
int r1t1, r2t1, r1t2, r2t2;
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) {
/*
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 = idx(ir, it);
r2t1 = idx(ir + 1, it);
r1t2 = idx(ir, it + 1);
r2t2 = idx(ir + 1, it + 1);
}
double brf = 0.0
bint = 0.0;
if((it >= 0) && (ir >= 0) && (it < Bfield.ntetS) && (ir < Bfield.nrad)) {
/* Br */
brf = (Bfield.dbr[r1t1] * wr2 * wt2 +
Bfield.dbr[r2t1] * wr1 * wt2 +
Bfield.dbr[r1t2] * wr2 * wt1 +
Bfield.dbr[r2t2] * wr1 * wt1) * z;
bint = - brf;
} else {
return true;
}
return false;
}
#endif
\ No newline at end of file
......@@ -228,7 +228,7 @@ void MagneticField<Value_type>::read(const Value_type& scaleFactor) {
}
fp1.close();
const gsl_interp2d_type *T = gsl_interp2d_bilinear;
const gsl_interp2d_type *T = gsl_interp2d_bicubic; //gsl_interp2d_bilinear;
spline = gsl_spline2d_alloc(T, Bfield.nrad, Bfield.ntet);
xacc = gsl_interp_accel_alloc();
yacc = gsl_interp_accel_alloc();
......
......@@ -5,7 +5,7 @@
* It has one template parameter that specifies the type of the variables and containers.
*
* @author Matthias Frey
* @version 1.0
* @version 1.1
*/
#ifndef MAPGENERATOR_H
#define MAPGENERATOR_H
......@@ -21,6 +21,9 @@
#include "FixedAlgebra/FTps.h"
#include "FixedAlgebra/FMatrix.h"
#include "Physics/Physics.h"
#include "Utilities/OpalException.h"
/// @brief This class generates the matrices for the one turn matrix of a cyclotron.
template<typename Value_type, typename Size_type, typename Series_type,
typename Map_type, typename Hamiltonian_type, typename Space_charge_type>
......@@ -45,7 +48,7 @@ class MapGenerator
/// Type for specifying vectors
typedef std::vector<value_type> vector_type;
/// Constructor
/// Initialize
/*!
* @param nMaps is the number of maps
*/
......@@ -59,16 +62,31 @@ class MapGenerator
*/
matrix_type generateMap(const series_type&, value_type, size_type);
/// Combines the space charge maps (for each angle one) and the cyclotron maps (for each angle one) to the ont turn map, taking lists of maps
/// Combine given maps
/*!
* Combines the space charge maps (for each angle one) and the cyclotron maps (for each angle one)
* to the ont turn map, taking lists of maps
* @param Mscs is a list of space charge maps (the higher the index, the higher the angle)
* @param Mcycs is a list of cyclotron maps (the higher the index, the higher the angle)
*/
void combine(std::vector<matrix_type>&, std::vector<matrix_type>&);
/*!
* Combine given container of maps.
* @param maps to be combined.
*/
matrix_type combine(std::vector<matrix_type>& maps);
/// Returns the one turn map
matrix_type getMap();
/*!
* Compute the radial and vertical tune from the map.
* @param map from where to compute the tunes.
* @returns the radial and vertical tunes (in this order)
*/
std::pair<value_type, value_type> computeTunes(const matrix_type& map);
private:
/// Number of maps
size_type nMaps_m;
......@@ -83,13 +101,41 @@ class MapGenerator
// -----------------------------------------------------------------------------------------------------------------------
template<typename Value_type, typename Size_type, typename Series_type, typename Map_type, typename Hamiltonian_type, typename Space_charge_type>
MapGenerator<Value_type, Size_type, Series_type, Map_type, Hamiltonian_type, Space_charge_type>::MapGenerator(size_type nMaps)
template<typename Value_type,
typename Size_type,
typename Series_type,
typename Map_type,
typename Hamiltonian_type,
typename Space_charge_type>
MapGenerator<Value_type,
Size_type,
Series_type,
Map_type,
Hamiltonian_type,
Space_charge_type>::MapGenerator(size_type nMaps)
: nMaps_m(nMaps)
{}
template<typename Value_type, typename Size_type, typename Series_type, typename Map_type, typename Hamiltonian_type, typename Space_charge_type>
typename MapGenerator<Value_type, Size_type, Series_type, Map_type, Hamiltonian_type, Space_charge_type>::matrix_type MapGenerator<Value_type, Size_type, Series_type, Map_type, Hamiltonian_type, Space_charge_type>::generateMap/*cyc*/(const series_type& H, value_type ds, size_type order) {
template<typename Value_type,
typename Size_type,
typename Series_type,
typename Map_type,
typename Hamiltonian_type,
typename Space_charge_type>
typename MapGenerator<Value_type,
Size_type,
Series_type,
Map_type,
Hamiltonian_type,
Space_charge_type>::matrix_type
MapGenerator<Value_type,
Size_type,
Series_type,
Map_type,
Hamiltonian_type,
Space_charge_type>::generateMap/*cyc*/(const series_type& H,
value_type ds,
size_type order) {
// expand
map_type M = ExpMap(-H * ds, order);
......@@ -108,11 +154,22 @@ typename MapGenerator<Value_type, Size_type, Series_type, Map_type, Hamiltonian_
return map;
}
template<typename Value_type, typename Size_type, typename Series_type, typename Map_type, typename Hamiltonian_type, typename Space_charge_type>
void MapGenerator<Value_type, Size_type, Series_type, Map_type, Hamiltonian_type, Space_charge_type>::combine(std::vector<matrix_type>& Mscs, std::vector<matrix_type>& Mcycs) {
template<typename Value_type,
typename Size_type,
typename Series_type,
typename Map_type,
typename Hamiltonian_type,
typename Space_charge_type>
void MapGenerator<Value_type,
Size_type,
Series_type,
Map_type,
Hamiltonian_type,
Space_charge_type>::combine(std::vector<matrix_type>& Mscs,
std::vector<matrix_type>& Mcycs) {
if (nMaps_m != Mscs.size() || nMaps_m != Mcycs.size())
throw std::length_error("Error in MapGenerator::combine: Wrong vector dimensions.");
throw OpalException("MapGenerator::combine()", "Wrong vector dimensions.");
Mturn_m = boost::numeric::ublas::identity_matrix<value_type>(6);
......@@ -120,9 +177,107 @@ void MapGenerator<Value_type, Size_type, Series_type, Map_type, Hamiltonian_type
Mturn_m = matt_boost::gemmm<matrix_type>(Mscs[i],Mcycs[i],Mturn_m);
}
template<typename Value_type, typename Size_type, typename Series_type, typename Map_type, typename Hamiltonian_type, typename Space_charge_type>
typename MapGenerator<Value_type, Size_type, Series_type, Map_type, Hamiltonian_type, Space_charge_type>::matrix_type MapGenerator<Value_type, Size_type, Series_type, Map_type, Hamiltonian_type, Space_charge_type>::getMap() {
template<typename Value_type,
typename Size_type,
typename Series_type,
typename Map_type,
typename Hamiltonian_type,
typename Space_charge_type>
typename MapGenerator<Value_type,
Size_type,
Series_type,
Map_type,
Hamiltonian_type,
Space_charge_type>::matrix_type
MapGenerator<Value_type,
Size_type,
Series_type,
Map_type,
Hamiltonian_type,
Space_charge_type>::combine(std::vector<matrix_type>& maps)
{
matrix_type map = boost::numeric::ublas::identity_matrix<value_type>(6);
for (std::size_t i = 0; i < maps.size(); ++i) {
matrix_type tmp = prod(maps[i], map);
map = tmp;
}
return map;
}
template<typename Value_type,
typename Size_type,
typename Series_type,
typename Map_type,
typename Hamiltonian_type,
typename Space_charge_type>
typename MapGenerator<Value_type,
Size_type,
Series_type,
Map_type,
Hamiltonian_type,
Space_charge_type>::matrix_type
MapGenerator<Value_type,
Size_type,
Series_type,
Map_type,
Hamiltonian_type,
Space_charge_type>::getMap() {
return Mturn_m;
}
template<typename Value_type,
typename Size_type,
typename Series_type,
typename Map_type,
typename Hamiltonian_type,
typename Space_charge_type>
std::pair<Value_type,
Value_type> MapGenerator<Value_type,
Size_type,
Series_type,
Map_type,
Hamiltonian_type,
Space_charge_type>::computeTunes(const matrix_type& map)
{
/*
* M = [ cos(mu) + alpha * sin(mu) beta * sin(mu)
* -gamma * sin(mu) cos(mu) - alpha * sin(mu)]
*
* i = 0, 2, 4
* --> cos(mu) = 0.5 * [ M(i, i) + M(i + 1, i + 1) ]
*/
value_type arg = 0.0;
// horizontal phase advance [rad]
arg = 0.5 * ( map(0, 0) + map(1, 1) );
if ( std::abs( arg ) > 1 )
throw OpalException("MapGenerator::computeTunes()",
"Horizontal phase advance: Acos argument " + std::to_string(arg) + " out of range.");
value_type mux = std::acos( arg );
// vertical phase advance [rad]
arg = 0.5 * ( map(2, 2) + map(3, 3) );
if ( std::abs( arg ) > 1 )
throw OpalException("MapGenerator::computeTunes()",
"Vertical phase advance: Acos argument " + std::to_string(arg) + " out of range.");
value_type muz = std::acos( arg );
return std::make_pair<value_type,
value_type>( mux * Physics::u_two_pi,
muz * Physics::u_two_pi );
}
#endif
This diff is collapsed.
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