Code indexing in gitaly is broken and leads to code not being visible to the user. We work on the issue with highest priority.

Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • OPAL/src
  • zheng_d/src
  • ext-rogers_c/src
  • ext-wang_c/src
  • cortes_c/src
  • ext-calvo_p/src
  • ext-edelen_a/src
  • albajacas_a/src
  • kraus/src
  • snuverink_j/OPAL-src
  • adelmann/src
  • muralikrishnan/src
  • wyssling_t/src
  • gsell/src
  • ext-piot_p/src
  • OPAL/opal-src-4-opalx-debug
  • winkle_m/src
17 results
Show changes
Commits on Source (45)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1659 additions and 1823 deletions
src/Algorithms/CavityAutophaser.cpp
View file @ 05397703
......@@ -72,7 +72,7 @@ double CavityAutophaser::getPhaseAtMaxEnergy(const Vector_t &R,
<< std::left << std::setw(68) << std::setfill('*') << ss.str()
<< std::setfill(' ') << endl);
if (!apVeto) {
double initialEnergy = Util::getEnergy(P, itsReference_m.getM()) * 1e-6;
double initialEnergy = Util::getKineticEnergy(P, itsReference_m.getM()) * 1e-6;
double AstraPhase = 0.0;
double designEnergy = element->getDesignEnergy();
......@@ -191,7 +191,7 @@ double CavityAutophaser::guessCavityPhase(double t) {
return orig_phi;
}
Phimax = element->getAutoPhaseEstimate(getEnergyMeV(refP),
Phimax = element->getAutoPhaseEstimate(Util::getKineticEnergy(refP, itsReference_m.getM()) * 1e-6,
t,
itsReference_m.getQ(),
itsReference_m.getM() * 1e-6);
......@@ -288,7 +288,7 @@ double CavityAutophaser::track(Vector_t /*R*/,
out);
rfc->setPhasem(initialPhase);
double finalKineticEnergy = Util::getEnergy(Vector_t(0.0, 0.0, pe.first), itsReference_m.getM() * 1e-6);
double finalKineticEnergy = Util::getKineticEnergy(Vector_t(0.0, 0.0, pe.first), itsReference_m.getM() * 1e-6);
return finalKineticEnergy;
}
\ No newline at end of file
src/Algorithms/CavityAutophaser.h
View file @ 05397703
......@@ -28,8 +28,6 @@ private:
const double dt,
const double phase,
std::ofstream *out = NULL) const;
double getEnergyMeV(const Vector_t &P);
double getMomentum(double kineticEnergyMeV);
const PartData &itsReference_m;
std::shared_ptr<Component> itsCavity_m;
......@@ -39,14 +37,4 @@ private:
};
inline
double CavityAutophaser::getEnergyMeV(const Vector_t &P) {
return itsReference_m.getM() * 1e-6 * (sqrt(dot(P,P) + 1) - 1);
}
inline
double CavityAutophaser::getMomentum(double kineticEnergyMeV) {
return sqrt(std::pow(kineticEnergyMeV / (itsReference_m.getM() * 1e-6) + 1, 2) - 1);
}
#endif
\ No newline at end of file
src/Algorithms/ParallelTTracker.cpp
View file @ 05397703
......@@ -1162,7 +1162,7 @@ void ParallelTTracker::findStartPosition(const BorisPusher &pusher) {
selectDT();
if ((back_track && itsOpalBeamline_m.containsSource()) ||
Util::getEnergy(itsBunch_m->RefPartP_m, itsBunch_m->getM()) < 1e-3) {
Util::getKineticEnergy(itsBunch_m->RefPartP_m, itsBunch_m->getM()) < 1e-3) {
double gamma = 0.1 / itsBunch_m->getM() + 1.0;
double beta = sqrt(1.0 - 1.0 / std::pow(gamma, 2));
itsBunch_m->RefPartP_m = itsBunch_m->toLabTrafo_m.rotateTo(beta * gamma * Vector_t(0, 0, 1));
......
src/Classic/AbsBeamline/RFCavity.cpp
View file @ 05397703
......@@ -494,22 +494,22 @@ ElementBase::ElementType RFCavity::getType() const {
double RFCavity::getAutoPhaseEstimateFallback(double E0, double t0, double q, double mass) {
const double dt = 1e-13;
const double p0 = Util::getP(E0, mass);
const double p0 = Util::getBetaGamma(E0, mass);
const double origPhase =getPhasem();
double dphi = Physics::pi / 18;
double phi = 0.0;
setPhasem(phi);
std::pair<double, double> ret = trackOnAxisParticle(E0 / mass, t0, dt, q, mass);
std::pair<double, double> ret = trackOnAxisParticle(p0, t0, dt, q, mass);
double phimax = 0.0;
double Emax = Util::getEnergy(Vector_t(0.0, 0.0, ret.first), mass);
double Emax = Util::getKineticEnergy(Vector_t(0.0, 0.0, ret.first), mass);
phi += dphi;
for (unsigned int j = 0; j < 2; ++ j) {
for (unsigned int i = 0; i < 36; ++ i, phi += dphi) {
setPhasem(phi);
ret = trackOnAxisParticle(p0, t0, dt, q, mass);
double Ekin = Util::getEnergy(Vector_t(0.0, 0.0, ret.first), mass);
double Ekin = Util::getKineticEnergy(Vector_t(0.0, 0.0, ret.first), mass);
if (Ekin > Emax) {
Emax = Ekin;
phimax = phi;
......@@ -635,7 +635,7 @@ double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const
}
double cosine_part = 0.0, sine_part = 0.0;
double p0 = std::sqrt((E0 / mass + 1) * (E0 / mass + 1) - 1);
double p0 = Util::getBetaGamma(E0, mass);
cosine_part += scale_m * std::cos(frequency_m * t0) * F[0];
sine_part += scale_m * std::sin(frequency_m * t0) * F[0];
......@@ -675,11 +675,11 @@ std::pair<double, double> RFCavity::trackOnAxisParticle(const double &p0,
const double zend = length_m + startField_m;
Vector_t z(0.0, 0.0, zbegin);
double dz = 0.5 * p(2) / std::sqrt(1.0 + dot(p, p)) * cdt;
double dz = 0.5 * p(2) / Util::getGamma(p) * cdt;
Vector_t Ef(0.0), Bf(0.0);
if (out) *out << std::setw(18) << z[2]
<< std::setw(18) << Util::getEnergy(p, mass)
<< std::setw(18) << Util::getKineticEnergy(p, mass)
<< std::endl;
while (z(2) + dz < zend && z(2) + dz > zbegin) {
z /= cdt;
......@@ -700,7 +700,7 @@ std::pair<double, double> RFCavity::trackOnAxisParticle(const double &p0,
t += dt;
if (out) *out << std::setw(18) << z[2]
<< std::setw(18) << Util::getEnergy(p, mass)
<< std::setw(18) << Util::getKineticEnergy(p, mass)
<< std::endl;
}
......
src/Classic/Algorithms/PartBunch.cpp
View file @ 05397703
......@@ -123,8 +123,7 @@ void PartBunch::computeSelfFields(int binNumber) {
if(fs_m->hasValidSolver()) {
/// Mesh the whole domain
if(fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
resizeMesh();
/// Scatter charge onto space charge grid.
this->Q *= this->dt;
......@@ -328,12 +327,14 @@ void PartBunch::computeSelfFields(int binNumber) {
}
void PartBunch::resizeMesh() {
if (fs_m->getFieldSolverType() != "SAAMG") {
return;
}
double xmin = fs_m->solver_m->getXRangeMin();
double xmax = fs_m->solver_m->getXRangeMax();
double ymin = fs_m->solver_m->getYRangeMin();
double ymax = fs_m->solver_m->getYRangeMax();
double zmin = fs_m->solver_m->getZRangeMin();
double zmax = fs_m->solver_m->getZRangeMax();
if(xmin > rmin_m[0] || xmax < rmax_m[0] ||
ymin > rmin_m[1] || ymax < rmax_m[1]) {
......@@ -354,14 +355,14 @@ void PartBunch::resizeMesh() {
boundp();
get_bounds(rmin_m, rmax_m);
}
Vector_t mymin = Vector_t(xmin, ymin , zmin);
Vector_t mymax = Vector_t(xmax, ymax , zmax);
for (int i = 0; i < 3; i++)
hr_m[i] = (mymax[i] - mymin[i])/nr_m[i];
Vector_t origin = Vector_t(0.0, 0.0, 0.0);
// update the mesh origin and mesh spacing hr_m
fs_m->solver_m->resizeMesh(origin, hr_m, rmin_m, rmax_m, dh_m);
getMesh().set_meshSpacing(&(hr_m[0]));
getMesh().set_origin(mymin);
getMesh().set_origin(origin);
rho_m.initialize(getMesh(),
getFieldLayout(),
......@@ -384,8 +385,7 @@ void PartBunch::computeSelfFields() {
if(fs_m->hasValidSolver()) {
//mesh the whole domain
if(fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
resizeMesh();
//scatter charges onto grid
this->Q *= this->dt;
......@@ -510,8 +510,7 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
if(fs_m->hasValidSolver()) {
/// mesh the whole domain
if(fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
resizeMesh();
/// scatter particles charge onto grid.
this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
......@@ -639,8 +638,7 @@ void PartBunch::computeSelfFields_cycl(int bin) {
if(fs_m->hasValidSolver()) {
/// mesh the whole domain
if(fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
resizeMesh();
/// scatter particles charge onto grid.
this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
......
src/Classic/Algorithms/PartBunchBase.h
View file @ 05397703
......@@ -411,6 +411,8 @@ public:
// virtual void setFieldLayout(FieldLayout_t* fLayout) = 0;
virtual FieldLayout_t &getFieldLayout() = 0;
virtual void resizeMesh() { };
/*
* Wrapped member functions of IpplParticleBase
*/
......
src/Classic/Utilities/Util.h
View file @ 05397703
......@@ -6,6 +6,7 @@
#include <string>
#include <cstring>
#include <limits>
#include <sstream>
#include <type_traits>
#include <functional>
......@@ -28,16 +29,23 @@ namespace Util {
}
inline
double getEnergy(Vector_t p, double mass) {
double getKineticEnergy(Vector_t p, double mass) {
return (getGamma(p) - 1.0) * mass;
}
inline
double getP(double E, double mass) {
double gamma = E / mass + 1;
return std::sqrt(std::pow(gamma, 2.0) - 1.0);
double getBetaGamma(double Ekin, double mass) {
double value = std::sqrt(std::pow(Ekin / mass + 1.0, 2) - 1.0);
if (value < std::numeric_limits<double>::epsilon())
value = std::sqrt(2 * Ekin / mass);
return value;
}
inline
double convertMomentumeVToBetaGamma(double p, double mass) {
return p / mass;
}
inline
std::string getTimeString(double time, unsigned int precision = 3) {
std::string timeUnit(" [ps]");
......@@ -162,7 +170,7 @@ namespace Util {
}
Vector_t getTaitBryantAngles(Quaternion rotation, const std::string &elementName = "");
std::string toUpper(const std::string &str);
std::string combineFilePath(std::initializer_list<std::string>);
......
src/Distribution/CMakeLists.txt
View file @ 05397703
set (_SRCS
Distribution.cpp
LaserProfile.cpp
RealDiracMatrix.cpp
SigmaGenerator.cpp
)
include_directories (
......@@ -12,10 +14,10 @@ add_opal_sources(${_SRCS})
set (HDRS
ClosedOrbitFinder.h
Distribution.h
Harmonics.h
LaserProfile.h
MapGenerator.h
matrix_vector_operation.h
RealDiracMatrix.h
SigmaGenerator.h
)
......
src/Distribution/ClosedOrbitFinder.h
View file @ 05397703
......@@ -40,6 +40,7 @@
#include "Utilities/Options.h"
#include "Utilities/Options.h"
#include "Utilities/OpalException.h"
#include "Physics/Physics.h"
#include "AbstractObjects/OpalData.h"
......
src/Distribution/Distribution.cpp
View file @ 05397703
......@@ -7,7 +7,7 @@
// OPAL is licensed under GNU GPL version 3.
#include "Distribution/Distribution.h"
#include "Distribution/SigmaGenerator.h"
#include "Distribution/ClosedOrbitFinder.h"
#include "AbsBeamline/SpecificElementVisitor.h"
#include <cmath>
......@@ -17,7 +17,6 @@
#include <string>
#include <vector>
#include <numeric>
#include <limits>
// IPPL
#include "DataSource/DataConnect.h"
......@@ -212,16 +211,6 @@ Distribution::~Distribution() {
delete laserProfile_m;
}
/*
void Distribution::printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out) {
for (int i=0; i<M.size1(); ++i) {
for (int j=0; j<M.size2(); ++j) {
*gmsg << M(i,j) << " ";
}
*gmsg << endl;
}
}
*/
/**
* Calculate the local number of particles evenly and adjust node 0
......@@ -885,14 +874,6 @@ void Distribution::chooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUn
}
double Distribution::converteVToBetaGamma(double valueIneV, double massIneV) {
double value = std::copysign(std::sqrt(std::pow(std::abs(valueIneV) / massIneV + 1.0, 2) - 1.0), valueIneV);
if (std::abs(value) < std::numeric_limits<double>::epsilon())
value = std::copysign(std::sqrt(2 * std::abs(valueIneV) / massIneV), valueIneV);
return value;
}
void Distribution::createDistributionBinomial(size_t numberOfParticles, double massIneV) {
setDistParametersBinomial(massIneV);
......@@ -1079,9 +1060,9 @@ void Distribution::createDistributionFromFile(size_t /*numberOfParticles*/, doub
saveProcessor = 0;
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
P(0) = converteVToBetaGamma(P(0), massIneV);
P(1) = converteVToBetaGamma(P(1), massIneV);
P(2) = converteVToBetaGamma(P(2), massIneV);
P(0) = Util::convertMomentumeVToBetaGamma(P(0), massIneV);
P(1) = Util::convertMomentumeVToBetaGamma(P(1), massIneV);
P(2) = Util::convertMomentumeVToBetaGamma(P(2), massIneV);
}
pmean_m += P;
......@@ -1278,18 +1259,18 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
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,
E_m*1E-6,
massIneV*1E-6,
CyclotronElement,
Nint,
Nsectors,
Attributes::getReal(itsAttr[Attrib::Distribution::ORDERMAPS]),
writeMap);
std::unique_ptr<SigmaGenerator> siggen = std::unique_ptr<SigmaGenerator>(
new SigmaGenerator(I_m,
Attributes::getReal(itsAttr[Attrib::Distribution::EX])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::EY])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::ET])*1E6,
E_m*1E-6,
massIneV*1E-6,
CyclotronElement,
Nint,
Nsectors,
Attributes::getReal(itsAttr[Attrib::Distribution::ORDERMAPS]),
writeMap));
if (siggen->match(accuracy,
Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSSI]),
......@@ -1297,7 +1278,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
CyclotronElement,
denergy,
rguess,
false, full)) {
full)) {
std::array<double,3> Emit = siggen->getEmittances();
......@@ -1311,62 +1292,18 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
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) {
if (siggen->getSigma()(i,j) < 10e-12){
if (std::abs(siggen->getSigma()(i,j)) < 1.0e-15) {
*gmsg << " & " << std::setprecision(4) << std::setw(11) << 0.0;
}
else{
*gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
*gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
}
}
*gmsg << " \\\\" << endl;
}
/*
Now setup the distribution generator
Units of the Sigma Matrix:
spatial: mm
moment: rad
*/
double gamma = E_m / massIneV + 1.0;
double beta = std::sqrt(1.0 - 1.0 / (gamma * gamma));
auto sigma = siggen->getSigma();
// change units from mm to m
for (unsigned int i = 0; i < 6; ++ i)
for (unsigned int j = 0; j < 6; ++ j) sigma(i, j) *= 1e-6;
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.");
}
sigmaR_m[0] = std::sqrt(sigma(0, 0));
sigmaP_m[0] = std::sqrt(sigma(1, 1))*beta*gamma;
sigmaR_m[2] = std::sqrt(sigma(2, 2));
sigmaP_m[2] = std::sqrt(sigma(3, 3))*beta*gamma;
sigmaR_m[1] = std::sqrt(sigma(4, 4));
//p_l^2 = [(delta+1)*beta*gamma]^2 - px^2 - pz^2
double pl2 = (std::sqrt(sigma(5,5)) + 1)*(std::sqrt(sigma(5,5)) + 1)*beta*gamma*beta*gamma -
sigmaP_m[0]*sigmaP_m[0] - sigmaP_m[2]*sigmaP_m[2];
double pl = std::sqrt(pl2);
sigmaP_m[1] = gamma*(pl - beta*gamma);
correlationMatrix_m(1, 0) = sigma(0, 1) / (std::sqrt(sigma(0, 0) * sigma(1, 1)));
correlationMatrix_m(3, 2) = sigma(2, 3) / (std::sqrt(sigma(2, 2) * sigma(3, 3)));
correlationMatrix_m(5, 4) = sigma(4, 5) / (std::sqrt(sigma(4, 4) * sigma(5, 5)));
correlationMatrix_m(4, 0) = sigma(0, 4) / (std::sqrt(sigma(0, 0) * sigma(4, 4)));
correlationMatrix_m(4, 1) = sigma(1, 4) / (std::sqrt(sigma(1, 1) * sigma(4, 4)));
correlationMatrix_m(5, 0) = sigma(0, 5) / (std::sqrt(sigma(0, 0) * sigma(5, 5)));
correlationMatrix_m(5, 1) = sigma(1, 5) / (std::sqrt(sigma(1, 1) * sigma(5, 5)));
createDistributionGauss(numberOfParticles, massIneV);
generateMatchedGauss(siggen->getSigma(), numberOfParticles, massIneV);
// update injection radius and radial momentum
CyclotronElement->setRinit(siggen->getInjectionRadius() * 1.0e3);
......@@ -1375,13 +1312,9 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
else {
*gmsg << "* Not converged for " << E_m*1E-6 << " MeV" << endl;
delete siggen;
throw OpalException("Distribution::CreateMatchedGaussDistribution",
"didn't find any matched distribution.");
}
delete siggen;
}
void Distribution::createDistributionGauss(size_t numberOfParticles, double massIneV) {
......@@ -1480,7 +1413,7 @@ void Distribution::createOpalT(PartBunchBase<double, 3> *beam,
// This is PC from BEAM
double deltaP = Attributes::getReal(itsAttr[Attrib::Distribution::OFFSETP]);
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
deltaP = converteVToBetaGamma(deltaP, beam->getM());
deltaP = Util::convertMomentumeVToBetaGamma(deltaP, beam->getM());
}
avrgpz_m = beam->getP()/beam->getM() + deltaP;
......@@ -2421,6 +2354,136 @@ void Distribution::generateGaussZ(size_t numberOfParticles) {
gsl_matrix_free(corMat);
}
void Distribution::generateMatchedGauss(const SigmaGenerator::matrix_t& sigma,
size_t numberOfParticles, double massIneV)
{
/* This particle distribution generation is based on a
* symplectic method described in
* https://arxiv.org/abs/1205.3601
*/
/* Units of the Sigma Matrix:
* spatial: m
* moment: rad
*
* Attention: The vertical and longitudinal direction must be interchanged!
*/
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::generateMatchedGauss()",
"Negative value on the diagonal of the sigma matrix.");
}
double bgam = Util::getBetaGamma(E_m, massIneV);
/*
* only used for printing
*/
// horizontal
sigmaR_m[0] = std::sqrt(sigma(0, 0));
sigmaP_m[0] = std::sqrt(sigma(1, 1)) * bgam;
// longitudinal
sigmaR_m[1] = std::sqrt(sigma(4, 4));
sigmaP_m[1] = std::sqrt(sigma(5, 5)) * bgam;
// vertical
sigmaR_m[2] = std::sqrt(sigma(2, 2));
sigmaP_m[2] = std::sqrt(sigma(3, 3)) * bgam;
correlationMatrix_m(1, 0) = sigma(0, 1) / (std::sqrt(sigma(0, 0) * sigma(1, 1)));
correlationMatrix_m(3, 2) = sigma(2, 3) / (std::sqrt(sigma(2, 2) * sigma(3, 3)));
correlationMatrix_m(5, 4) = sigma(4, 5) / (std::sqrt(sigma(4, 4) * sigma(5, 5)));
correlationMatrix_m(4, 0) = sigma(0, 4) / (std::sqrt(sigma(0, 0) * sigma(4, 4)));
correlationMatrix_m(4, 1) = sigma(1, 4) / (std::sqrt(sigma(1, 1) * sigma(4, 4)));
correlationMatrix_m(5, 0) = sigma(0, 5) / (std::sqrt(sigma(0, 0) * sigma(5, 5)));
correlationMatrix_m(5, 1) = sigma(1, 5) / (std::sqrt(sigma(1, 1) * sigma(5, 5)));
inputMoUnits_m = InputMomentumUnitsT::NONE;
/*
* decouple horizontal and longitudinal direction
*/
// extract horizontal and longitudinal directions
RealDiracMatrix::matrix_t A(4, 4);
A(0, 0) = sigma(0, 0);
A(1, 1) = sigma(1, 1);
A(2, 2) = sigma(4, 4);
A(3, 3) = sigma(5, 5);
A(0, 1) = sigma(0, 1);
A(0, 2) = sigma(0, 4);
A(0, 3) = sigma(0, 5);
A(1, 0) = sigma(1, 0);
A(2, 0) = sigma(4, 0);
A(3, 0) = sigma(5, 0);
A(1, 2) = sigma(1, 4);
A(2, 1) = sigma(4, 1);
A(1, 3) = sigma(1, 5);
A(3, 1) = sigma(5, 1);
A(2, 3) = sigma(4, 5);
A(3, 2) = sigma(5, 4);
RealDiracMatrix rdm;
RealDiracMatrix::sparse_matrix_t R1 = rdm.diagonalize(A);
RealDiracMatrix::vector_t variances(8);
for (int i = 0; i < 4; ++i) {
variances(i) = std::sqrt(A(i, i));
}
/*
* decouple vertical direction
*/
A *= 0.0;
A(0, 0) = 1;
A(1, 1) = 1;
A(2, 2) = sigma(2, 2);
A(3, 3) = sigma(3, 3);
A(2, 3) = sigma(2, 3);
A(3, 2) = sigma(3, 2);
RealDiracMatrix::sparse_matrix_t R2 = rdm.diagonalize(A);
for (int i = 0; i < 4; ++i) {
variances(4 + i) = std::sqrt(A(i, i));
}
int saveProcessor = -1;
const int myNode = Ippl::myNode();
const int numNodes = Ippl::getNodes();
const bool scalable = Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE]);
RealDiracMatrix::vector_t p1(4), p2(4);
for (size_t i = 0; i < numberOfParticles; i++) {
for (int j = 0; j < 4; j++) {
p1(j) = gsl_ran_gaussian(randGen_m, 1.0) * variances(j);
p2(j) = gsl_ran_gaussian(randGen_m, 1.0) * variances(4 + j);
}
p1 = boost::numeric::ublas::prod(R1, p1);
p2 = boost::numeric::ublas::prod(R2, p2);
// Save to each processor in turn.
saveProcessor++;
if (saveProcessor >= numNodes)
saveProcessor = 0;
if (scalable || myNode == saveProcessor) {
xDist_m.push_back(p1(0));
pxDist_m.push_back(p1(1) * bgam);
yDist_m.push_back(p1(2));
pyDist_m.push_back(p1(3) * bgam);
tOrZDist_m.push_back(p2(2));
pzDist_m.push_back(p2(3) * bgam);
}
}
}
void Distribution::generateLongFlattopT(size_t numberOfParticles) {
double flattopTime = tPulseLengthFWHM_m
......@@ -2980,7 +3043,7 @@ void Distribution::printDistMatchedGauss(Inform &os) const {
os << "* SIGMAPX = " << sigmaP_m[0] << " [Beta Gamma]" << endl;
os << "* SIGMAPY = " << sigmaP_m[1] << " [Beta Gamma]" << endl;
os << "* SIGMAPZ = " << sigmaP_m[2] << " [Beta Gamma]" << endl;
os << "* AVRGPZ = " << avrgpz_m << " [Beta Gamma]" << endl;
// os << "* AVRGPZ = " << avrgpz_m << " [Beta Gamma]" << endl;
os << "* CORRX = " << correlationMatrix_m(1, 0) << endl;
os << "* CORRY = " << correlationMatrix_m(3, 2) << endl;
......@@ -2989,12 +3052,12 @@ void Distribution::printDistMatchedGauss(Inform &os) const {
os << "* R62 = " << correlationMatrix_m(5, 1) << endl;
os << "* R51 = " << correlationMatrix_m(4, 0) << endl;
os << "* R52 = " << correlationMatrix_m(4, 1) << endl;
os << "* CUTOFFX = " << cutoffR_m[0] << " [units of SIGMAX]" << endl;
os << "* CUTOFFY = " << cutoffR_m[1] << " [units of SIGMAY]" << endl;
os << "* CUTOFFLONG = " << cutoffR_m[2] << " [units of SIGMAZ]" << endl;
os << "* CUTOFFPX = " << cutoffP_m[0] << " [units of SIGMAPX]" << endl;
os << "* CUTOFFPY = " << cutoffP_m[1] << " [units of SIGMAPY]" << endl;
os << "* CUTOFFPZ = " << cutoffP_m[2] << " [units of SIGMAPY]" << endl;
// os << "* CUTOFFX = " << cutoffR_m[0] << " [units of SIGMAX]" << endl;
// os << "* CUTOFFY = " << cutoffR_m[1] << " [units of SIGMAY]" << endl;
// os << "* CUTOFFLONG = " << cutoffR_m[2] << " [units of SIGMAZ]" << endl;
// os << "* CUTOFFPX = " << cutoffP_m[0] << " [units of SIGMAPX]" << endl;
// os << "* CUTOFFPY = " << cutoffP_m[1] << " [units of SIGMAPY]" << endl;
// os << "* CUTOFFPZ = " << cutoffP_m[2] << " [units of SIGMAPY]" << endl;
}
void Distribution::printDistGauss(Inform &os) const {
......@@ -3577,9 +3640,9 @@ void Distribution::setSigmaP_m(double massIneV) {
// Check what input units we are using for momentum.
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
sigmaP_m[0] = converteVToBetaGamma(sigmaP_m[0], massIneV);
sigmaP_m[1] = converteVToBetaGamma(sigmaP_m[1], massIneV);
sigmaP_m[2] = converteVToBetaGamma(sigmaP_m[2], massIneV);
sigmaP_m[0] = Util::convertMomentumeVToBetaGamma(sigmaP_m[0], massIneV);
sigmaP_m[1] = Util::convertMomentumeVToBetaGamma(sigmaP_m[1], massIneV);
sigmaP_m[2] = Util::convertMomentumeVToBetaGamma(sigmaP_m[2], massIneV);
}
}
......@@ -3912,14 +3975,14 @@ void Distribution::setupEmissionModel(PartBunchBase<double, 3> *beam) {
void Distribution::setupEmissionModelAstra(PartBunchBase<double, 3> *beam) {
double wThermal = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]));
pTotThermal_m = converteVToBetaGamma(wThermal, beam->getM());
pTotThermal_m = Util::getBetaGamma(wThermal, beam->getM());
pmean_m = Vector_t(0.0, 0.0, 0.5 * pTotThermal_m);
}
void Distribution::setupEmissionModelNone(PartBunchBase<double, 3> *beam) {
double wThermal = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]));
pTotThermal_m = converteVToBetaGamma(wThermal, beam->getM());
pTotThermal_m = Util::getBetaGamma(wThermal, beam->getM());
double avgPz = std::accumulate(pzDist_m.begin(), pzDist_m.end(), 0.0);
size_t numParticles = pzDist_m.size();
reduce(avgPz, avgPz, OpAddAssign());
......@@ -4073,9 +4136,9 @@ void Distribution::shiftDistCoordinates(double massIneV) {
// Check input momentum units.
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
deltaPx = converteVToBetaGamma(deltaPx, massIneV);
deltaPy = converteVToBetaGamma(deltaPy, massIneV);
deltaPz = converteVToBetaGamma(deltaPz, massIneV);
deltaPx = Util::convertMomentumeVToBetaGamma(deltaPx, massIneV);
deltaPy = Util::convertMomentumeVToBetaGamma(deltaPy, massIneV);
deltaPz = Util::convertMomentumeVToBetaGamma(deltaPz, massIneV);
}
size_t endIdx = startIdx + particlesPerDist_m[i];
......@@ -4356,8 +4419,8 @@ void Distribution::adjustPhaseSpace(double massIneV) {
double deltaPy = Attributes::getReal(itsAttr[Attrib::Distribution::OFFSETPY]);
// Check input momentum units.
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
deltaPx = converteVToBetaGamma(deltaPx, massIneV);
deltaPy = converteVToBetaGamma(deltaPy, massIneV);
deltaPx = Util::convertMomentumeVToBetaGamma(deltaPx, massIneV);
deltaPy = Util::convertMomentumeVToBetaGamma(deltaPy, massIneV);
}
double avrg[6];
......
src/Distribution/Distribution.h
View file @ 05397703
......@@ -19,6 +19,8 @@
#include "AppTypes/SymTenzor.h"
#include "Attributes/Attributes.h"
#include "Distribution/SigmaGenerator.h"
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_qrng.h>
#include <gsl/gsl_rng.h>
......@@ -279,7 +281,6 @@ private:
void eraseTOrZDist();
void eraseBGzDist();
// void printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out);
void addDistributions();
void applyEmissionModel(double lowEnergyLimit, double &px, double &py, double &pz, std::vector<double> &additionalRNs);
void applyEmissModelAstra(double &px, double &py, double &pz, std::vector<double> &additionalRNs);
......@@ -291,7 +292,6 @@ private:
void checkIfEmitted();
void checkParticleNumber(size_t &numberOfParticles);
void chooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUnitsT inputMoUnits);
double converteVToBetaGamma(double valueIneV, double massIneV);
size_t getNumberOfParticlesInFile(std::ifstream &inputFile);
class BinomialBehaviorSplitter {
......@@ -337,6 +337,8 @@ private:
void generateFlattopT(size_t numberOfParticles);
void generateFlattopZ(size_t numberOfParticles);
void generateGaussZ(size_t numberOfParticles);
void generateMatchedGauss(const SigmaGenerator::matrix_t&,
size_t numberOfParticles, double massIneV);
void generateLongFlattopT(size_t numberOfParticles);
void generateTransverseGauss(size_t numberOfParticles);
void initializeBeam(PartBunchBase<double, 3> *beam);
......
src/Distribution/Harmonics.h deleted 100644 → 0
View file @ ab56b03a
//
// Class Harmonics
// This class computes the cyclotron map based on harmonics.
// All functions are copied and translated to C++ from the original program inj2_ana.c of Dr. C. Baumgarten.
//
// Copyright (c) 2014 - 2015, Christian Baumgarten, Paul Scherrer Institut, Villigen PSI, Switzerland
// Matthias Frey, ETH Zürich
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef HARMONICS_H
#define HARMONICS_H
#include <cmath>
#include <fstream>
#include <string>
#include <utility>
#include <vector>
#include "Physics/Physics.h"
#include <boost/numeric/ublas/matrix.hpp>
#include "matrix_vector_operation.h"
template<typename Value_type, typename Size_type>
class Harmonics
{
public:
/// Type of variable
typedef Value_type value_type;
/// Type of size
typedef Size_type size_type;
/// Type for specifying matrices
typedef boost::numeric::ublas::matrix<value_type> matrix_type;
/// Defines starting point of computation
enum START { SECTOR_ENTRANCE, SECTOR_CENTER, SECTOR_EXIT, VALLEY_CENTER };
/// Constructor
/*!
* @param wo is the nominal orbital frequency in Hz
* @param Emin is the minimal energy in MeV
* @param Emax is the maximal energy in MeV
* @param nth is the number of angle steps
* @param nr is the number of radial steps
* @param nSector is the number of sectors
* @param E is the kinetic energy [MeV]
* @param E0 is the potential energy [MeV]
*/
Harmonics(value_type, value_type, value_type, size_type, size_type, size_type, value_type, value_type);
/// Compute all maps of the cyclotron using harmonics
/*!
* @param filename1 name of file 1
* @param filename2 name of file 2
* @param startmod (in center of valley, start of sector, ...)
*/
std::vector<matrix_type> computeMap(const std::string, const std::string, const int);
/// Returns the radial and vertical tunes
std::pair<value_type, value_type> getTunes();
/// Returns the radius
value_type getRadius();
/// Returns step size
value_type getPathLength();
private:
/// Stores the nominal orbital frequency in Hz
value_type wo_m;
/// Stores the minimum energy in MeV
value_type Emin_m;
/// Stores the maximum energy in MeV
value_type Emax_m;
/// Stores the number of angle splits
size_type nth_m;
/// Stores the number of radial splits
size_type nr_m;
/// Stores the number of sectors
size_type nSector_m;
/// Stores the tunes (radial, vertical)
std::pair<value_type, value_type> tunes_m;
/// Stores the radius
value_type radius_m;
/// Stores the path length
value_type ds_m;
/// Stores the energy for which we perform the computation
value_type E_m;
/// Stores the potential energy [MeV]
value_type E0_m;
/// Compute some matrix (ask Dr. C. Baumgarten)
matrix_type __Mix6(value_type, value_type, value_type);
/// Compute some matrix (ask Dr. C. Baumgarten)
matrix_type __Md6(value_type, value_type);
/// Compute some matrix (ask Dr. C. Baumgarten)
matrix_type __Mb6(value_type, value_type, value_type);
/// Compute some matrix (ask Dr. C. Baumgarten)
matrix_type __Mb6k(value_type, value_type, value_type, value_type);
};
// -----------------------------------------------------------------------------------------------------------------------
// PUBLIC MEMBER FUNCTIONS
// -----------------------------------------------------------------------------------------------------------------------
template<typename Value_type, typename Size_type>
Harmonics<Value_type, Size_type>::Harmonics(value_type wo, value_type Emin, value_type Emax,
size_type nth, size_type nr, size_type nSector, value_type E, value_type E0)
: wo_m(wo), Emin_m(Emin), Emax_m(Emax), nth_m(nth), nr_m(nr), nSector_m(nSector), E_m(E), E0_m(E0)
{}
template<typename Value_type, typename Size_type>
std::vector<typename Harmonics<Value_type, Size_type>::matrix_type> Harmonics<Value_type, Size_type>::computeMap(const std::string filename1, const std::string filename2, const int startmod) {
// i --> different energies
// j --> different angles
// ---------------
// unsigned int nth_m = 360; //270;
unsigned int N = 4;
value_type two_pi = 2.0 * M_PI;
value_type tpiN = two_pi / value_type(N);
value_type piN = M_PI / value_type(N);
// unsigned int nr_m = 180;
unsigned int cnt = 0;
bool set = false;
value_type gap = 0.05;
value_type K1 = 0.45;
std::vector<matrix_type> Mcyc(nth_m);
value_type beta_m;
std::vector<value_type> gamma(nr_m), E(nr_m), PC(nr_m), nur(nr_m), nuz(nr_m);
std::vector<value_type> R(nr_m), r(nr_m);
std::vector<value_type> Bmag(nr_m), k(nr_m);
std::vector<value_type> alpha(nr_m), beta(nr_m);
//----------------------------------------------------------------------
// read files
std::ifstream infile2(filename2);
std::ifstream infile1(filename1);
if (!infile1.is_open() || !infile2.is_open()) {
std::cerr << "Couldn't open files!" << std::endl;
std::exit(0);
}
unsigned int n = 0;
value_type rN1, cN1, sN1, aN1, phN1, rN2, cN2, sN2, aN2, phN2;
while (infile1 >> rN1 >> cN1 >> sN1 >> aN1 >> phN1 &&
infile2 >> rN2 >> cN2 >> sN2 >> aN2 >> phN2) {
R[n] = rN1 * 0.01;
alpha[n] = 2.0 * std::acos(aN2 / aN1) / value_type(N);
beta[n] = std::atan(sN1 / cN1) / value_type(N);
value_type f4 = aN1 * M_PI * 0.5 / std::sin(0.5 * alpha[n] * value_type(N));
value_type f8 = aN2 * M_PI / std::sin(alpha[n] * value_type(N));
Bmag[n] = 0.5 * (f4 + f8);
n++;
}
infile1.close();
infile2.close();
//----------------------------------------------------------------------
value_type s = std::sin(piN);
value_type c = std::cos(piN);
std::vector<value_type> dalp(nr_m), dbet(nr_m), len2(nr_m), len(nr_m), t1eff(nr_m), t2eff(nr_m);
std::vector<value_type> t1(nr_m), t2(nr_m), t3(nr_m);
dalp[0] = (alpha[1] - alpha[0]) / (R[1] - R[0]);
dalp[n-1] = (alpha[n-1] - alpha[n-2]) / (R[n-1] - R[n-2]);
dbet[0] = (beta[1] - alpha[0]) / (R[1] - R[0]);
dbet[n-1] = (beta[n-1] - beta[n-2]) / (R[n-1] - R[n-2]);
for (unsigned int i = 1; i < n - 1; ++i) {
dalp[i] = R[i] * (alpha[i+1] - alpha[i-1]) / (R[i+1] - R[i-1]);
dbet[i] = R[i] * (beta[i+1] - beta[i-1]) / (R[i+1] - R[i-1]);
}
for (unsigned int i = 0; i < n; ++i) {
value_type cot = 1.0 / std::tan(0.5 * alpha[i]);
value_type eps1 = std::atan(dbet[i] - 0.5 * dalp[i]);
value_type eps2 = std::atan(dbet[i] + 0.5 * dalp[i]);
value_type phi = piN - 0.5 * alpha[i];
// bending radius
r[i] = R[i] * sin(0.5 * alpha[i]) / s;
len2[i] = R[i] * std::sin(phi);
len[i] = two_pi * r[i] + 2.0 * len2[i] * value_type(N);
value_type g1 = phi + eps1;
value_type g2 = - phi + eps2;
t1[i] = std::tan(g1);
t2[i] = std::tan(g2);
t3[i] = std::tan(phi);
value_type fac = gap / r[i] * K1;
value_type psi = fac * (1.0 + std::sin(g1) * std::sin(g1)) / std::cos(g1);
t1eff[i] = std::tan(g1 - psi);
psi = fac * (1.0 + std::sin(g2) * std::sin(g2)) / std::cos(g2);
t2eff[i] = std::tan(g2 + psi);
beta_m = wo_m / Physics::c * len[i] / two_pi;
gamma[i] = 1.0 / std::sqrt(1.0 - beta_m * beta_m);
E[i] = E0_m * (gamma[i] - 1.0);
PC[i] = E0_m * gamma[i] * beta_m;
Bmag[i] = E0_m * 1.0e6 / Physics::c * beta_m * gamma[i] / r[i] * 10.0;
if (!set && E[i] >= E_m) {
set = true;
cnt = i;
}
}
// (average) field gradient
for (unsigned int i = 0; i < n; ++i) {
value_type dBdr;
if (i == 0)
dBdr = (Bmag[1] - Bmag[0]) / (R[1] - R[0]);
else if (i == n - 1)
dBdr = (Bmag[n-1] - Bmag[n-2]) / (R[n-1] - R[n-2]);
else
dBdr = (Bmag[i+1] - Bmag[i-1]) / (R[i+1] - R[i-1]);
value_type bb = std::sin(piN - 0.5 * alpha[i]) / std::sin(0.5 * alpha[i]);
value_type eps = 2.0 * bb / (1.0 + bb * bb);
value_type BaV = Bmag[i];
k[i] = r[i] * r[i] / (R[i] * BaV) * dBdr * (1.0 + bb * value_type(N) / M_PI * s);
}
for (unsigned int i = 0; i < n; ++i) {
if ((k[1] > -0.999) && (E[i] > Emin_m) && (E[i] < Emax_m + 1.0)) {
value_type kx = std::sqrt(1.0 + k[i]);
value_type Cx = std::cos(tpiN * kx);
value_type Sx = std::sin(tpiN * kx);
value_type pm, ky, Cy, Sy;
if (k[i] >= 0.0) {
pm = 1.0;
ky = std::sqrt(std::abs(k[i]));
Cy = std::cosh(tpiN * ky);
Sy = std::sinh(tpiN * ky);
} else {
pm = -1.0;
ky = std::sqrt(std::abs(k[i]));
Cy = std::cos(tpiN * ky);
Sy = std::sin(tpiN * ky);
}
value_type tav = 0.5 * (t1[i] - t2[i]);
value_type tmul = t1[i] * t2[i];
value_type lrho = 2.0 * len2[i] / r[i];
nur[i] = std::acos(Cx * (1.0 + lrho * tav) + Sx / kx * (tav - lrho * 0.5 * (kx * kx + tmul))) / tpiN;
tav = 0.5 * (t1eff[i] - t2eff[i]);
tmul = t1eff[i] * t2eff[i];
nuz[i] = std::acos(Cy * (1.0 - lrho * tav) + Sy / ky * (0.5 * lrho * (pm * ky * ky - tmul) - tav)) / tpiN;
// std::cout << E[i] << " " << R[i] << " " << nur[i] << " " << nuz[i] << std::endl;
}
}
// -------------------------------------------------------------------------
int i = cnt;
tunes_m = { nur[i], nuz[i] };
radius_m = R[i];
// for (int i = 9; i < 10; ++i) { // n = 1 --> only for one energy
value_type smax, slim, ss, gam2;
smax = len[i] / value_type(N);
slim = r[i] * tpiN;
ds_m = smax / value_type(nth_m);
ss = 0.0;
gam2 = gamma[i] * gamma[i];
matrix_type Mi = __Mix6( t1[i], t1eff[i], r[i]);
matrix_type Mx = __Mix6(- t2[i], - t2eff[i], r[i]);
matrix_type Md = __Md6(ds_m, gam2);
matrix_type Mb = __Mb6k(ds_m / r[i], r[i], k[i], gam2);
unsigned int j = 0;
matrix_type M0 = boost::numeric::ublas::zero_matrix<value_type>(6);
matrix_type M1 = boost::numeric::ublas::zero_matrix<value_type>(6);
switch (startmod) {
case 0: // SECTOR_ENTRANCE
ss = slim - ds_m;
Mcyc[j++] = boost::numeric::ublas::prod(Mb, Mi);
while (ss > ds_m) {
Mcyc[j++] = Mb;
ss -= ds_m;
}
M0 = __Mb6k(ss / r[i], r[i], k[i], gam2);
M1 = __Md6(ds_m - ss, gam2);
Mcyc[j++] = matt_boost::gemmm<matrix_type>(M1, Mx, M0);
ss = smax - slim - (ds_m - ss);
while ((ss > ds_m) && (j < nth_m)) {
Mcyc[j++] = Md;
ss -= ds_m;
}
if (j < nth_m)
Mcyc[j++] = __Md6(ss, gam2);
break;
case 1: // SECTOR_CENTER
j = 0;
ss = slim * 0.5;
while (ss > ds_m) {
Mcyc[j++] = Mb;
ss -= ds_m;
}
M0 = __Mb6k(ss / r[i], r[i], k[i], gam2);
M1 = __Md6(ds_m - ss, gam2);
Mcyc[j++] = matt_boost::gemmm<matrix_type>(M1, Mx, M0);
ss = smax - slim - (ds_m - ss);
while (ss > ds_m) {
Mcyc[j++] = Md;
ss -= ds_m;
}
M0 = __Md6(ss, gam2);
M1 = __Mb6k((ds_m - ss) / r[i], r[i], k[i], gam2);
Mcyc[j++] = matt_boost::gemmm<matrix_type>(M1, Mi, M0);
ss = slim * 0.5 - (ds_m - ss);
while ((ss > ds_m) && (j < nth_m)) {
Mcyc[j++] = Mb;
ss -= ds_m;
}
if (j < nth_m)
Mcyc[j++] = __Mb6k(ss / r[i], r[i], k[i], gam2);
break;
case 2: // SECTOR_EXIT
j = 0;
Mcyc[j++] = boost::numeric::ublas::prod(Md,Mx);
ss = smax - slim - ds_m;
while (ss > ds_m) {
Mcyc[j++] = Md;
ss -= ds_m;
}
M0 = __Md6(ss, gam2);
M1 = __Mb6k((ds_m - ss) / r[i], r[i], k[i], gam2);
Mcyc[j++] = matt_boost::gemmm<matrix_type>(M1, Mi, M0);
ss = slim - (ds_m - ss);
while ((ss > ds_m) && (j < nth_m)) {
Mcyc[j++] = Mb;
ss -= ds_m;
}
if (j < nth_m)
Mcyc[j++] = __Mb6k(ss / r[i], r[i], k[i], gam2);
break;
case 3: // VALLEY_CENTER
default:
j = 0;
ss = (smax - slim) * 0.5;
while (ss > ds_m) {
Mcyc[j++] = Md;
ss -= ds_m;
}
M0 = __Md6(ss, gam2);
M1 = __Mb6k((ds_m - ss) / r[i], r[i], k[i], gam2);
Mcyc[j++] = matt_boost::gemmm<matrix_type>(M1, Mi, M0);
ss = slim - (ds_m - ss);
while (ss > ds_m) {
Mcyc[j++] = Mb;
ss -= ds_m;
}
M0 = __Mb6k(ss / r[i], r[i], k[i], gam2);
M1 = __Md6(ds_m - ss, gam2);
Mcyc[j++] = matt_boost::gemmm<matrix_type>(M1, Mx, M0);
ss = (smax - slim) * 0.5 - (ds_m - ss);
while ((ss > ds_m) && (j < nth_m)) {
Mcyc[j++] = Md;
ss -= ds_m;
}
if (j < nth_m)
Mcyc[j++] = __Md6(ss, gam2);
break;
} // end of switch
// } // end of for
return std::vector<matrix_type>(Mcyc);
}
template<typename Value_type, typename Size_type>
inline std::pair<Value_type, Value_type> Harmonics<Value_type, Size_type>::getTunes() {
return tunes_m;
}
template<typename Value_type, typename Size_type>
inline typename Harmonics<Value_type, Size_type>::value_type Harmonics<Value_type, Size_type>::getRadius() {
return radius_m;
}
template<typename Value_type, typename Size_type>
inline typename Harmonics<Value_type, Size_type>::value_type Harmonics<Value_type, Size_type>::getPathLength() {
return ds_m;
}
// -----------------------------------------------------------------------------------------------------------------------
// PRIVATE MEMBER FUNCTIONS
// -----------------------------------------------------------------------------------------------------------------------
template<typename Value_type, typename Size_type>
typename Harmonics<Value_type, Size_type>::matrix_type Harmonics<Value_type, Size_type>::__Mix6(value_type tx, value_type ty, value_type r) {
matrix_type M = boost::numeric::ublas::identity_matrix<value_type>(6);
M(1,0) = tx / r;
M(3,2) = - ty / r;
return M;
}
template<typename Value_type, typename Size_type>
typename Harmonics<Value_type, Size_type>::matrix_type Harmonics<Value_type, Size_type>::__Md6(value_type l, value_type gam2) {
matrix_type M = boost::numeric::ublas::identity_matrix<value_type>(6);
M(0,1) = l;
M(2,3) = l;
M(4,5) = l / gam2;
return M;
}
template<typename Value_type, typename Size_type>
typename Harmonics<Value_type, Size_type>::matrix_type Harmonics<Value_type, Size_type>::__Mb6(value_type phi, value_type r, value_type gam2) {
matrix_type M = boost::numeric::ublas::identity_matrix<value_type>(6);
value_type C = std::cos(phi);
value_type S = std::sin(phi);
value_type l =r * phi;
M(0,0) = M(1,1) = C;
M(0,1) = S * r;
M(1,0) = - S / r;
M(2,3) = l;
M(4,5) = l / gam2 - l + r * S;
M(4,0) = - (M(1,5) = S);
M(4,1) = - (M(0,5) = r * (1.0 - C));
return M;
}
template<typename Value_type, typename Size_type>
typename Harmonics<Value_type, Size_type>::matrix_type Harmonics<Value_type, Size_type>::__Mb6k(value_type phi, value_type r, value_type k, value_type gam2) {
matrix_type M = boost::numeric::ublas::identity_matrix<value_type>(6);
value_type C,S;
if (k == 0.0) return __Mb6(phi,r,gam2);
value_type fx = std::sqrt(1.0 + k);
value_type fy = std::sqrt(std::abs(k));
value_type c = std::cos(phi * fx);
value_type s = std::sin(phi * fx);
value_type l = r * phi;
if (k > 0.0) {
C = std::cosh(phi * fy);
S = std::sinh(phi * fy);
} else {
C = std::cos(phi * fy);
S = std::sin(phi * fy);
}
M(0,0) = M(1,1) = c;
M(0,1) = s * r / fx;
M(1,0) = - s * fx / r;
M(2,2) = M(3,3) = C;
M(2,3) = S * r / fy;
value_type sign = (std::signbit(k)) ? value_type(-1) : value_type(1);
M(3,2) = sign * S * fy / r;
M(4,5) = l / gam2 - r / (1.0 + k) * (phi - s / fx);
M(4,0)= - (M(1,5) = s / fx);
M(4,1)= - (M(0,5) = r * (1.0 - c) / (1.0 + k));
return M;
}
#endif
\ No newline at end of file
src/Distribution/RealDiracMatrix.cpp 0 → 100644
View file @ 05397703
//
// Class: RealDiracMatrix
// Real Dirac matrix class.
// They're ordered after the paper of Dr. C. Baumgarten: "Use of real Dirac matrices in two-dimensional coupled linear optics".
// The diagonalizing method is based on the paper "Geometrical method of decoupling" (2012) of Dr. C. Baumgarten.
//
// Copyright (c) 2014, 2020 Christian Baumgarten, Paul Scherrer Institut, Villigen PSI, Switzerland
// Matthias Frey, ETH Zürich
// All rights reserved
//
// Implemented as part of the Semester thesis by Matthias Frey
// "Matched Distributions in Cyclotrons"
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "RealDiracMatrix.h"
#include "Utilities/OpalException.h"
#include <cmath>
#include <string>
#include "matrix_vector_operation.h"
RealDiracMatrix::RealDiracMatrix() {};
typename RealDiracMatrix::sparse_matrix_t
RealDiracMatrix::getRDM(short i) {
sparse_matrix_t rdm(4,4,4); // #nrows, #ncols, #non-zeros
switch(i) {
case 0: rdm(0,1) = rdm(2,3) = 1; rdm(1,0) = rdm(3,2) = -1; break;
case 1: rdm(0,1) = rdm(1,0) = -1; rdm(2,3) = rdm(3,2) = 1; break;
case 2: rdm(0,3) = rdm(1,2) = rdm(2,1) = rdm(3,0) = 1; break;
case 3: rdm(0,0) = rdm(2,2) = -1; rdm(1,1) = rdm(3,3) = 1; break;
case 4: rdm(0,0) = rdm(3,3) = -1; rdm(1,1) = rdm(2,2) = 1; break;
case 5: rdm(0,2) = rdm(2,0) = 1; rdm(1,3) = rdm(3,1) = -1; break;
case 6: rdm(0,1) = rdm(1,0) = rdm(2,3) = rdm(3,2) = 1; break;
case 7: rdm(0,3) = rdm(2,1) = 1; rdm(1,2) = rdm(3,0) = -1; break;
case 8: rdm(0,1) = rdm(3,2) = 1; rdm(1,0) = rdm(2,3) = -1; break;
case 9: rdm(0,2) = rdm(1,3) = -1; rdm(2,0) = rdm(3,1) = 1; break;
case 10: rdm(0,2) = rdm(3,1) = 1; rdm(1,3) = rdm(2,0) = -1; break;
case 11: rdm(0,2) = rdm(1,3) = rdm(2,0) = rdm(3,1) = -1; break;
case 12: rdm(0,0) = rdm(1,1) = -1; rdm(2,2) = rdm(3,3) = 1; break;
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: throw OpalException("RealDiracMatrix::getRDM()",
"Index (i = " + std::to_string(i)
+ " out of range: 0 <= i <= 15"); break;
}
return rdm;
}
void RealDiracMatrix::diagonalize(matrix_t& Ms, sparse_matrix_t& R, sparse_matrix_t& invR) {
// R and invR store the total transformation
R = boost::numeric::ublas::identity_matrix<double>(4);
invR = R;
vector_t P(3), E(3), B(3), b;
double mr, mg, mb, eps;
// Lambda function to compute vectors E, P, B and scalar eps (it takes the current Ms as reference argument (--> [&])
auto mult = [&](short i) {
/*
* For computing E, P, B, eps according to formula (C4) from paper:
* Geometrical method of decoupling
*/
matrix_t tmp = boost::numeric::ublas::prod(Ms,getRDM(i))+boost::numeric::ublas::prod(getRDM(i),Ms);
return 0.125*matt_boost::trace(tmp);
};
// 1. Transformation with \gamma_{0}
P = vector_t({ mult(1), mult(2), mult(3)});
E = vector_t({ mult(4), mult(5), mult(6)});
B = vector_t({-mult(7), -mult(8), -mult(9)});
mr = boost::numeric::ublas::inner_prod(E, B); // formula (31), paper: Geometrical method of decoupling
mg = boost::numeric::ublas::inner_prod(B, P); // formula (31), paper: Geometrical method of decoupling
transform(Ms, 0, 0.5 * std::atan2(mg,mr), R, invR);
// 2. Transformation with \gamma_{7}
eps = -mult(0);
P = vector_t({ mult(1), mult(2), mult(3)});
E = vector_t({ mult(4), mult(5), mult(6)});
B = vector_t({-mult(7), -mult(8), -mult(9)});
b = eps * B + matt_boost::cross_prod(E, P); // formula (32), paper: Geometrical method of decoupling
transform(Ms, 7, 0.5 * std::atan2(b(2), b(1)), R, invR);
// 3. Transformation with \gamma_{9}
eps = -mult(0);
P = vector_t({ mult(1), mult(2), mult(3)});
E = vector_t({ mult(4), mult(5), mult(6)});
B = vector_t({-mult(7), -mult(8), -mult(9)});
b = eps * B + matt_boost::cross_prod(E, P);
transform(Ms, 9, -0.5 * std::atan2(b(0), b(1)), R, invR);
// 4. Transformation with \gamma_{2}
eps = -mult(0);
P = vector_t({ mult(1), mult(2), mult(3)});
E = vector_t({ mult(4), mult(5), mult(6)});
B = vector_t({-mult(7), -mult(8), -mult(9)});
mr = boost::numeric::ublas::inner_prod(E, B);
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, 2, 0.5 * std::atanh(mr / b(1)), R, invR);
} else {
transform(Ms, 2, 0.5 * std::atanh(b(1) / mr), R, invR);
}
eps = -mult(0);
P = vector_t({ mult(1), mult(2), mult(3)});
E = vector_t({ mult(4), mult(5), mult(6)});
B = vector_t({-mult(7), -mult(8), -mult(9)});
// formula (31), paper: Geometrical method of decoupling
mr = boost::numeric::ublas::inner_prod(E, B);
mg = boost::numeric::ublas::inner_prod(B, P);
mb = boost::numeric::ublas::inner_prod(E, P);
double P2 = boost::numeric::ublas::inner_prod(P, P);
double E2 = boost::numeric::ublas::inner_prod(E, E);
// 5. Transform with \gamma_{0}
transform(Ms, 0, 0.25 * std::atan2(mb,0.5 * (E2 - P2)), R, invR);
// 6. Transformation with \gamma_{8}
P(0) = mult(1); P(2) = mult(3);
transform(Ms, 8, -0.5 * std::atan2(P(2), P(0)), R, invR);
}
RealDiracMatrix::sparse_matrix_t
RealDiracMatrix::diagonalize(matrix_t& sigma) {
matrix_t S = boost::numeric::ublas::prod(sigma, getRDM(0));
sparse_matrix_t R = boost::numeric::ublas::identity_matrix<double>(4);
sparse_matrix_t invR = boost::numeric::ublas::identity_matrix<double>(4);
sparse_matrix_t iRtot = boost::numeric::ublas::identity_matrix<double>(4);
for (int i = 0; i < 6; ++i) {
update(sigma, i, R, invR);
iRtot = boost::numeric::ublas::prod(iRtot, invR);
S = matt_boost::gemmm<matrix_t>(R, S, invR);
sigma = - boost::numeric::ublas::prod(S, getRDM(0));
}
return iRtot;
}
typename RealDiracMatrix::matrix_t
RealDiracMatrix::symplex(const matrix_t& M) {
/*
* formula(16), p. 3
*/
sparse_matrix_t rdm0 = getRDM(0);
return 0.5 * (M + matt_boost::gemmm<matrix_t>(rdm0,boost::numeric::ublas::trans(M),rdm0));
}
void RealDiracMatrix::transform(matrix_t& M, short i, double phi,
sparse_matrix_t& Rtot, sparse_matrix_t& invRtot)
{
if (phi) { // if phi == 0 --> nothing happens, since R and invR would be identity_matrix matrix
sparse_matrix_t R(4,4), invR(4,4);
transform(i, phi, R, invR);
// update matrices
M = matt_boost::gemmm<matrix_t>(R,M,invR);
Rtot = boost::numeric::ublas::prod(R,Rtot);
invRtot = boost::numeric::ublas::prod(invRtot,invR);
}
}
void RealDiracMatrix::transform(short i, double phi,
sparse_matrix_t& R, sparse_matrix_t& invR)
{
if (phi) { // if phi == 0 --> nothing happens, since R and invR would be identity_matrix matrix
sparse_matrix_t I = boost::numeric::ublas::identity_matrix<double>(4);
if ((i < 7 && i != 0) || (i > 10 && i != 14)) {
R = I * std::cosh(phi) + getRDM(i) * std::sinh(phi);
invR = I * std::cosh(phi) - getRDM(i) * std::sinh(phi);
} else {
R = I * std::cos(phi) + getRDM(i) * std::sin(phi);
invR = I * std::cos(phi) - getRDM(i) * std::sin(phi);
}
}
}
void RealDiracMatrix::update(matrix_t& sigma, short i, sparse_matrix_t& R,
sparse_matrix_t& invR)
{
double s0 = 0.25 * ( sigma(0, 0) + sigma(1, 1) + sigma(2, 2) + sigma(3, 3) );
double s1 = 0.25 * (-sigma(0, 0) + sigma(1, 1) + sigma(2, 2) - sigma(3, 3) );
double s2 = 0.5 * ( sigma(0, 2) - sigma(1, 3) );
double s3 = 0.5 * ( sigma(0, 1) + sigma(2, 3) );
double s4 = 0.5 * ( sigma(0, 1) - sigma(2, 3) );
double s5 = -0.5 * ( sigma(0, 3) + sigma(1, 2) );
double s6 = 0.25 * ( sigma(0, 0) - sigma(1, 1) + sigma(2, 2) - sigma(3, 3) );
double s7 = 0.5 * ( sigma(0, 2) + sigma(1, 3) );
double s8 = 0.25 * ( sigma(0, 0) + sigma(1, 1) - sigma(2, 2) - sigma(3, 3) );
double s9 = 0.5 * ( sigma(0, 3) - sigma(1, 2) );
vector_t P(3); P(0) = s1; P(1) = s2; P(2) = s3;
vector_t E(3); E(0) = s4; E(1) = s5; E(2) = s6;
vector_t B(3); B(0) = s7; B(1) = s8; B(2) = s9;
double mr = boost::numeric::ublas::inner_prod(E, B);
double mg = boost::numeric::ublas::inner_prod(B, P);
double mb = boost::numeric::ublas::inner_prod(E, P);
vector_t b = -s0 * B + matt_boost::cross_prod(E, P);
switch (i) {
case 0:
{
double eps = std::atan2(mg, mr);
transform(0, 0.5 * eps, R, invR);
break;
}
case 1:
{
double eps = std::atan2(b(2), b(1));
transform(7, 0.5 * eps, R, invR);
break;
}
case 2:
{
double eps = - std::atan2(b(0), b(1));
transform(9, 0.5 * eps, R, invR);
break;
}
case 3:
{
double eps = - std::atanh(mr / b(1));
transform(2, 0.5 * eps, R, invR);
break;
}
case 4:
{
double eps = 0.5 * std::atan2(2.0 * mb,
boost::numeric::ublas::inner_prod(E, E) -
boost::numeric::ublas::inner_prod(P, P));
transform(0, 0.5 * eps, R, invR);
break;
}
case 5:
{
double eps = - std::atan2(P(2), P(0));
transform(8, 0.5 * eps, R, invR);
break;
}
default:
{
throw OpalException("RealDiracMatrix::update()",
"Case " + std::to_string(i) +
" not available.");
}
}
}
\ No newline at end of file
src/Distribution/RealDiracMatrix.h 0 → 100644
View file @ 05397703
//
// Class: RealDiracMatrix
// Real Dirac matrix class.
// They're ordered after the paper of Dr. C. Baumgarten: "Use of real Dirac matrices in two-dimensional coupled linear optics".
// The diagonalizing method is based on the paper "Geometrical method of decoupling" (2012) of Dr. C. Baumgarten.
//
// Copyright (c) 2014, 2020 Christian Baumgarten, Paul Scherrer Institut, Villigen PSI, Switzerland
// Matthias Frey, ETH Zürich
// All rights reserved
//
// Implemented as part of the Semester thesis by Matthias Frey
// "Matched Distributions in Cyclotrons"
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef RDM_H
#define RDM_H
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
class RealDiracMatrix
{
public:
/// Sparse matrix type definition
typedef boost::numeric::ublas::compressed_matrix<
double,
boost::numeric::ublas::row_major
> sparse_matrix_t;
/// Dense matrix type definition
typedef boost::numeric::ublas::matrix<double> matrix_t;
/// Dense vector type definition
typedef boost::numeric::ublas::vector<
double,
std::vector<double>
> vector_t;
/// Default constructor (sets only NumOfRDMs and DimOfRDMs)
RealDiracMatrix();
/*!
* @param i specifying the matrix (has to be in the range from 0 to 15)
* @returns the i-th Real Dirac matrix
*/
sparse_matrix_t getRDM(short);
/*!
* Brings a 4x4 symplex matrix into Hamilton form and
* computes the transformation matrix and its inverse
*
* @param Ms is a 4x4 symplex matrix
* @param R is the 4x4 transformation matrix (gets computed)
* @param invR is the 4x4 inverse transformation matrix (gets computed)
*/
void diagonalize(matrix_t&, sparse_matrix_t&, sparse_matrix_t&);
/*!
* Diagonalizes (in-place) the 4x4 sigma matrix. This algorithm
* is Implemented according to https://arxiv.org/abs/1205.3601
*
* @param sigma is the 4x4 sigma matrix
* @returns the inverse of the total transformation
*/
sparse_matrix_t diagonalize(matrix_t&);
/*!
* @param M 4x4 real-valued matrix
* @returns the symplex part of a 4x4 real-valued matrix
*/
matrix_t symplex(const matrix_t&);
private:
/*!
* Applies a rotation to the matrix M by a given angle
*
* @param M is the matrix to be transformed
* @param i is the i-th RDM used for transformation
* @param phi is the angle of rotation
* @param Rtot is a reference to the current transformation matrix
* @param invRtot is a reference to the inverse of the current transformation matrix
*/
void transform(matrix_t&, short, double, sparse_matrix_t&, sparse_matrix_t&);
/*!
* Obtain transformation matrices.
*
* @param i is the i-th RDM used for transformation
* @param phi is the angle of rotation
* @param R is a reference to the transformation matrix
* @param invR is a reference to the inverse of the transformation matrix
*/
void transform(short, double, sparse_matrix_t&, sparse_matrix_t&);
/*!
* Update quantites to decouple the sigma-matrix
*
* @param sigma current state of sigma-matrix
* @param i is the i-th step
* @param R is a reference to the transformation matrix
* @param invR is a reference to the inverse of the transformation matrix
*/
void update(matrix_t& sigma, short i, sparse_matrix_t& R, sparse_matrix_t& invR);
};
#endif
\ No newline at end of file
src/Distribution/SigmaGenerator.cpp 0 → 100644
View file @ 05397703
This diff is collapsed.
src/Distribution/SigmaGenerator.h
View file @ 05397703
This diff is collapsed.
src/Distribution/matrix_vector_operation.h
View file @ 05397703
......@@ -55,10 +55,10 @@ namespace matt_boost {
}
/// Computes the cross product \f$ v_{1}\times v_{2}\f$ of two vectors in \f$ \mathbb{R}^{3} \f$
template<class V>
template<class V, class A>
BOOST_UBLAS_INLINE
ublas::vector<V> cross_prod(ublas::vector<V>& v1, ublas::vector<V>& v2) {
ublas::vector<V> v(v1.size());
ublas::vector<V, A> cross_prod(ublas::vector<V, A>& v1, ublas::vector<V, A>& v2) {
ublas::vector<V, A> 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);
......
src/Elements/OpalElement.cpp
View file @ 05397703
......@@ -568,7 +568,8 @@ void OpalElement::update() {
rotation = rotTheta * (rotPhi * rotPsi);
} else {
if (itsAttr[ORIENTATION]) {
throw OpalException("Line::parse","Parameter orientation is array of 3 values (theta, phi, psi);\n" +
throw OpalException("OpalElement::update",
"Parameter orientation is array of 3 values (theta, phi, psi);\n" +
std::to_string(dir.size()) + " values provided");
}
}
......@@ -577,7 +578,8 @@ void OpalElement::update() {
origin = Vector_t(ori[0], ori[1], ori[2]);
} else {
if (itsAttr[ORIGIN]) {
throw OpalException("Line::parse","Parameter origin is array of 3 values (x, y, z);\n" +
throw OpalException("OpalElement::update",
"Parameter origin is array of 3 values (x, y, z);\n" +
std::to_string(ori.size()) + " values provided");
}
}
......
src/Elements/OpalRBend.cpp
View file @ 05397703
......@@ -24,7 +24,6 @@
#include "Structure/OpalWake.h"
#include "Structure/ParticleMatterInteraction.h"
#include "Utilities/OpalException.h"
#include <cmath>
OpalRBend::OpalRBend():
OpalBend("RBEND",
......@@ -47,10 +46,8 @@ OpalRBend::OpalRBend(const std::string &name, OpalRBend *parent):
OpalRBend::~OpalRBend() {
if(owk_m)
delete owk_m;
if(parmatint_m)
delete parmatint_m;
delete owk_m;
delete parmatint_m;
}
......@@ -131,7 +128,7 @@ void OpalRBend::update() {
double k0s = itsAttr[K0S] ? Attributes::getReal(itsAttr[K0S]) : 0.0;
//JMJ 4/10/2000: above line replaced
// length ? angle / length : angle;
// to avoid closed orbit created by RBEND with defalt K0.
// to avoid closed orbit created by RBEND with default K0.
field.setNormalComponent(1, factor * k0);
field.setSkewComponent(1, factor * Attributes::getReal(itsAttr[K0S]));
field.setNormalComponent(2, factor * Attributes::getReal(itsAttr[K1]));
......@@ -177,8 +174,11 @@ void OpalRBend::update() {
}
// Energy in eV.
if(itsAttr[DESIGNENERGY]) {
if(itsAttr[DESIGNENERGY] && Attributes::getReal(itsAttr[DESIGNENERGY]) != 0.0) {
bend->setDesignEnergy(Attributes::getReal(itsAttr[DESIGNENERGY]), false);
} else if (bend->getName() != "RBEND") {
throw OpalException("OpalRBend::update",
"RBend requires non-zero DESIGNENERGY");
}
bend->setFullGap(Attributes::getReal(itsAttr[GAP]));
......
src/Elements/OpalRBend3D.cpp
View file @ 05397703
......@@ -72,10 +72,8 @@ OpalRBend3D::OpalRBend3D(const std::string &name, OpalRBend3D *parent):
}
OpalRBend3D::~OpalRBend3D() {
if(owk_m)
delete owk_m;
if(parmatint_m)
delete parmatint_m;
delete owk_m;
delete parmatint_m;
}
OpalRBend3D *OpalRBend3D::clone(const std::string &name) {
......@@ -130,8 +128,11 @@ void OpalRBend3D::update() {
bend->setEntranceAngle(e1);
// Energy in eV.
if(itsAttr[DESIGNENERGY]) {
if(itsAttr[DESIGNENERGY] && Attributes::getReal(itsAttr[DESIGNENERGY]) != 0.0) {
bend->setDesignEnergy(Attributes::getReal(itsAttr[DESIGNENERGY]), false);
} else if (bend->getName() != "RBEND3D") {
throw OpalException("OpalRBend3D::update",
"RBend3D requires non-zero DESIGNENERGY");
}
bend->setFullGap(Attributes::getReal(itsAttr[GAP]));
......