14 #ifndef SIGMAGENERATOR_H
15 #define SIGMAGENERATOR_H
36 #include <boost/numeric/ublas/matrix.hpp>
37 #include <boost/numeric/ublas/matrix_sparse.hpp>
38 #include <boost/numeric/ublas/vector.hpp>
40 #include <boost/numeric/ublas/vector_proxy.hpp>
41 #include <boost/numeric/ublas/triangular.hpp>
42 #include <boost/numeric/ublas/lu.hpp>
43 #include <boost/numeric/ublas/io.hpp>
45 #include <gsl/gsl_matrix.h>
46 #include <gsl/gsl_math.h>
47 #include <gsl/gsl_eigen.h>
49 #include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
50 #if BOOST_VERSION >= 106000
51 #include <boost/numeric/odeint/integrate/check_adapter.hpp>
63 template<
typename Value_type,
typename Size_type>
74 typedef boost::numeric::ublas::matrix<value_type>
matrix_type;
82 typedef boost::numeric::ublas::compressed_matrix<
complex_t,
83 boost::numeric::ublas::row_major
86 typedef boost::numeric::ublas::vector<value_type>
vector_type;
94 typedef std::function<Series(value_type,value_type,value_type)>
Hamiltonian;
96 typedef std::function<Series(value_type,value_type,value_type)>
SpaceCharge;
275 const std::vector<matrix_type>&);
322 template<
typename Value_type,
typename Size_type>
335 , wo_m(cycl->getRfFrequ()*1E6/cycl->getCyclHarm()*2.0*Physics::
pi)
338 , gamma2_m(gamma_m*gamma_m)
339 , nh_m(cycl->getCyclHarm())
340 , beta_m(std::
sqrt(1.0-1.0/gamma2_m))
344 , Emin_m(cycl->getFMLowE())
345 , Emax_m(cycl->getFMHighE())
347 , nSectors_m(Nsectors)
348 , nStepsPerSector_m(N/cycl->getSymmetry())
351 , truncOrder_m(truncOrder)
353 , sigmas_m(nStepsPerSector_m)
417 return -0.5 * Kx *
x_m *
x_m
423 template<
typename Value_type,
typename Size_type>
430 bool harmonic,
bool full)
441 nSteps_m = nStepsPerSector_m;
452 std::vector<matrix_type> Mcycs(nSteps_m), Mscs(nSteps_m);
454 container_type h(nSteps_m), r(nSteps_m), ds(nSteps_m), fidx(nSteps_m);
455 value_type ravg = 0.0, const_ds = 0.0;
456 std::pair<value_type,value_type> tunes;
460 boost::numeric::odeint::runge_kutta4<container_type> > cof(m_m, N_m, cycl,
false, nSectors_m);
462 if ( !cof.findOrbit(accuracy, maxitOrbit, E_m, denergy, rguess) ) {
464 "Closed orbit finder didn't converge.");
467 cof.computeOrbitProperties(E_m);
470 tunes = cof.getTunes();
471 ravg = cof.getAverageRadius();
483 writeOrbitOutput_m(tunes, ravg, cof.getFrequencyError(),
484 r_turn, peo, h_turn, fidx_turn, ds_turn);
487 *gmsg <<
"* ----------------------------" <<
endl
488 <<
"* Closed orbit info:" <<
endl
490 <<
"* average radius: " << ravg <<
" [m]" <<
endl
491 <<
"* initial radius: " << r_turn[0] <<
" [m]" <<
endl
492 <<
"* initial momentum: " << peo[0] <<
" [Beta Gamma]" <<
endl
493 <<
"* frequency error: " << cof.getFrequencyError()*100 <<
" [ % ] "<<
endl
494 <<
"* horizontal tune: " << tunes.first <<
endl
495 <<
"* vertical tune: " << tunes.second <<
endl
496 <<
"* ----------------------------" <<
endl <<
endl;
499 std::copy_n(r_turn.begin(), nSteps_m, r.begin());
500 std::copy_n(h_turn.begin(), nSteps_m, h.begin());
501 std::copy_n(fidx_turn.begin(), nSteps_m, fidx.begin());
502 std::copy_n(ds_turn.begin(), nSteps_m, ds.begin());
508 *gmsg <<
"Not yet supported." <<
endl;
513 initialize(tunes.second,ravg);
516 std::ofstream writeMturn, writeMcyc, writeMsc;
520 std::string energy = float2string(E_m);
522 writeMturn.open(
"data/OneTurnMapForEnergy"+energy+
"MeV.dat",std::ios::app);
523 writeMsc.open(
"data/SpaceChargeMapPerAngleForEnergy"+energy+
"MeV.dat",std::ios::app);
524 writeMcyc.open(
"data/CyclotronMapPerAngleForEnergy"+energy+
"MeV.dat",std::ios::app);
526 writeMturn <<
"--------------------------------" <<
std::endl;
527 writeMturn <<
"Iteration: 0 " <<
std::endl;
528 writeMturn <<
"--------------------------------" <<
std::endl;
530 writeMsc <<
"--------------------------------" <<
std::endl;
531 writeMsc <<
"Iteration: 0 " <<
std::endl;
532 writeMsc <<
"--------------------------------" <<
std::endl;
536 for (size_type i = 0; i < nSteps_m; ++i) {
538 Mcycs[i] = mapgen.generateMap(H_m(h[i],
543 Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),
548 Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),
551 const_ds,truncOrder_m);
561 mapgen.combine(Mscs,Mcycs);
576 value_type weight = 0.5;
578 while (error_m > accuracy && !stop) {
580 decouple(Mturn,
R,invR);
583 newSigma = updateInitialSigma(Mturn,
R, invR);
586 updateSigma(Mscs,Mcycs);
589 error_m = L2ErrorNorm(sigmas_m[0],newSigma);
592 sigmas_m[0] = weight*newSigma + (1.0-weight)*sigmas_m[0];
595 writeMsc <<
"--------------------------------" <<
std::endl;
596 writeMsc <<
"Iteration: " << niterations_m + 1 <<
std::endl;
597 writeMsc <<
"--------------------------------" <<
std::endl;
601 for (size_type i = 0; i < nSteps_m; ++i) {
603 Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),
608 Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),
611 const_ds,truncOrder_m);
619 mapgen.combine(Mscs,Mcycs);
620 Mturn = mapgen.getMap();
623 writeMturn <<
"--------------------------------" <<
std::endl;
624 writeMturn <<
"Iteration: " << niterations_m + 1 <<
std::endl;
625 writeMturn <<
"--------------------------------" <<
std::endl;
630 stop = (niterations_m++ > maxit);
634 sigma_m.resize(6,6,
false);
635 sigma_m.swap(newSigma);
638 converged_m = error_m < accuracy;
647 }
catch(
const std::exception&
e) {
651 if ( converged_m && write_m ) {
653 std::ofstream writeSigmaMatched(
"data/MatchedDistributions.dat", std::ios::app);
655 std::array<double,3> emit = this->getEmittances();
657 writeSigmaMatched <<
"* Converged (Ex, Ey, Ez) = (" << emit[0]
658 <<
", " << emit[1] <<
", " << emit[2]
659 <<
") pi mm mrad for E= " << E_m <<
" (MeV)"
662 for(
unsigned int i = 0; i < sigma_m.size1(); ++ i) {
663 writeSigmaMatched << std::setprecision(4) << std::setw(11)
665 for(
unsigned int j = 1; j < sigma_m.size2(); ++ j) {
666 writeSigmaMatched <<
" & " << std::setprecision(4)
667 << std::setw(11) << sigma_m(i,j);
669 writeSigmaMatched <<
" \\\\" <<
std::endl;
673 writeSigmaMatched.close();
680 template<
typename Value_type,
typename Size_type>
684 typedef gsl_matrix* gsl_matrix_t;
685 typedef gsl_vector_complex* gsl_cvector_t;
686 typedef gsl_matrix_complex* gsl_cmatrix_t;
687 typedef gsl_eigen_nonsymmv_workspace* gsl_wspace_t;
688 typedef boost::numeric::ublas::vector<complex_t> complex_vector_type;
690 gsl_cvector_t evals = gsl_vector_complex_alloc(6);
691 gsl_cmatrix_t evecs = gsl_matrix_complex_alloc(6, 6);
692 gsl_wspace_t wspace = gsl_eigen_nonsymmv_alloc(6);
693 gsl_matrix_t M = gsl_matrix_alloc(6, 6);
698 gsl_matrix_set(M, i, j, Mturn(i,j));
702 gsl_eigen_nonsymmv(M, evals, evecs, wspace);
716 gsl_vector_complex_view evec_i = gsl_matrix_complex_column(evecs, i);
719 gsl_complex zgsl = gsl_vector_complex_get(&evec_i.vector, j);
720 complex_t z(GSL_REAL(zgsl), GSL_IMAG(zgsl));
728 bool isZdirection =
false;
729 std::vector<complex_vector_type> zVectors{};
730 std::vector<complex_vector_type> xyVectors{};
735 if(z == 0.) isZdirection =
true;
736 complex_vector_type v(6);
742 zVectors.push_back(v);
749 xyVectors.push_back(v);
751 isZdirection =
false;
755 if(zVectors.size() != 2)
757 "Couldn't find the vertical Eigenvectors.");
762 for(
size_type j = 0; j < 6; j++) norm += std::norm(xyVectors[i](j));
767 for(
size_type j = 0; j < 6; j++) norm += std::norm(zVectors[i](j));
772 R(i,0) = xyVectors[0](i);
773 R(i,1) = xyVectors[1](i);
774 R(i,2) = zVectors[0](i);
775 R(i,3) = zVectors[1](i);
776 R(i,4) = xyVectors[2](i);
777 R(i,5) = xyVectors[3](i);
780 gsl_vector_complex_free(evals);
781 gsl_matrix_complex_free(evecs);
782 gsl_eigen_nonsymmv_free(wspace);
787 template<
typename Value_type,
typename Size_type>
791 typedef boost::numeric::ublas::permutation_matrix<size_type> pmatrix_t;
797 pmatrix_t pm(A.size1());
800 int res = lu_factorize(A,pm);
806 invR.assign(boost::numeric::ublas::identity_matrix<complex_t>(A.size1()));
809 boost::numeric::ublas::lu_substitute(A, pm, invR);
815 template<
typename Value_type,
typename Size_type>
820 this->eigsolve_m(Mturn, R);
822 if ( !this->invertMatrix_m(R, invR) )
824 "Couldn't compute inverse matrix!");
827 template<
typename Value_type,
typename Size_type>
833 matrix_type newSigma = matt_boost::gemmm<matrix_type>(Mturn,
835 boost::numeric::ublas::trans(Mturn));
838 return L2ErrorNorm(sigma,newSigma);
841 template<
typename Value_type,
typename Size_type>
848 template<
typename Value_type,
typename Size_type>
852 return (converged_m) ? niterations_m :
size_type(0);
855 template<
typename Value_type,
typename Size_type>
862 template<
typename Value_type,
typename Size_type>
863 inline std::array<Value_type,3>
867 return std::array<value_type,3>{{
878 template<
typename Value_type,
typename Size_type>
945 gamma_m * nh_m) *
std::sqrt(rcyc * rcyc * rcyc / (e * e * e));
952 sig = sig0 * std::cbrt(1.0 + alpha);
953 }
else if (alpha >= 0) {
954 sig = sig0 * (1 + alpha * (0.25 - 0.03125 *
alpha));
957 "Negative alpha value: " + std::to_string(alpha) +
" < 0");
961 value_type K = K3 * gamma_m / (3.0 * sig * sig * sig);
971 "Negative value --> No real-valued frequency.");
976 "Square root of negative number.");
982 "Square root of negative number.");
984 if (h * h * nuz * nuz <= K)
986 "h^{2} * nu_{z}^{2} <= K (i.e. square root of negative number)");
998 matrix_type sigma = boost::numeric::ublas::zero_matrix<value_type>(6);
1001 sigma(0,0) = invAB * (B * ex / Omega + A * ez / omega) * mega;
1004 sigma(0,5) = sigma(5,0) = invAB * (ex / Omega + ez / omega) * mega;
1007 sigma(1,1) = invAB * (B * ex * Omega + A * ez * omega) * mega;
1010 sigma(1,4) = sigma(4,1) = invAB * (ex * Omega+ez * omega) / (K * gamma2_m) * mega;
1013 sigma(2,2) = ey / (
std::sqrt(h * h * nuz * nuz - K)) * mega;
1015 sigma(3,3) = (ey * mega) * (ey * mega) / sigma(2,2);
1018 sigma(4,4) = invAB * (A * ex * Omega + B * ez * omega) / (K * gamma2_m) * mega;
1021 sigma(5,5) = invAB * (ex / (B * Omega) + ez / (A * omega)) * mega;
1024 sigmas_m.resize(nSteps_m);
1030 std::string energy = float2string(E_m);
1031 std::ofstream writeInit(
"data/maps/InitialSigmaPerAngleForEnergy" +
1032 energy +
"MeV.dat", std::ios::app);
1039 template<
typename Value_type,
typename Size_type>
1058 cmatrix_type S = boost::numeric::ublas::zero_matrix<complex_t>(6,6);
1059 S(0,1) = S(2,3) = S(4,5) = 1;
1060 S(1,0) = S(3,2) = S(5,4) = -1;
1067 cmatrix_type D = boost::numeric::ublas::zero_matrix<complex_t>(6,6);
1071 D(2*i, 2*i) = emittance_m[i] * invbg *
im;
1072 D(2*i+1, 2*i+1) = -emittance_m[i] * invbg *
im;
1085 sigma(i,j) = csigma(i,j).real();
1095 std::string energy = float2string(E_m);
1096 std::ofstream writeSigma(
"data/maps/SigmaPerAngleForEnergy" +
1097 energy +
"MeV.dat", std::ios::app);
1099 writeSigma <<
"--------------------------------" <<
std::endl;
1100 writeSigma <<
"Iteration: " << niterations_m + 1 <<
std::endl;
1101 writeSigma <<
"--------------------------------" <<
std::endl;
1110 template<
typename Value_type,
typename Size_type>
1112 const std::vector<matrix_type>& Mcycs)
1114 matrix_type M = boost::numeric::ublas::matrix<value_type>(6,6);
1116 std::ofstream writeSigma;
1119 std::string energy = float2string(E_m);
1120 writeSigma.open(
"data/maps/SigmaPerAngleForEnergy"+energy+
"MeV.dat",std::ios::app);
1124 for (
size_type i = 1; i < nSteps_m; ++i) {
1128 sigmas_m[i] = matt_boost::gemmm<matrix_type>(M,sigmas_m[i - 1],
1129 boost::numeric::ublas::trans(M));
1141 template<
typename Value_type,
typename Size_type>
1150 return std::sqrt(std::inner_product(diff.data().begin(),
1152 diff.data().begin(),
1156 template<
typename Value_type,
typename Size_type>
1170 return std::accumulate(diff.data().begin(), diff.data().end(), 0.0);
1174 template<
typename Value_type,
typename Size_type>
1176 std::ostringstream out;
1177 out << std::setprecision(6) << val;
1182 template<
typename Value_type,
typename Size_type>
1184 const std::pair<value_type,value_type>& tunes,
1194 std::ofstream writeTunes(
"data/Tunes.dat", std::ios::app);
1196 if(writeTunes.tellp() == 0)
1197 writeTunes <<
"energy [MeV]" << std::setw(15)
1198 <<
"nur" << std::setw(25)
1201 writeTunes << E_m << std::setw(30) << std::setprecision(10)
1202 << tunes.first << std::setw(25)
1206 std::ofstream writeAvgRadius(
"data/AverageValues.dat", std::ios::app);
1208 if (writeAvgRadius.tellp() == 0)
1209 writeAvgRadius <<
"energy [MeV]" << std::setw(15)
1210 <<
"avg. radius [m]" << std::setw(15)
1211 <<
"r [m]" << std::setw(15)
1214 writeAvgRadius << E_m << std::setw(25) << std::setprecision(10)
1215 << ravg << std::setw(25) << std::setprecision(10)
1216 << r_turn[0] << std::setw(25) << std::setprecision(10)
1220 std::ofstream writePhase(
"data/FrequencyError.dat",std::ios::app);
1222 if(writePhase.tellp() == 0)
1223 writePhase <<
"energy [MeV]" << std::setw(15)
1226 writePhase << E_m << std::setw(30) << std::setprecision(10)
1230 std::string energy = float2string(E_m);
1231 std::ofstream writeProperties(
"data/PropertiesForEnergy"+energy+
"MeV.dat", std::ios::out);
1232 writeProperties << std::left
1233 << std::setw(25) <<
"orbit radius"
1234 << std::setw(25) <<
"orbit momentum"
1235 << std::setw(25) <<
"inverse bending radius"
1236 << std::setw(25) <<
"field index"
1237 << std::setw(25) <<
"path length" <<
std::endl;
1239 for (
size_type i = 0; i < r_turn.size(); ++i) {
1240 writeProperties << std::setprecision(10) << std::left
1241 << std::setw(25) << r_turn[i]
1242 << std::setw(25) << peo[i]
1243 << std::setw(25) << h_turn[i]
1244 << std::setw(25) << fidx_turn[i]
1245 << std::setw(25) << ds_turn[i] <<
std::endl;
1250 writeAvgRadius.close();
1252 writeProperties.close();
value_type Emin_m
Minimum energy needed in cyclotron, .
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
bool converged_m
Is true if converged, false otherwise.
void decouple(const matrix_type &Mturn, sparse_matrix_type &R, sparse_matrix_type &invR)
Block diagonalizes the symplex part of the one turn transfer matrix.
boost::numeric::ublas::matrix< value_type > matrix_type
Type for storing maps.
size_type nSteps_m
Number of integration steps in total.
constexpr double e
The value of .
std::array< value_type, 3 > getEmittances() const
Returns the emittances (ex,ey,ez) in for which the code converged (since the whole simulation is don...
std::complex< value_type > complex_t
Finds a closed orbit of a cyclotron for a given energy.
bool match(value_type accuracy, size_type maxit, size_type maxitOrbit, Cyclotron *cycl, value_type denergy, value_type rguess, bool harmonic, bool full)
Searches for a matched distribution.
Hamiltonian H_m
Stores the Hamiltonian of the cyclotron.
Interface for a Cyclotron.
value_type getError() const
Returns the error (if the program didn't converged, one can call this function to check the error) ...
PETE_TUTree< FnFabs, typename T::PETE_Expr_t > fabs(const PETE_Expr< T > &l)
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
virtual double getPHIinit() const
The base class for all OPAL exceptions.
void initialize(value_type, value_type)
size_type N_m
Number of integration steps for closed orbit computation.
size_type nSectors_m
Number of (symmetric) sectors the field is averaged over.
std::array< value_type, 3 > emittance_m
Stores the desired emittances, .
value_type L2ErrorNorm(const matrix_type &, const matrix_type &)
Returns the L2-error norm between the old and new sigma-matrix.
matrix_type updateInitialSigma(const matrix_type &, sparse_matrix_type &, sparse_matrix_type &)
Computes the new initial sigma matrix.
value_type E_m
Stores the user-define energy, .
std::function< Series(value_type, value_type, value_type)> Hamiltonian
Type of the Hamiltonian for the cyclotron.
matrix_type sigma_m
Stores the stationary distribution (sigma matrix)
Value_type value_type
Type of variables.
This class generates the matrices for the one turn matrix of a cyclotron.
const double & getInjectionRadius() const
constexpr double alpha
The fine structure constant, no dimension.
Size_type size_type
Type for specifying sizes.
value_type gamma_m
Relativistic factor (which can be computed out ot the kinetic energy and rest mass (potential energy)...
void updateSigma(const std::vector< matrix_type > &, const std::vector< matrix_type > &)
Computes new sigma matrices (one for each angle)
size_type niterations_m
Is the number of iterations needed for convergence.
constexpr double pi
The value of .
value_type I_m
Stores the value of the current, .
boost::numeric::ublas::matrix< complex_t > cmatrix_type
Type for storing complex matrices.
boost::numeric::ublas::compressed_matrix< complex_t, boost::numeric::ublas::row_major > sparse_matrix_type
Type for storing the sparse maps.
constexpr double mu_0
The permeability of vacuum in Vs/Am.
value_type isEigenEllipse(const matrix_type &Mturn, const matrix_type &sigma)
Checks if the sigma-matrix is an eigenellipse and returns the L2 error.
value_type m_m
Is the mass of the particles, .
constexpr double c
The velocity of light in m/s.
const double & getInjectionMomentum() const
SigmaGenerator(value_type I, value_type ex, value_type ey, value_type ez, value_type E, value_type m, const Cyclotron *cycl, size_type N, size_type Nsectors, size_type truncOrder, bool write=true)
Constructs an object of type SigmaGenerator.
value_type L1ErrorNorm(const matrix_type &, const matrix_type &)
Returns the Lp-error norm between the old and new sigma-matrix.
void writeOrbitOutput_m(const std::pair< value_type, value_type > &tunes, const value_type &ravg, const value_type &freqError, const container_type &r_turn, const container_type &peo, const container_type &h_turn, const container_type &fidx_turn, const container_type &ds_turn)
Called within SigmaGenerator::match().
value_type error_m
Error of computation.
std::function< Series(value_type, value_type, value_type)> SpaceCharge
Type of the Hamiltonian for the space charge.
MMatrix< double > im(MMatrix< m_complex > mc)
value_type wo_m
Is the orbital frequency, .
T * value_type(const SliceIterator< T > &)
Tps< T > sqrt(const Tps< T > &x)
Square root.
boost::numeric::ublas::vector< value_type > vector_type
Type for storing vectors.
std::vector< value_type > container_type
Container for storing the properties for each angle.
bool invertMatrix_m(const sparse_matrix_type &R, sparse_matrix_type &invR)
const value_type const_value_type
Type of constant variables.
void eigsolve_m(const matrix_type &Mturn, sparse_matrix_type &R)
SpaceCharge Hsc_m
Stores the Hamiltonian for the space charge.
static void setGlobalTruncOrder(int order)
Set the global truncation order.
size_type getIterations() const
Returns the number of iterations needed for convergence (if not converged, it returns zero) ...
T::PETE_Expr_t::PETE_Return_t prod(const PETE_Expr< T > &expr)
value_type beta_m
Velocity (c/v), .
FTps< value_type, 2 *3 > Series
Type of the truncated powere series.
Series x_m
All variables x, px, y, py, z, delta.
constexpr double epsilon_0
The permittivity of vacuum in As/Vm.
size_type truncOrder_m
Truncation order of the power series for the Hamiltonian (cyclotron and space charge) ...
std::string float2string(value_type val)
Transforms a floating point value to a string.
value_type gamma2_m
Relativistic factor squared.
std::string::iterator iterator
bool for_each(const BareFieldIterator< T, D > &p, SameFieldID s, C)
size_type nStepsPerSector_m
Number of integration steps per sector (–> also: number of maps per sector)
bool write_m
Decides for writing output (default: true)
std::vector< matrix_type > sigmas_m
Vector storing the second moment matrix for each angle.
matrix_type & getSigma()
Returns the converged stationary distribution.
value_type nh_m
Harmonic number, .
Vector truncated power series in n variables.
static FTps makeVariable(int var)
Make variable.
value_type Emax_m
Maximum energy reached in cyclotron, .
Inform & endl(Inform &inf)
This class computes the matched distribution.
FVps< value_type, 2 *3 > Map
Type of a map.