Commit d533d65a authored by adelmann's avatar adelmann 🎗

Matched distribution working for the Ring

parent af2ff192
This diff is collapsed.
......@@ -159,6 +159,7 @@ namespace AttributesT
MAXSTEPSCO,
MAXSTEPSSI,
ORDERMAPS,
E2,
SIZE
};
}
......@@ -1707,24 +1708,10 @@ void Distribution::CreateDistributionFromFile(size_t numberOfParticles, double m
void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, double massIneV) {
/* ADA
/*
ToDo:
- fix M_m and wo
- if an energy range is specified, only calculate the tunes
- display tunes, r0, p0 etc.
- add flag in order to calculate tunes from FMLOWE to FMHIGHE
- eliminate physics and error
Note that you should use the following macros to access the Inform objects,
so that these messages may be left out at compile time:
INFOMSG("This is some information " << 34 << endl);
WARNMSG("This is a warning " << 34 << endl);
ERRORMSG("This is an error message " << 34 << endl);
DEBUGMSG("This is some debugging info " << 34 << endl);
*/
std::string LineName = Attributes::getString(itsAttr[AttributesT::LINE]);
......@@ -1756,17 +1743,27 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
<< " HN = " << CyclotronElement->getCyclHarm()
<< " INTSTEPS = " << Attributes::getReal(itsAttr[AttributesT::INTSTEPS]) << endl;
*gmsg << "----------------------------------------------------" << endl;
double wo = 8.44167E6*2.0*Physics::pi; // where is 8.44167E6 comming from
double M_m = 1; // units?
const double wo = CyclotronElement->getRfFrequ()*1E6/CyclotronElement->getCyclHarm()*2.0*Physics::pi;
const double fmLowE = CyclotronElement->getFMLowE();
const double fmHighE = CyclotronElement->getFMHighE();
if ((fmLowE>fmHighE) || (fmLowE<0) || (fmHighE<0)) {
ERRORMSG("FMHIGHE of FMLOWE not set propperly" << endl);
exit(1);
}
SigmaGenerator<double,unsigned int> siggen(I_m,
Attributes::getReal(itsAttr[AttributesT::EX])*1E6,
Attributes::getReal(itsAttr[AttributesT::EY])*1E6,
Attributes::getReal(itsAttr[AttributesT::EX])*1E6,
wo,
E_m*1E-6,
CyclotronElement->getCyclHarm(),M_m,72,590,
CyclotronElement->getCyclHarm(),
massIneV*1E-6,
fmLowE,
fmHighE,
CyclotronElement->getSymmetry(),
Attributes::getReal(itsAttr[AttributesT::INTSTEPS]),
Attributes::getString(itsAttr[AttributesT::FMAPFN]));
......@@ -1778,9 +1775,9 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
DEBUGMSG("Converged: Sigma-Matrix for " << E_m*1E-6 << " MeV" << endl);
for(unsigned int i=0; i<siggen.getSigma().size1(); ++i) {
for(unsigned int j=0; j<siggen.getSigma().size2(); ++j) {
*gmsg << std::setprecision(4) << siggen.getSigma()(i,j) << "\t";
INFOMSG(std::setprecision(4) << siggen.getSigma()(i,j) << "\t");
}
*gmsg << endl;
INFOMSG(endl);
}
/*
......@@ -3902,24 +3899,21 @@ void Distribution::SetAttributes() {
= Attributes::makeReal("ET", "Projected normalized emittance EY (m-rad) used to create matched distibution ", 1E-6);
itsAttr[AttributesT::INTSTEPS]
= Attributes::makeReal("INTSTEPS", "Integration steps used to create matched distibution ", 1440);
itsAttr[AttributesT::E2]
= Attributes::makeReal("E2", "If E2<Eb, we compute the tunes from the beams energy Eb to E2 with dE=0.25 MeV ", 0.0);
itsAttr[AttributesT::RESIDUUM]
= Attributes::makeReal("RESIDUUM", "Residuum for the closed orbit finder and sigma matrix generatro ", 1e-8);
itsAttr[AttributesT::MAXSTEPSCO]
= Attributes::makeReal("MAXSTEPSCO", "Maximum steps used to find closed orbit ", 100);
itsAttr[AttributesT::MAXSTEPSSI]
= Attributes::makeReal("MAXSTEPSSI", "Maximum steps used to find matched distribution ",1000);
itsAttr[AttributesT::ORDERMAPS]
= Attributes::makeReal("ORDERMAPS", "Oder used in the field expansion ", 7);
itsAttr[AttributesT::FNAME]
= Attributes::makeString("FNAME", "File for reading in 6D particle "
"coordinates.", "");
itsAttr[AttributesT::WRITETOFILE]
= Attributes::makeBool("WRITETOFILE", "Write initial distribution to file.",
false);
......
......@@ -10,6 +10,8 @@
#include <string>
#include <vector>
#include <Ippl.h>
#define MALLOC(n,a) ((a*)malloc((n)*sizeof(a)))
/// Temporary class for comparison with reference program of Dr. Christian Baumgarten.
......@@ -227,9 +229,9 @@ int MagneticField::ReadSectorMap(float** b,int nr,int nth,int nsc, std::string f
if (!err) j++;
}
fclose(f);
if (!err) fprintf(stderr,"Sector map %s read in.\n",fname.c_str());
if (!err) INFOMSG("Sector map " << fname << " read in" << endl);
return j;
} else fprintf(stderr,"Error opening %s.\n",fname.c_str());
} else ERRORMSG("Can not open " << fname << endl);
return 0;
}
......@@ -330,4 +332,4 @@ void MagneticField::MakeNFoldSymmetric(float** bmag, int N, int nr, int nth, int
}
}
#endif
\ No newline at end of file
#endif
This diff is collapsed.
This diff is collapsed.
/**
* @file math.h
* Mathematical functions that aren't part of the standard of C++
*
* @author Matthias Frey
* @version 1.0
*/
#ifndef MATH_H
#define MATH_H
#include <cmath>
/*!
* \addtogroup matt
* @{
*/
/// Defines additional functions
/// @brief Defines additional mathematical functions
namespace matt {
/// Computes the sign of the input argument
template<typename T>
T sign(const T val) {
return (val<0) ? T(-1) : T(1);
}
/// Computes the sign of the input argument
template<typename T>
T sign(const T val) {
return (std::signbit(val)) ? T(-1) : T(1);
}
}
/*! @} End of Doxygen Groups*/
#endif
\ No newline at end of file
#endif
/**
* @file matrix_vector_operation.h
* This file provides additional functions for the matrix classes of uBLAS that is part of the
* BOOST library (http://www.boost.org/).
*
* @author Matthias Frey
* @version 1.0
*/
#ifndef MATRIX_OPERATION_H
#define MATRIX_OPERATION_H
......@@ -14,90 +23,64 @@
* @{
*/
/// Expands the existing functions of the boost library uBLAS (http://www.boost.org/).
/// @brief Expands the existing functions of the boost library uBLAS (http://www.boost.org/).
namespace matt_boost {
using namespace boost::numeric::ublas;
/// Computes the trace of a square matrix
template<class V>
BOOST_UBLAS_INLINE
V trace(matrix<V> &e) {
V tr = 0;
if(e.size1()==e.size2()) {
matrix_vector_range<matrix<V> > diag(e,range(0,e.size1()),range(0,e.size2()));
tr = sum(diag);
}
else
Error::message("matt_boost::trace",Error::dim);
return tr;
}
/// Computes the cross product \f$ v_{1}\times v_{2}\f$ of two vectors in \f$ \mathbb{R}^{3} \f$
template<class V>
BOOST_UBLAS_INLINE
vector<V> cross_prod(vector<V> &v1, vector<V> &v2) {
vector<V> v(v1.size());
if(v1.size() == v2.size() && v1.size() == 3) {
v(0) = v1(1)*v2(2)-v1(2)*v2(1);
v(1) = v1(2)*v2(0)-v1(0)*v2(2);
v(2) = v1(0)*v2(1)-v1(1)*v2(0);
}
else
Error::message("matt_boost::cross_prod",Error::dim);
return v;
}
/// Computes Taylor-Series of M(s) = exp(F*s)
template<class V>
BOOST_UBLAS_INLINE
matrix<V> taylor_exp(const matrix<V>& F, const V ds, const unsigned int order) {
double fac = 1.0;
matrix<V> Fn = identity_matrix<V>(6);
matrix<V> M = Fn;
for(unsigned int k=1; k<order; ++k) {
fac *= ds/V(k);
Fn = prod(Fn,F);
M = M + fac*Fn;
}
return M;
}
/// Generalized matrix-matrix-matrix multiplication \f$ e_{1}\cdot e_{2}\cdot e_{3} \f$
template<class M, class E1, class E2, class E3>
BOOST_UBLAS_INLINE
M gemmm(const matrix_expression<E1>& e1, const matrix_expression<E2>& e2, const matrix_expression<E3>& e3) {
M tmp = prod(e2,e3);
return prod(e1,tmp);
}
// template<class V>
// BOOST_UBLAS_INLINE
// matrix<V> gemmm(const matrix<V>& A, const matrix<V>& B, const matrix<V>& C) {
// matrix<V> tmp(A.size1(),C.size2());
// if(A.size2() == B.size1() && B.size2() == C.size1()) {
// matrix<V> tmp1 = prod(A,B);
// tmp = prod(tmp1,C);
// } else
// Error::message("matt_boos::gemmm",Error::dim);
// return tmp;
// }
//
// /// Compute a matrix-matrix-matrix multiplication A*B*C
// template<class V>
// BOOST_UBLAS_INLINE
// matrix<V> gemmm(const compressed_matrix<V,row_major>& A, const matrix<V>& B, const compressed_matrix<V,row_major>& C) {
// matrix<V> tmp(A.size1(),C.size2());
// if(A.size2() == B.size1() && B.size2() == C.size1()) {
// matrix<V> tmp1 = prod(A,B);
// tmp = prod(tmp1,C);
// } else
// Error::message("matt_boost::gemmm",Error::dim);
// return tmp;
// }
using namespace boost::numeric::ublas;
/// Computes the trace of a square matrix
template<class V>
BOOST_UBLAS_INLINE
V trace(matrix<V>& e) {
V tr = 0;
if (e.size1() == e.size2()) {
matrix_vector_range<matrix<V> > diag(e,range(0,e.size1()),range(0,e.size2()));
tr = sum(diag);
}
else
Error::message("matt_boost::trace",Error::dim);
return tr;
}
/// Computes the cross product \f$ v_{1}\times v_{2}\f$ of two vectors in \f$ \mathbb{R}^{3} \f$
template<class V>
BOOST_UBLAS_INLINE
vector<V> cross_prod(vector<V>& v1, vector<V>& v2) {
vector<V> v(v1.size());
if (v1.size() == v2.size() && v1.size() == 3) {
v(0) = v1(1) * v2(2) - v1(2) * v2(1);
v(1) = v1(2) * v2(0) - v1(0) * v2(2);
v(2) = v1(0) * v2(1) - v1(1) * v2(0);
}
else
Error::message("matt_boost::cross_prod",Error::dim);
return v;
}
/// Computes Taylor-Series of M(s) = exp(F*s)
template<class V>
BOOST_UBLAS_INLINE
matrix<V> taylor_exp(const matrix<V>& F, const V ds, const unsigned int order) {
double fac = 1.0;
matrix<V> Fn = identity_matrix<V>(6);
matrix<V> M = Fn;
for (unsigned int k = 1; k < order; ++k) {
fac *= ds / V(k);
Fn = prod(Fn,F);
M = M + fac * Fn;
}
return M;
}
/// Generalized matrix-matrix-matrix multiplication \f$ e_{1}\cdot e_{2}\cdot e_{3} \f$
template<class M, class E1, class E2, class E3>
BOOST_UBLAS_INLINE
M gemmm(const matrix_expression<E1>& e1, const matrix_expression<E2>& e2, const matrix_expression<E3>& e3) {
M tmp = prod(e2,e3);
return prod(e1,tmp);
}
}
/*! @} End of Doxygen Groups*/
#endif
\ No newline at end of file
#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