diff --git a/src/Algorithms/Matrix.h b/src/Algorithms/Matrix.h
index 9f9892fe33168c145e9eebd68f977df4c9c6f521..9b8fee37b45e45099ed3ecef2daa1003c59c7d27 100644
--- a/src/Algorithms/Matrix.h
+++ b/src/Algorithms/Matrix.h
@@ -1,3 +1,20 @@
+//
+// This is header is made to facilitate boost's matrix and vector_t operations in OPAL.
+//
+// Copyright (c) 2023 - 2024, Mohsen Sadr, Paul Scherrer Institute, Villigen PSI, Switzerland
+// 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 OPAL_MATRIX_HH
#define OPAL_MATRIX_HH
diff --git a/src/Algorithms/ParallelCyclotronTracker.cpp b/src/Algorithms/ParallelCyclotronTracker.cpp
index 5eb350f4b15ab53057e728c8c71a72782e472d13..ba28c603e36f5e0b784d05c4ff9acd552354229a 100644
--- a/src/Algorithms/ParallelCyclotronTracker.cpp
+++ b/src/Algorithms/ParallelCyclotronTracker.cpp
@@ -1667,7 +1667,7 @@ void ParallelCyclotronTracker::globalToLocal(ParticleAttrib<Vector_t>& particleV
IpplTimings::startTimer(TransformTimer_m);
particleVectors -= translationToGlobal;
- boost::numeric::ublas::matrix<double> rotation(3,3);
+ matrix_t rotation(3,3);
rotation(0,0) = std::cos(phi);
rotation(0,1) = std::sin(phi);
rotation(0,2) = 0;
@@ -1678,7 +1678,6 @@ void ParallelCyclotronTracker::globalToLocal(ParticleAttrib<Vector_t>& particleV
rotation(2,1) = 0;
rotation(2,2) = 1;
- boost::numeric::ublas::vector<double> particleVector(3);
for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
particleVectors[i] = prod_boost_vector(rotation, particleVectors[i]);
}
@@ -1689,7 +1688,7 @@ void ParallelCyclotronTracker::localToGlobal(ParticleAttrib<Vector_t>& particleV
double phi, Vector_t const translationToGlobal) {
IpplTimings::startTimer(TransformTimer_m);
- boost::numeric::ublas::matrix<double> rotation(3,3);
+ matrix_t rotation(3,3);
rotation(0,0) = std::cos(phi);
rotation(0,1) = -std::sin(phi);
rotation(0,2) = 0;
@@ -1700,7 +1699,6 @@ void ParallelCyclotronTracker::localToGlobal(ParticleAttrib<Vector_t>& particleV
rotation(2,1) = 0;
rotation(2,2) = 1;
- boost::numeric::ublas::vector<double> particleVector(3);
for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
particleVectors[i] = prod_boost_vector(rotation, particleVectors[i]);
}
@@ -1864,7 +1862,7 @@ inline void ParallelCyclotronTracker::normalizeVector(Vector_t& vector) {
inline void ParallelCyclotronTracker::rotateAroundZ(ParticleAttrib<Vector_t>& particleVectors, double const phi) {
// Clockwise rotation of particles 'particleVectors' by 'phi' around Z axis
- boost::numeric::ublas::matrix<double> rotation(3,3);
+ matrix_t rotation(3,3);
rotation(0,0) = std::cos(phi);
rotation(0,1) = std::sin(phi);
rotation(0,2) = 0;
@@ -1875,7 +1873,6 @@ inline void ParallelCyclotronTracker::rotateAroundZ(ParticleAttrib<Vector_t>& pa
rotation(2,1) = 0;
rotation(2,2) = 1;
- boost::numeric::ublas::vector<double> particleVector(3);
for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
particleVectors[i] = prod_boost_vector(rotation, particleVectors[i]);
}
@@ -1884,7 +1881,7 @@ inline void ParallelCyclotronTracker::rotateAroundZ(ParticleAttrib<Vector_t>& pa
inline void ParallelCyclotronTracker::rotateAroundZ(Vector_t& myVector, double const phi) {
// Clockwise rotation of single vector 'myVector' by 'phi' around Z axis
- boost::numeric::ublas::matrix<double> rotation(3,3);
+ matrix_t rotation(3,3);
rotation(0,0) = std::cos(phi);
rotation(0,1) = std::sin(phi);
rotation(0,2) = 0;
@@ -1901,7 +1898,7 @@ inline void ParallelCyclotronTracker::rotateAroundZ(Vector_t& myVector, double c
inline void ParallelCyclotronTracker::rotateAroundX(ParticleAttrib<Vector_t>& particleVectors, double const psi) {
// Clockwise rotation of particles 'particleVectors' by 'psi' around X axis
- boost::numeric::ublas::matrix<double> rotation(3,3);
+ matrix_t rotation(3,3);
rotation(0,0) = 1;
rotation(0,1) = 0;
rotation(0,2) = 0;
@@ -1912,7 +1909,6 @@ inline void ParallelCyclotronTracker::rotateAroundX(ParticleAttrib<Vector_t>& pa
rotation(2,1) = -std::sin(psi);
rotation(2,2) = std::cos(psi);
- boost::numeric::ublas::vector<double> particleVector(3);
for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
particleVectors[i] = prod_boost_vector(rotation, particleVectors[i]);
}
@@ -1921,7 +1917,7 @@ inline void ParallelCyclotronTracker::rotateAroundX(ParticleAttrib<Vector_t>& pa
inline void ParallelCyclotronTracker::rotateAroundX(Vector_t& myVector, double const psi) {
// Clockwise rotation of single vector 'myVector' by 'psi' around X axis
- boost::numeric::ublas::matrix<double> rotation(3,3);
+ matrix_t rotation(3,3);
rotation(0,0) = 1;
rotation(0,1) = 0;
rotation(0,2) = 0;
diff --git a/src/Classic/Algorithms/CoordinateSystemTrafo.h b/src/Classic/Algorithms/CoordinateSystemTrafo.h
index c9a8acc651c9f75185ef2441bfba4ceab521fd93..99e30b1a0d257abc3a9c24307e52569f29ad72be 100644
--- a/src/Classic/Algorithms/CoordinateSystemTrafo.h
+++ b/src/Classic/Algorithms/CoordinateSystemTrafo.h
@@ -4,7 +4,6 @@
#include "Algorithms/Vektor.h"
#include "Algorithms/Quaternion.h"
#include "Algorithms/Matrix.h"
-#include <boost/numeric/ublas/matrix.hpp>
class CoordinateSystemTrafo {
public:
@@ -76,13 +75,7 @@ CoordinateSystemTrafo CoordinateSystemTrafo::inverted() const {
}
inline Vector_t CoordinateSystemTrafo::transformTo(const Vector_t& r) const {
- boost::numeric::ublas::vector<double> result(3);
- for (size_t i = 0; i < 3; ++i) {
- result(i) = r[i]-origin_m[i];
- }
- result = boost::numeric::ublas::prod(rotationMatrix_m, result);
- Vector_t transformedVector(result(0), result(1), result(2));
- return transformedVector;
+ return prod_boost_vector(rotationMatrix_m, r-origin_m);
}
inline Vector_t CoordinateSystemTrafo::transformFrom(const Vector_t& r) const {
@@ -90,24 +83,11 @@ inline Vector_t CoordinateSystemTrafo::transformFrom(const Vector_t& r) const {
}
inline Vector_t CoordinateSystemTrafo::rotateTo(const Vector_t& r) const {
- boost::numeric::ublas::vector<double> result(3);
- for (size_t i = 0; i < 3; ++i) {
- result(i) = r[i];
- }
- result = boost::numeric::ublas::prod(rotationMatrix_m, result);
- Vector_t rotatedVector(result(0), result(1), result(2));
-
- return rotatedVector;
+ return prod_boost_vector(rotationMatrix_m, r);
}
inline Vector_t CoordinateSystemTrafo::rotateFrom(const Vector_t& r) const {
- boost::numeric::ublas::vector<double> result(3);
- for (size_t i = 0; i < 3; ++i) {
- result(i) = r[i];
- }
- result = boost::numeric::ublas::prod(boost::numeric::ublas::trans(rotationMatrix_m), result);
- Vector_t rotatedVector(result(0), result(1), result(2));
- return rotatedVector;
+ return prod_boost_vector(boost::numeric::ublas::trans(rotationMatrix_m), r);
}
#endif
diff --git a/src/Classic/Algorithms/Quaternion.cpp b/src/Classic/Algorithms/Quaternion.cpp
index 828b108d85a556d6593d01c0c3a590faa37d7785..768c752ae9a5a7c93316b7c0853b1152e580e279 100644
--- a/src/Classic/Algorithms/Quaternion.cpp
+++ b/src/Classic/Algorithms/Quaternion.cpp
@@ -2,6 +2,7 @@
#include "Physics/Physics.h"
#include "Utility/RandomNumberGen.h"
#include "Utilities/GeneralClassicException.h"
+#include "Algorithms/Matrix.h"
namespace {
Vector_t normalize(const Vector_t & vec)
@@ -131,7 +132,7 @@ Vector_t Quaternion::rotate(const Vector_t & vec) const
return ((*this) * (quat * (*this).conjugate())).imag();
}
-boost::numeric::ublas::matrix<double> Quaternion::getRotationMatrix() const
+matrix_t Quaternion::getRotationMatrix() const
{
Quaternion rot(*this);
rot.normalize();
diff --git a/src/Classic/Algorithms/Quaternion.h b/src/Classic/Algorithms/Quaternion.h
index bad10a49ff34106b009d39f42a65a2dee798926d..171625377198ce8ed1f4f8fc4d05c3f18ce096d7 100644
--- a/src/Classic/Algorithms/Quaternion.h
+++ b/src/Classic/Algorithms/Quaternion.h
@@ -36,7 +36,7 @@ public:
Vector_t rotate(const Vector_t &) const;
- boost::numeric::ublas::matrix<double> getRotationMatrix() const;
+ matrix_t getRotationMatrix() const;
};
typedef Quaternion Quaternion_t;
diff --git a/src/Distribution/Distribution.h b/src/Distribution/Distribution.h
index 489229f2aee2d79d4c218b0ec41505e4f8cdec5f..a6c90201c15081263d10ef3c34039a6444a89bc9 100644
--- a/src/Distribution/Distribution.h
+++ b/src/Distribution/Distribution.h
@@ -23,7 +23,7 @@
#include "Algorithms/Vektor.h"
#include "Attributes/Attributes.h"
#include "Distribution/SigmaGenerator.h"
-#include <boost/numeric/ublas/matrix.hpp>
+#include "Algorithms/Matrix.h"
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_qrng.h>
@@ -467,7 +467,7 @@ private:
Vector_t cutoffR_m;
Vector_t cutoffP_m;
Vector_t mBinomial_m;
- boost::numeric::ublas::matrix<double> correlationMatrix_m;
+ matrix_t correlationMatrix_m;
// Parameters multiGauss distribution.
double sepPeaks_m;