Commit ba29a483 authored by adelmann's avatar adelmann 🎗

latest changes

parent d712e9ed
This diff is collapsed.
......@@ -154,6 +154,7 @@ namespace AttributesT
EX, // below is for the matched distribution
EY,
ET,
MAGSYM, // number of sector magnets
LINE,
FMAPFN,
RESIDUUM,
......@@ -1739,9 +1740,10 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
*gmsg << "EX= " << Attributes::getReal(itsAttr[AttributesT::EX])
<< " EY= " << Attributes::getReal(itsAttr[AttributesT::EY])
<< " ET= " << Attributes::getReal(itsAttr[AttributesT::ET])
<< " FMAPFN " << CyclotronElement->getFieldMapFN()
<< " FMSYM= " << CyclotronElement->getSymmetry()
<< " HN = " << CyclotronElement->getCyclHarm() << endl;
<< " FMAPFN " << Attributes::getString(itsAttr[AttributesT::FMAPFN]) << endl //CyclotronElement->getFieldMapFN() << endl
<< "FMSYM= " << (int)Attributes::getReal(itsAttr[AttributesT::MAGSYM])
<< " HN= " << CyclotronElement->getCyclHarm()
<< " PHIINIT= " << CyclotronElement->getPHIinit() << endl;
*gmsg << "----------------------------------------------------" << endl;
const double wo = CyclotronElement->getRfFrequ()*1E6/CyclotronElement->getCyclHarm()*2.0*Physics::pi;
......@@ -1751,14 +1753,16 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
if ((fmLowE>fmHighE) || (fmLowE<0) || (fmHighE<0)) {
throw OpalException("Distribution::CreateMatchedGaussDistribution",
"FMHIGHE or FMLOWE not set properly");
"FMHIGHE or FMLOW not set properly");
}
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;
SigmaGenerator<double,unsigned int> siggen(I_m,
Attributes::getReal(itsAttr[AttributesT::EX])*1E6,
......@@ -1770,14 +1774,21 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
massIneV*1E-6,
fmLowE,
fmHighE,
nsc,
(int)Attributes::getReal(itsAttr[AttributesT::MAGSYM]),
rmin,
Nint,
nth,
Attributes::getString(itsAttr[AttributesT::FMAPFN]));
nr,
dr,
Attributes::getString(itsAttr[AttributesT::FMAPFN]),
Attributes::getReal(itsAttr[AttributesT::ORDERMAPS]),
writeMap);
if(siggen.match(Attributes::getReal(itsAttr[AttributesT::RESIDUUM]),
Attributes::getReal(itsAttr[AttributesT::MAXSTEPSCO]),
Attributes::getReal(itsAttr[AttributesT::MAXSTEPSSI]),
Attributes::getReal(itsAttr[AttributesT::ORDERMAPS]))) {
CyclotronElement->getPHIinit(),
false)) {
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) {
......@@ -1785,6 +1796,8 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
}
INFOMSG(endl);
}
}
/*
......@@ -1795,29 +1808,26 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
moment: rad
*/
sigmaR_m[0] = std::sqrt(1E-3*siggen.getSigma()(0,0)); // rms x
sigmaR_m[1] = std::sqrt(1E-3*siggen.getSigma()(2,2)); // rms y
sigmaR_m[2] = std::sqrt(1E-3*siggen.getSigma()(4,4)); // rms z
sigmaP_m[0] = std::sqrt(1E-3*siggen.getSigma()(1,1)); // UNITS ...
sigmaP_m[1] = std::sqrt(1E-3*siggen.getSigma()(3,3));
sigmaP_m[2] = std::sqrt(1E-3*siggen.getSigma()(5,5));
distCorr_m.push_back(1E-3*siggen.getSigma()(0,1)); // R12
distCorr_m.push_back(1E-3*siggen.getSigma()(2,3)); // R34
distCorr_m.push_back(1E-3*siggen.getSigma()(4,5)); // R56
distCorr_m.push_back(1E-3*siggen.getSigma()(0,5)); // R16
distCorr_m.push_back(1E-3*siggen.getSigma()(1,5)); // R26
distCorr_m.push_back(1E-3*siggen.getSigma()(0,4)); // R15
distCorr_m.push_back(1E-3*siggen.getSigma()(1,4)); // R25
} else {
*gmsg << "Not converged." << endl;
sigmaR_m[0] = std::sqrt(1E-3*siggen.getSigma()(0,0)); // rms x
sigmaR_m[1] = std::sqrt(1E-3*siggen.getSigma()(2,2)); // rms y
sigmaR_m[2] = std::sqrt(1E-3*siggen.getSigma()(4,4)); // rms z
sigmaP_m[0] = std::sqrt(1E-3*siggen.getSigma()(1,1)); // UNITS ...
sigmaP_m[1] = std::sqrt(1E-3*siggen.getSigma()(3,3));
sigmaP_m[2] = std::sqrt(1E-3*siggen.getSigma()(5,5));
distCorr_m.push_back(1E-3*siggen.getSigma()(0,1)); // R12
distCorr_m.push_back(1E-3*siggen.getSigma()(2,3)); // R34
distCorr_m.push_back(1E-3*siggen.getSigma()(4,5)); // R56
distCorr_m.push_back(1E-3*siggen.getSigma()(0,5)); // R16
distCorr_m.push_back(1E-3*siggen.getSigma()(1,5)); // R26
distCorr_m.push_back(1E-3*siggen.getSigma()(0,4)); // R15
distCorr_m.push_back(1E-3*siggen.getSigma()(1,4)); // R25
} else {
*gmsg << "Not converged." << endl;
}
CreateDistributionGauss(numberOfParticles, massIneV);
}
CreateDistributionGauss(numberOfParticles, massIneV);
}
else
throw OpalException("Distribution::CreateMatchedGaussDistribution",
......@@ -3912,7 +3922,13 @@ void Distribution::SetAttributes() {
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);
= Attributes::makeReal("ORDERMAPS", "Order used in the field expansion ", 7);
itsAttr[AttributesT::MAGSYM]
= Attributes::makeReal("MAGSYM", "Number of sector magnets ", 0);
itsAttr[AttributesT::FNAME]
= Attributes::makeString("FNAME", "File for reading in 6D particle "
"coordinates.", "");
......@@ -5051,4 +5067,4 @@ void Distribution::WriteOutFileInjection() {
reduce(numberOfParticles, numberOfParticles, OpAddAssign());
}
}
}
\ No newline at end of file
}
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -16,7 +16,7 @@
#include <boost/numeric/ublas/vector.hpp>
#include "error.h"
#include <stdexcept>
/*!
* \addtogroup matt_boost
......@@ -37,7 +37,8 @@ namespace matt_boost {
tr = sum(diag);
}
else
Error::message("matt_boost::trace",Error::dim);
throw std::length_error("Error in function trace() of matrix_vector_operation.h: Wrong matrix dimensions.");
return tr;
}
......@@ -52,7 +53,8 @@ namespace matt_boost {
v(2) = v1(0) * v2(1) - v1(1) * v2(0);
}
else
Error::message("matt_boost::cross_prod",Error::dim);
throw std::length_error("Error in function cross_prod() of matrix_vector_operation.h: Wrong vector dimensions.");
return v;
}
......
......@@ -11,30 +11,30 @@
/// Defines the typical constants in physics
namespace physics {
/// Vacuum permittivity \f$ \left[\frac{As}{Vm} = \frac{C^{2}}{Nm^{2}} = \frac{F}{m}\right] \f$ (deu: elektrische Feldkonstante)
long double eps0 = 8.854187817e-12;
/// Vacuum permeability \f$ \left[\frac{N}{A^{2}} = \frac{Vs}{Am}\right] \f$ (deu: magnetische Feldkonstante)
long double mu0 = 4.0e-7*M_PI;
/// Elementary charge \f$ \left[C\right] \f$
long double q0 = 1.0; //1.60217733e-19;
/// Speed of light \f$ \left[\frac{m}{s}\right] \f$
long double c = 2.99792458e8;
/// Proton mass \f$ \left[\frac{MeV}{c^{2}}\right] \f$
long double E0=938.272;
/// Cyclotron unit length \f$ \left[m\right] \f$
/*!
* This quantity is defined in the paper "Transverse-Longitudinal Coupling by Space Charge in Cyclotrons"
* of Dr. Christian Baumgarten (2012)
* The lambda function takes the orbital frequency \f$ \omega_{o} \f$ (also defined in paper) as input argument.
*/
std::function<double(double)> acon = [](double wo) { return c/wo; };
/// Cyclotron unit \f$ \left[T\right] \f$ (Tesla)
/*!
* The lambda function takes the orbital frequency \f$ \omega_{o} \f$ as input argument.
*/
std::function<double(double)> bcon = [](double wo) { return E0*1.0e7/(q0*c*acon(wo)); };
/// Vacuum permittivity \f$ \left[\frac{As}{Vm} = \frac{C^{2}}{Nm^{2}} = \frac{F}{m}\right] \f$ (deu: elektrische Feldkonstante)
long double eps0 = 8.854187817e-12;
/// Vacuum permeability \f$ \left[\frac{N}{A^{2}} = \frac{Vs}{Am}\right] \f$ (deu: magnetische Feldkonstante)
long double mu0 = 4.0e-7*M_PI;
/// Elementary charge \f$ \left[e\right] \f$
long double q0 = 1.0;
/// Speed of light \f$ \left[\frac{m}{s}\right] \f$
long double c = 2.99792458e8;
/// Kinetic energy of a proton \f$ \left[MeV\right] \f$
long double E0=938.272;
/// Cyclotron unit length \f$ \left[m\right] \f$
/*!
* This quantity is defined in the paper "Transverse-Longitudinal Coupling by Space Charge in Cyclotrons"
* of Dr. Christian Baumgarten (2012)
* The lambda function takes the orbital frequency \f$ \omega_{o} \f$ (also defined in paper) as input argument.
*/
std::function<double(double)> acon = [](double wo) { return c/wo; };
/// Cyclotron unit \f$ \left[T\right] \f$ (Tesla)
/*!
* The lambda function takes the orbital frequency \f$ \omega_{o} \f$ as input argument.
*/
std::function<double(double)> bcon = [](double wo) { return E0*1.0e7/(q0*c*acon(wo)); };
};
/*! @} End of Doxygen Groups*/
#endif
\ No newline at end of file
#endif
......@@ -12,13 +12,12 @@
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include "error.h"
#include "physical_error.h"
#include "matrix_vector_operation.h"
// Remark: Comments starting with "///" are doxygen comments, comments with "//" are normal comments.
......@@ -38,7 +37,10 @@ class RDM
typedef boost::numeric::ublas::matrix<value_type> matrix_type;
/// Dense vector type definition
typedef boost::numeric::ublas::vector<value_type> vector_type;
/// Default constructor (sets only NumOfRDMs and DimOfRDMs)
RDM();
/// Returns the i-th Real Dirac matrix
/*!
* @param i specifying the matrix (has to be in the range from 0 to 15)
......@@ -94,17 +96,13 @@ class RDM
void transform(matrix_type&, short, value_type, sparse_matrix_type&, sparse_matrix_type&);
};
// -----------------------------------------------------------------------------------------------------------------------
// PUBLIC MEMBER VARIABLES
// -----------------------------------------------------------------------------------------------------------------------
short NumOfRDMs = 16;
short DimOfRDMs = 4;
// -----------------------------------------------------------------------------------------------------------------------
// PUBLIC MEMBER FUNCTIONS
// -----------------------------------------------------------------------------------------------------------------------
template<typename Value_type, typename Size_type>
RDM<Value_type, Size_type>::RDM() : NumOfRDMs(16), DimOfRDMs(4) {};
template<typename Value_type, typename Size_type>
typename RDM<Value_type, Size_type>::sparse_matrix_type RDM<Value_type, Size_type>::getRDM(short i) {
sparse_matrix_type rdm(4,4,4); // #nrows, #ncols, #non-zeros
......@@ -125,7 +123,7 @@ typename RDM<Value_type, Size_type>::sparse_matrix_type RDM<Value_type, Size_typ
case 13: rdm(0,3) = rdm(3,0) = -1; rdm(1,2) = rdm(2,1) = 1; break;
case 14: rdm(0,3) = rdm(1,2) = -1; rdm(2,1) = rdm(3,0) = 1; break;
case 15: rdm(0,0) = rdm(1,1) = rdm(2,2) = rdm(3,3) = 1; break;
default: Error::message("RDM<Value_type, Size_type>::generateRDM(short& i)",Error::notdefined); break;
default: throw std::out_of_range("Error in RDM::generate: 0 <= i <= 15"); break;
}
return rdm;
}
......@@ -136,8 +134,8 @@ typename RDM<Value_type, Size_type>::vector_type RDM<Value_type, Size_type>::dec
* Formula (11) from paper:
* Geometrical method of decoupling
*/
if(M.size1() != 4 && M.size1() != M.size2())
Error::message("RDM<Value_type, Size_type>::decompose(const matrix_type&)", Error::dim);
if (M.size1() != 4 && M.size1() != M.size2())
throw std::length_error("Error in RDM::decompose: Wrong matrix dimensions.");
vector_type coeffs(16);
......@@ -162,7 +160,7 @@ typename RDM<Value_type, Size_type>::vector_type RDM<Value_type, Size_type>::dec
template<typename Value_type, typename Size_type>
typename RDM<Value_type, Size_type>::matrix_type RDM<Value_type, Size_type>::combine(const vector_type& coeffs) {
if (coeffs.size() > 16)
Error::message("RDM<Value_type, Size_type>::combine(const vector_type&)",Error::size);
throw std::length_error("Error in RDM::combine: Wrong size of coefficient vector.");
// initialize a 4x4 zero matrix
matrix_type M = boost::numeric::ublas::zero_matrix<value_type>(4,4);
......@@ -227,7 +225,9 @@ void RDM<Value_type, Size_type>::diagonalize(matrix_type& Ms, sparse_matrix_type
E(0) = mult(4); E(1) = mult(5); E(2) = mult(6);
B(0) = -mult(7); B(1) = -mult(8); B(2) = -mult(9);
mr = boost::numeric::ublas::inner_prod(E,B);
b = eps*B+matt_boost::cross_prod(E,P);
b = eps*B+matt_boost::cross_prod(E,P);
// Transformation distinction made according to function "rdm_Decouple_F" in rdm.c of Dr. Christian Baumgarten
if (std::fabs(mr) < std::fabs(b(1))) {
transform(Ms,short(2),0.5 * std::atanh(mr / b(1)),R,invR);
......@@ -248,9 +248,10 @@ void RDM<Value_type, Size_type>::diagonalize(matrix_type& Ms, sparse_matrix_type
value_type P2 = boost::numeric::ublas::inner_prod(P,P);
value_type E2 = boost::numeric::ublas::inner_prod(E,E);
// 5. Transform with \gamma_{0}
transform(Ms,short(0),0.25 * std::atan2(mb,0.5 * (E2 - P2)),R,invR);
// 5. Transformation with \gamma_{8}
// 6. Transformation with \gamma_{8}
P(0) = mult(1); P(2) = mult(3);
transform(Ms,short(8),-0.5 * std::atan2(P(2),P(0)),R,invR);
......
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