Commit 36087904 authored by kraus's avatar kraus
Browse files

- adding small class Quaternion (RegTest RingCyclotron runs smoothly)

- updating .gitignore
parent 5b274724
......@@ -154,6 +154,7 @@ classic/5.0/src/Algorithms/TrackIntegrator.cpp -text
classic/5.0/src/Algorithms/TrackIntegrator.h -text
classic/5.0/src/Algorithms/Tracker.cpp -text
classic/5.0/src/Algorithms/Tracker.h -text
classic/5.0/src/Algorithms/Vektor.cpp -text
classic/5.0/src/Algorithms/Vektor.h -text
classic/5.0/src/Algorithms/rbendmap.h -text
classic/5.0/src/BeamlineCore/BeamBeamRep.cpp -text
......
......@@ -20,6 +20,7 @@ set (_SRCS
ThinTracker.cpp
TrackIntegrator.cpp
Tracker.cpp
Vektor.cpp
)
include_directories (
......
#include "Algorithms/Vektor.h"
void normalize(Vector_t & vec)
{
double length = sqrt(dot(vec, vec));
PAssert(length > 1e-12);
vec /= length;
}
Quaternion getQuaternion(Vector_t u, Vector_t v)
{
const double tol1 = 1e-5;
const double tol2 = 1e-8;
Vector_t imaginaryPart(0.0);
normalize(u);
normalize(v);
double k = 1.0;
double cosTheta = dot(u, v);
if (std::abs(cosTheta + 1.0) < tol1) {
Vector_t xAxis(1.0, 0, 0);
Vector_t zAxis(0, 0, 1.0);
if (std::abs(1.0 - dot(u, xAxis)) > tol2)
imaginaryPart = cross(u, xAxis);
else
imaginaryPart = cross(u, zAxis);
k = 0.0;
cosTheta = 0.0;
} else {
imaginaryPart = cross(u, v);
}
Quaternion result(cosTheta + k, imaginaryPart);
result.normalize();
return result;
}
Quaternion Quaternion::operator*(const double & d) const
{
Quaternion result(d * this->real(), d * this->imag());
return result;
}
Quaternion Quaternion::operator*(const Quaternion & other) const
{
Quaternion result(*this);
return result *= other;
}
Quaternion & Quaternion::operator*=(const Quaternion & other)
{
Vector_t imagThis = this->imag();
Vector_t imagOther = other.imag();
*this = Quaternion((*this)(0) * other(0) - dot(imagThis, imagOther),
(*this)(0) * imagOther + other(0) * imagThis + cross(imagThis, imagOther));
return *this;
}
Quaternion Quaternion::operator/(const double & d) const
{
Quaternion result(this->real() / d, this->imag() / d);
return result;
}
Quaternion & Quaternion::normalize()
{
PAssert(this->Norm() > 1e-12);
(*this) /= this->length();
return (*this);
}
Quaternion Quaternion::inverse() const
{
PAssert(this->Norm() > 1e-12);
Quaternion returnValue = conjugate();
return returnValue.normalize();
}
Vector_t Quaternion::rotate(const Vector_t & vec) const
{
PAssert(this->isUnit());
Quaternion quat(vec);
return ((*this) * (quat * (*this).conjugate())).imag();
}
Tenzor<double, 3> Quaternion::getRotationMatrix() const
{
Tenzor<double, 3> mat(1 - 2 * ((*this)(2)*(*this)(2) + (*this)(3)*(*this)(3)),
2 * ((*this)(0) * (*this)(3) + (*this)(1) * (*this)(2)),
2 * (-(*this)(0) * (*this)(2) + (*this)(1) * (*this)(3)),
2 * (-(*this)(0) * (*this)(3) + (*this)(1) * (*this)(2)),
1 - 2 * ((*this)(1)*(*this)(1) + (*this)(3)*(*this)(3)),
2 * ((*this)(0) * (*this)(1) + (*this)(2) * (*this)(3)),
2 * ((*this)(0) * (*this)(2) + (*this)(1) * (*this)(3)),
2 * (-(*this)(0) * (*this)(1) + (*this)(2) * (*this)(3)),
1 - 2 * ((*this)(1)*(*this)(1) + (*this)(2)*(*this)(2)));
return mat;
}
#ifndef OPAL_VEKTOR_HH
#define OPAL_VEKTOR_HH
#include <limits>
#include "AppTypes/Vektor.h"
#include "AppTypes/Tenzor.h"
typedef Vektor<double, 3> Vector_t;
typedef Vektor<double, 4> Quaternion_t; // used in ParallelCyclotron.{cpp,h}
void normalize(Vector_t & vec);
class Quaternion: public Vektor<double, 4> {
public:
Quaternion();
Quaternion(const Quaternion &);
Quaternion(const double &, const double &, const double &, const double &);
Quaternion(const Vector_t &);
Quaternion(const double &, const Vector_t &);
Quaternion operator*(const double &) const;
Quaternion operator*(const Quaternion &) const;
Quaternion& operator*=(const Quaternion &);
Quaternion operator/(const double &) const;
double Norm() const;
double length() const;
Quaternion & normalize();
bool isUnit() const;
bool isPure() const;
bool isPureUnit() const;
Quaternion inverse() const;
Quaternion conjugate() const;
double real() const;
Vector_t imag() const;
Vector_t rotate(const Vector_t &) const;
Tenzor<double, 3> getRotationMatrix() const;
};
typedef Quaternion Quaternion_t;
//typedef Vektor<double, 4> Quaternion_t; // used in ParallelCyclotron.{cpp,h}
Quaternion getQuaternion(Vector_t u, Vector_t v);
inline
Quaternion::Quaternion():
Vektor<double, 4>(0.0)
{}
inline
Quaternion::Quaternion(const Quaternion & quat):
Vektor<double, 4>(quat)
{}
inline
Quaternion::Quaternion(const double & x0, const double & x1, const double & x2, const double & x3):
Vektor<double, 4>(x0, x1, x2, x3)
{}
inline
Quaternion::Quaternion(const Vector_t & vec):
Quaternion(0.0, vec(0), vec(1), vec(2))
{}
inline
Quaternion::Quaternion(const double & realPart, const Vector_t & vec):
Quaternion(realPart, vec(0), vec(1), vec(2))
{}
inline
double Quaternion::Norm() const
{
return dot(*this, *this);
}
inline
double Quaternion::length() const
{
return sqrt(this->Norm());
}
inline
bool Quaternion::isUnit() const
{
return (std::abs(this->Norm() - 1.0) < 1e-12);
}
inline
bool Quaternion::isPure() const
{
return (std::abs((*this)(0)) < 1e-12);
}
inline
bool Quaternion::isPureUnit() const
{
return (this->isPure() && this->isUnit());
}
inline
Quaternion Quaternion::conjugate() const
{
Quaternion quat(this->real(), -this->imag());
return quat;
}
inline
double Quaternion::real() const
{
return (*this)(0);
}
inline
Vector_t Quaternion::imag() const
{
Vector_t vec((*this)(1), (*this)(2), (*this)(3));
return vec;
}
#endif
......@@ -1526,6 +1526,19 @@ void ParallelTTracker::Tracker_Default() {
secondaryFlg_m = false;
dtTrack_m = itsBunch->getdT();
double testPhi = 22.0 / 180.0 * Physics::pi;
double testTheta = 112 / 180.0 * Physics::pi;
Vector_t testA(1.0, 0,0);
Vector_t testB(cos(testPhi)*sin(testTheta), sin(testPhi)*sin(testTheta), cos(testTheta));
Quaternion testQ = getQuaternion(testA, testB);
Vector_t testC = testQ.rotate(testA);
msg << __FILE__ << ": " << __LINE__ << "\n"
<< testC << "\n"
<< testB << endl;
// msg << testB << "\n"
// << testC << endl;
// upper limit of particle number when we do field emission and secondary emission
// simulation. Could be reset to another value in input file with MAXPARTSNUM.
maxNparts_m = 100000000;
......
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