Skip to content

GitLab

Commit 5c235b25 authored by snuverink_j's avatar snuverink_j
Browse files

for src/Distribution: reduce/fix indentation, spelling, small fixes

parent ca46c451
Hide whitespace changes
Inline Side-by-side
Showing with 161 additions and 166 deletions
src/Distribution/ClosedOrbitFinder.h
View file @ 5c235b25
......@@ -29,13 +29,8 @@
// #include "physics.h"
// #include "MagneticField.h" // ONLY FOR STAND-ALONE PROGRAM
#include "MagneticField.h"
#include <fstream>
// include headers for integration
#include <boost/numeric/odeint/integrate/integrate_n_steps.hpp>
......@@ -72,11 +67,11 @@ class ClosedOrbitFinder
* @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,
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 guess,
const std::string& type, value_type scaleFactor = 1.0,
bool = true);
bool domain = true);
/// Returns the inverse bending radius (size of container N+1)
container_type& getInverseBendingRadius();
......@@ -505,8 +500,9 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
ptheta = std::sqrt(p2 - pr2);
invptheta = 1.0 / ptheta;
// intepolate values of magnetic field
// interpolate values of magnetic field
bField_m.interpolate(bint, brint, btint, y[0], theta * 180.0 / Physics::pi);
bint *= invbcon;
brint *= invbcon;
......@@ -581,8 +577,8 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
// store initial values for updating values for higher energies
container_type previous_init = {0.0, 0.0};
do {
do {
// (re-)set inital values for r and pr
r_m[0] = init[0];
pr_m[0] = init[1];
......@@ -647,7 +643,7 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
invgamma4 = 1.0 / (gamma2 * gamma2);
} while (E != E_m);
} while (E != E_m);
/* store last entry, since it is needed in computeVerticalOscillations(), because we have to do the same
* number of integrations steps there.
......
src/Distribution/Distribution.cpp
View file @ 5c235b25
......@@ -53,7 +53,7 @@
extern Inform *gmsg;
#define SMALLESTCUTOFF 1e-12
constexpr double SMALLESTCUTOFF = 1e-12;
namespace {
SymTenzor<double, 6> getUnit6x6() {
......@@ -966,7 +966,7 @@ void Distribution::calcPartPerDist(size_t numberOfParticles) {
nPartFromFiles.insert(std::make_pair(i, nPart));
if (nPart > numberOfParticles) {
throw OpalException("Distribution::calcPartPerDist",
"Number of particles is to small");
"Number of particles is too small");
}
numberOfParticles -= nPart;
......@@ -1320,160 +1320,159 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
*/
std::string LineName = Attributes::getString(itsAttr[Attrib::Distribution::LINE]);
if (LineName != "") {
const BeamSequence* LineSequence = BeamSequence::find(LineName);
if (LineSequence != NULL) {
SpecificElementVisitor<Cyclotron> CyclotronVisitor(*LineSequence->fetchLine());
CyclotronVisitor.execute();
size_t NumberOfCyclotrons = CyclotronVisitor.size();
if (NumberOfCyclotrons > 1) {
throw OpalException("Distribution::createMatchedGaussDistribution",
"I am confused, found more than one Cyclotron element in line");
}
if (NumberOfCyclotrons == 0) {
throw OpalException("Distribution::createMatchedGaussDistribution",
"didn't find any Cyclotron element in line");
}
const Cyclotron* CyclotronElement = CyclotronVisitor.front();
*gmsg << "* ----------------------------------------------------" << endl;
*gmsg << "* About to find closed orbit and matched distribution " << endl;
*gmsg << "* I= " << I_m*1E3 << " (mA) E= " << E_m*1E-6 << " (MeV)" << endl;
*gmsg << "* EX= " << Attributes::getReal(itsAttr[Attrib::Distribution::EX])
<< "* EY= " << Attributes::getReal(itsAttr[Attrib::Distribution::EY])
<< "* ET= " << Attributes::getReal(itsAttr[Attrib::Distribution::ET])
<< "* FMAPFN " << Attributes::getString(itsAttr[Attrib::Distribution::FMAPFN]) << endl //CyclotronElement->getFieldMapFN() << endl
<< "* FMSYM= " << (int)Attributes::getReal(itsAttr[Attrib::Distribution::MAGSYM])
<< "* HN= " << CyclotronElement->getCyclHarm()
<< "* PHIINIT= " << CyclotronElement->getPHIinit() << endl;
*gmsg << "* ----------------------------------------------------" << endl;
const double wo = CyclotronElement->getRfFrequ()*1E6/CyclotronElement->getCyclHarm()*2.0*Physics::pi;
const double fmLowE = CyclotronElement->getFMLowE();
const double fmHighE = CyclotronElement->getFMHighE();
double lE,hE;
lE = fmLowE;
hE = fmHighE;
if ((lE<0) || (hE<0)) {
lE = E_m*1E-6;
hE = E_m*1E-6;
}
if (LineName == "") return;
int Nint = 1000;
bool writeMap = true;
typedef SigmaGenerator<double, unsigned int> sGenerator_t;
sGenerator_t *siggen = new sGenerator_t(I_m,
Attributes::getReal(itsAttr[Attrib::Distribution::EX])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::EY])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::ET])*1E6,
wo,
E_m*1E-6,
CyclotronElement->getCyclHarm(),
massIneV*1E-6,
lE,
hE,
(int)Attributes::getReal(itsAttr[Attrib::Distribution::MAGSYM]),
Nint,
Attributes::getString(itsAttr[Attrib::Distribution::FMAPFN]),
Attributes::getReal(itsAttr[Attrib::Distribution::ORDERMAPS]),
writeMap);
if (siggen->match(Attributes::getReal(itsAttr[Attrib::Distribution::RESIDUUM]),
Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSSI]),
Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSCO]),
CyclotronElement->getPHIinit(),
Attributes::getReal(itsAttr[Attrib::Distribution::RGUESS]),
Attributes::getString(itsAttr[Attrib::Distribution::FMTYPE]),
false)) {
std::array<double,3> Emit = siggen->getEmittances();
if (Attributes::getReal(itsAttr[Attrib::Distribution::RGUESS]) > 0)
*gmsg << "* RGUESS " << Attributes::getReal(itsAttr[Attrib::Distribution::RGUESS])/1000.0 << " (m) " << endl;
*gmsg << "* Converged (Ex, Ey, Ez) = (" << Emit[0] << ", " << Emit[1] << ", "
<< Emit[2] << ") pi mm mrad for E= " << E_m*1E-6 << " (MeV)" << endl;
*gmsg << "* Sigma-Matrix " << endl;
for (unsigned int i = 0; i < siggen->getSigma().size1(); ++ i) {
*gmsg << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,0);
for (unsigned int j = 1; j < siggen->getSigma().size2(); ++ j) {
*gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
}
*gmsg << " \\\\" << endl;
}
const BeamSequence* LineSequence = BeamSequence::find(LineName);
if (LineSequence == NULL)
throw OpalException("Distribution::CreateMatchedGaussDistribution",
"didn't find any Cyclotron element in line");
/*
SpecificElementVisitor<Cyclotron> CyclotronVisitor(*LineSequence->fetchLine());
CyclotronVisitor.execute();
size_t NumberOfCyclotrons = CyclotronVisitor.size();
Now setup the distribution generator
Units of the Sigma Matrix:
if (NumberOfCyclotrons > 1) {
throw OpalException("Distribution::createMatchedGaussDistribution",
"I am confused, found more than one Cyclotron element in line");
}
if (NumberOfCyclotrons == 0) {
throw OpalException("Distribution::createMatchedGaussDistribution",
"didn't find any Cyclotron element in line");
}
const Cyclotron* CyclotronElement = CyclotronVisitor.front();
*gmsg << "* ----------------------------------------------------" << endl;
*gmsg << "* About to find closed orbit and matched distribution " << endl;
*gmsg << "* I= " << I_m*1E3 << " (mA) E= " << E_m*1E-6 << " (MeV)" << endl;
*gmsg << "* EX= " << Attributes::getReal(itsAttr[Attrib::Distribution::EX])
<< " EY= " << Attributes::getReal(itsAttr[Attrib::Distribution::EY])
<< " ET= " << Attributes::getReal(itsAttr[Attrib::Distribution::ET]) << endl
<< "* FMAPFN= " << Attributes::getString(itsAttr[Attrib::Distribution::FMAPFN]) << endl //CyclotronElement->getFieldMapFN() << endl
<< "* FMSYM= " << (int)Attributes::getReal(itsAttr[Attrib::Distribution::MAGSYM])
<< " HN= " << CyclotronElement->getCyclHarm()
<< " PHIINIT= " << CyclotronElement->getPHIinit() << endl;
*gmsg << "* ----------------------------------------------------" << endl;
const double wo = CyclotronElement->getRfFrequ()*1E6/CyclotronElement->getCyclHarm()*2.0*Physics::pi;
const double fmLowE = CyclotronElement->getFMLowE();
const double fmHighE = CyclotronElement->getFMHighE();
double lE,hE;
lE = fmLowE;
hE = fmHighE;
if ((lE<0) || (hE<0)) {
lE = E_m*1E-6;
hE = E_m*1E-6;
}
spatial: mm
moment: rad
int Nint = 1000;
bool writeMap = true;
typedef SigmaGenerator<double, unsigned int> sGenerator_t;
sGenerator_t *siggen = new sGenerator_t(I_m,
Attributes::getReal(itsAttr[Attrib::Distribution::EX])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::EY])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::ET])*1E6,
wo,
E_m*1E-6,
CyclotronElement->getCyclHarm(),
massIneV*1E-6,
lE,
hE,
(int)Attributes::getReal(itsAttr[Attrib::Distribution::MAGSYM]),
Nint,
Attributes::getString(itsAttr[Attrib::Distribution::FMAPFN]),
Attributes::getReal(itsAttr[Attrib::Distribution::ORDERMAPS]),
writeMap);
if (siggen->match(Attributes::getReal(itsAttr[Attrib::Distribution::RESIDUUM]),
Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSSI]),
Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSCO]),
CyclotronElement->getPHIinit(),
Attributes::getReal(itsAttr[Attrib::Distribution::RGUESS]),
Attributes::getString(itsAttr[Attrib::Distribution::FMTYPE]),
false)) {
std::array<double,3> Emit = siggen->getEmittances();
if (Attributes::getReal(itsAttr[Attrib::Distribution::RGUESS]) > 0)
*gmsg << "* RGUESS " << Attributes::getReal(itsAttr[Attrib::Distribution::RGUESS])/1000.0 << " (m) " << endl;
*gmsg << "* Converged (Ex, Ey, Ez) = (" << Emit[0] << ", " << Emit[1] << ", "
<< Emit[2] << ") pi mm mrad for E= " << E_m*1E-6 << " (MeV)" << endl;
*gmsg << "* Sigma-Matrix " << endl;
for (unsigned int i = 0; i < siggen->getSigma().size1(); ++ i) {
*gmsg << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,0);
for (unsigned int j = 1; j < siggen->getSigma().size2(); ++ j) {
*gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
}
*gmsg << " \\\\" << endl;
}
*/
/*
if (Options::cloTuneOnly)
throw OpalException("Do only CLO and tune calculation","... ");
Now setup the distribution generator
Units of the Sigma Matrix:
spatial: mm
moment: rad
auto sigma = siggen->getSigma();
// change units from mm to m
for (unsigned int i = 0; i < 3; ++ i) {
for (unsigned int j = 0; j < 6; ++ j) {
sigma(2 * i, j) *= 1e-3;
sigma(j, 2 * i) *= 1e-3;
}
}
*/
for (unsigned int i = 0; i < 3; ++ i) {
if ( sigma(2 * i, 2 * i) < 0 || sigma(2 * i + 1, 2 * i + 1) < 0 )
throw OpalException("Distribution::CreateMatchedGaussDistribution()",
"Negative value on the diagonal of the sigma matrix.");
if (Options::cloTuneOnly)
throw OpalException("Do only CLO and tune calculation","... ");
sigmaR_m[i] = std::sqrt(sigma(2 * i, 2 * i));
sigmaP_m[i] = std::sqrt(sigma(2 * i + 1, 2 * i + 1));
}
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
for (unsigned int i = 0; i < 3; ++ i) {
sigmaP_m[i] = converteVToBetaGamma(sigmaP_m[i], massIneV);
}
}
auto sigma = siggen->getSigma();
// change units from mm to m
for (unsigned int i = 0; i < 3; ++ i) {
for (unsigned int j = 0; j < 6; ++ j) {
sigma(2 * i, j) *= 1e-3;
sigma(j, 2 * i) *= 1e-3;
}
}
correlationMatrix_m(1, 0) = sigma(0, 1) / (sqrt(sigma(0, 0) * sigma(1, 1)));
correlationMatrix_m(3, 2) = sigma(2, 3) / (sqrt(sigma(2, 2) * sigma(3, 3)));
correlationMatrix_m(5, 4) = sigma(4, 5) / (sqrt(sigma(4, 4) * sigma(5, 5)));
correlationMatrix_m(4, 0) = sigma(0, 4) / (sqrt(sigma(0, 0) * sigma(4, 4)));
correlationMatrix_m(4, 1) = sigma(1, 4) / (sqrt(sigma(1, 1) * sigma(4, 4)));
correlationMatrix_m(5, 0) = sigma(0, 5) / (sqrt(sigma(0, 0) * sigma(5, 5)));
correlationMatrix_m(5, 1) = sigma(1, 5) / (sqrt(sigma(1, 1) * sigma(5, 5)));
for (unsigned int i = 0; i < 3; ++ i) {
if ( sigma(2 * i, 2 * i) < 0 || sigma(2 * i + 1, 2 * i + 1) < 0 )
throw OpalException("Distribution::CreateMatchedGaussDistribution()",
"Negative value on the diagonal of the sigma matrix.");
createDistributionGauss(numberOfParticles, massIneV);
sigmaR_m[i] = std::sqrt(sigma(2 * i, 2 * i));
sigmaP_m[i] = std::sqrt(sigma(2 * i + 1, 2 * i + 1));
}
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
for (unsigned int i = 0; i < 3; ++ i) {
sigmaP_m[i] = converteVToBetaGamma(sigmaP_m[i], massIneV);
}
else {
*gmsg << "* Not converged for " << E_m*1E-6 << " MeV" << endl;
}
if (siggen)
delete siggen;
correlationMatrix_m(1, 0) = sigma(0, 1) / (sqrt(sigma(0, 0) * sigma(1, 1)));
correlationMatrix_m(3, 2) = sigma(2, 3) / (sqrt(sigma(2, 2) * sigma(3, 3)));
correlationMatrix_m(5, 4) = sigma(4, 5) / (sqrt(sigma(4, 4) * sigma(5, 5)));
correlationMatrix_m(4, 0) = sigma(0, 4) / (sqrt(sigma(0, 0) * sigma(4, 4)));
correlationMatrix_m(4, 1) = sigma(1, 4) / (sqrt(sigma(1, 1) * sigma(4, 4)));
correlationMatrix_m(5, 0) = sigma(0, 5) / (sqrt(sigma(0, 0) * sigma(5, 5)));
correlationMatrix_m(5, 1) = sigma(1, 5) / (sqrt(sigma(1, 1) * sigma(5, 5)));
throw OpalException("Distribution::CreateMatchedGaussDistribution",
"didn't find any matched distribution.");
}
createDistributionGauss(numberOfParticles, massIneV);
}
else {
*gmsg << "* Not converged for " << E_m*1E-6 << " MeV" << endl;
if (siggen)
delete siggen;
if (siggen)
delete siggen;
}
else
throw OpalException("Distribution::CreateMatchedGaussDistribution",
"didn't find any Cyclotron element in line");
throw OpalException("Distribution::CreateMatchedGaussDistribution",
"didn't find any matched distribution.");
}
if (siggen)
delete siggen;
}
void Distribution::createDistributionGauss(size_t numberOfParticles, double massIneV) {
......@@ -3556,15 +3555,15 @@ void Distribution::setAttributes() {
itsAttr[Attrib::Distribution::LINE]
= Attributes::makeString("LINE", "Beamline that contains a cyclotron or ring ", "");
itsAttr[Attrib::Distribution::FMAPFN]
= Attributes::makeString("FMAPFN", "File for reading fieldmap used to create matched distibution ", "");
= Attributes::makeString("FMAPFN", "File for reading fieldmap used to create matched distribution ", "");
itsAttr[Attrib::Distribution::FMTYPE]
= Attributes::makeString("FMTYPE", "File format for reading fieldmap used to create matched distribution ", "");
itsAttr[Attrib::Distribution::EX]
= Attributes::makeReal("EX", "Projected normalized emittance EX (m-rad), used to create matched distibution ", 1E-6);
= Attributes::makeReal("EX", "Projected normalized emittance EX (m-rad), used to create matched distribution ", 1E-6);
itsAttr[Attrib::Distribution::EY]
= Attributes::makeReal("EY", "Projected normalized emittance EY (m-rad) used to create matched distibution ", 1E-6);
= Attributes::makeReal("EY", "Projected normalized emittance EY (m-rad) used to create matched distribution ", 1E-6);
itsAttr[Attrib::Distribution::ET]
= Attributes::makeReal("ET", "Projected normalized emittance EY (m-rad) used to create matched distibution ", 1E-6);
= Attributes::makeReal("ET", "Projected normalized emittance ET (m-rad) used to create matched distribution ", 1E-6);
itsAttr[Attrib::Distribution::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[Attrib::Distribution::RESIDUUM]
......
src/Distribution/MagneticField.cpp
View file @ 5c235b25
......@@ -11,7 +11,7 @@ MagneticField::MagneticField(const std::string fmapfn,
}
void MagneticField::read(const std::string& type,
const double& scaleFactor) {
const double& scaleFactor) {
if ( type == "CARBONCYCL" )
this->getFieldFromFile_Carbon(scaleFactor);
......@@ -41,18 +41,18 @@ void MagneticField::read(const std::string& type,
void MagneticField::interpolate(double& bint,
double& brint,
double& btint,
const double& rad,
const double& theta
)
double& brint,
double& btint,
const double& rad,
const double& theta
)
{
// 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.
// ir : the number of path whose radius is less than the 4 points of cell which surround the particle.
const int ir = (int)xir;
// wr1 : the relative distance to the inner path radius
......@@ -77,7 +77,7 @@ void MagneticField::interpolate(double& bint,
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.
// it : the number of point on the inner path whose angle is less than the particle' corresponding angle.
// include zero degree point
it = it + 1;
......
src/Distribution/SigmaGenerator.h
View file @ 5c235b25
......@@ -120,8 +120,8 @@ class SigmaGenerator
* @param type specifies the magnetic field format (e.g. CARBONCYCL)
* @param harmonic is a boolean. If "true" the harmonics are used instead of the closed orbit finder.
*/
bool match(value_type, size_type, size_type, value_type,
value_type, const std::string&, bool);
bool match(value_type accuracy, size_type maxit, size_type maxitOrbit, value_type angle,
value_type guess, const std::string& type, bool harmonic);
/// 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>.
......@@ -133,7 +133,7 @@ class SigmaGenerator
* @param R is the 4x4 dimensional transformation matrix (gets computed)
* @param invR is the 4x4 dimensional inverse transformation (gets computed)
*/
vector_type decouple(const matrix_type&, sparse_matrix_type&, sparse_matrix_type&);
vector_type decouple(const matrix_type& Mturn, sparse_matrix_type& R, sparse_matrix_type& invR);
/// Checks if the sigma-matrix is an eigenellipse and returns the L2 error.
/*!
......@@ -141,7 +141,7 @@ class SigmaGenerator
* @param Mturn is the one turn transfer matrix
* @param sigma is the sigma matrix to be tested
*/
value_type isEigenEllipse(const matrix_type&, const matrix_type&);
value_type isEigenEllipse(const matrix_type& Mturn, const matrix_type& sigma);
/// Returns the converged stationary distribution
matrix_type& getSigma();
......
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