36#ifndef VIGRA_WIGNER_MATRIX_HXX
37#define VIGRA_WIGNER_MATRIX_HXX
43#include "mathutil.hxx"
44#include "array_vector.hxx"
46#include "tinyvector.hxx"
47#include "quaternion.hxx"
48#include "clebsch-gordan.hxx"
68 typedef std::complex<Real> Complex;
107 std::string message = std::string(
"WignerMatrix::get_D(): index out of bounds: l=");
108 message <<
l <<
" l_max=" << D.
size() <<
" m=" <<
m <<
" n=" << n <<
"\n";
110 vigra_precondition(
l < D.
size() &&
m+
l <= 2*
l+1 &&
113 return D[
l](n+
l,
m+
l);
117 return Complex(Real(1.0));
127 std::string message = std::string(
"WignerMatrix::get_D(): index out of bounds: l=");
128 message <<
l <<
" l_max=" << l_max <<
"\n";
130 vigra_precondition(
l > 0 &&
l <= l_max, message.c_str());
160 ArrayVector<double> CGcoeff;
161 ArrayVector<Matrix<Complex> > D;
170 for (
int l = 2;
l <=
band+1;
l++)
172 for(
int m = -
l;
m <=
l ;
m++)
174 for(
int n = -
l; n <=
l ; n++)
176 for (
int m2 = -1;
m2 <= 1;
m2++)
178 for (
int n2 = -1;
n2 <= 1;
n2++)
184 CGcoeff.push_back((clebschGordan(
l-1,
m1,1,
m2,
l,
m))*(clebschGordan(
l-1,
n1,1,
n2,
l,n)));
197 double s = std::sin(
theta);
198 double c = std::cos(
theta);
206 if (D.size() < (std::size_t)(
l+1))
210 D[1](0,0) =
emipsi * Complex(Real(0.5*(1.0+c))) *
emiphi;
211 D[1](0,1) = Complex(Real(-s/M_SQRT2)) *
emiphi;
212 D[1](0,2) =
eipsi * Complex(Real(0.5*(1.0-c))) *
emiphi;
213 D[1](1,0) =
emipsi * Complex(Real(s/M_SQRT2));
214 D[1](1,1) = Complex(Real(c));
215 D[1](1,2) =
eipsi * Complex(Real(-s/M_SQRT2));
216 D[1](2,0) =
emipsi * Complex(Real(0.5*(1.0-c))) *
eiphi;
217 D[1](2,1) = Complex(Real(s/M_SQRT2)) *
eiphi;
218 D[1](2,2) =
eipsi * Complex(Real(0.5*(1.0+c))) *
eiphi;
231 if (D.size() < (std::size_t)(
l+1) )
243 D[1](0,0) = Real(0.5);
244 D[1](0,1) = Real(0.5*M_SQRT2);
245 D[1](0,2) = Real(0.5);
246 D[1](1,0) = Real(-0.5*M_SQRT2);
247 D[1](1,1) = Real(0.0);
248 D[1](1,2) = Real(0.5*M_SQRT2);
249 D[1](2,0) = Real(0.5);
250 D[1](2,1) = Real(-0.5*M_SQRT2);
251 D[1](2,2) = Real(0.5);
264 D[
l].init(Real(0.0));
266 for(
int m = -
l;
m <=
l ;
m++)
268 for(
int n = -
l; n <=
l ; n++)
270 for (
int m2 = -1;
m2 <= 1;
m2++)
272 for (
int n2 = -1;
n2 <= 1;
n2++)
278 D[
l](
m+
l,n+
l) += D[1](
m2+1,
n2+1) * D[
l-1](
m1+
l-1,
n1+
l-1) * Real(CGcoeff[cg1cnt++]);
299 for(
int n=1; n<=
band; n++)
305 for(
int m=-
l;
m<=
l;
m++)
308 for (
int h=-
l;
h<=
l;
h++)
Definition array_vector.hxx:514
Class for a single RGB value.
Definition rgbvalue.hxx:128
RGBValue()
Definition rgbvalue.hxx:209
size_type size() const
Definition tinyvector.hxx:913
Class for fixed size vectors.
Definition tinyvector.hxx:1008
computation of Wigner D matrix + rotation functions in SH,VH and R³
Definition wigner-matrix.hxx:64
Matrix< Complex > const & get_D(int l) const
Return the rotation matrix D for the lth band.
Definition wigner-matrix.hxx:125
WignerMatrix(int l_max)
constructor
Definition wigner-matrix.hxx:165
void rotatePH(NestedArray const &PH, Real phi, Real theta, Real psi, NestedArray &PHresult)
Rotate in PH.
Definition wigner-matrix.hxx:291
Complex get_D(int l, int n, int m) const
Get the (n,m) entry of D.
Definition wigner-matrix.hxx:103
void compute_D(int band)
Compute D with fixed theta = pi/2, phi=0, psi=0.
Definition wigner-matrix.hxx:229
FFTWComplex< R > conj(const FFTWComplex< R > &a)
complex conjugate
Definition fftw3.hxx:1030