#ifndef ATOOLS_Math_Tensor_Contractions_Epsilon_H
#define ATOOLS_Math_Tensor_Contractions_Epsilon_H
/********************************************************************************************
****** ******
****** epsilon tensor contractions explicitely ******
****** ******
********************************************************************************************/
template<typename Scalar>
Lorentz_Ten3<Scalar>
ATOOLS::EpsilonTensorContraction(const Vec4<Scalar>& v) {
// T^{mu,nu,rho} = eps^{mu,nu,rho,alpha} g_{alpha,alpha'} v^{alpha}
return Lorentz_Ten3<Scalar>(
0. , 0. , 0. , 0. ,
0. , 0. ,-v[3], v[2],
0. , v[3], 0. ,-v[1],
0. ,-v[2], v[1], 0. ,
0. , 0. , v[3],-v[2],
0. , 0. , 0. , 0. ,
-v[3], 0. , 0. , v[0],
v[2], 0. ,-v[0], 0. ,
0. ,-v[3], 0. , v[1],
v[3], 0. , 0. ,-v[0],
0. , 0. , 0. , 0. ,
-v[1], v[0], 0. , 0. ,
0. , v[2],-v[1], 0. ,
-v[2], 0. , v[0], 0. ,
v[1],-v[0], 0. , 0. ,
0. , 0. , 0. , 0. );
}
template<typename Scal1, typename Scal2>
Lorentz_Ten2<PROMOTE(Scal1,Scal2)>
ATOOLS::EpsilonTensorContraction(const Vec4<Scal1>& p, const Vec4<Scal2>& q) {
// returns epsilon tensor contractions
// T^{mu,nu} = eps^{mu,nu,alpha,beta}g_{alpha,alpha'}g_{beta,beta'}p^alpha' q^beta'
return Lorentz_Ten2<PROMOTE(Scal1,Scal2)>(
0. , -p[2]*q[3]+p[3]*q[2] , p[1]*q[3]-p[3]*q[1] , -p[1]*q[2]+p[2]*q[1] ,
p[2]*q[3]-p[3]*q[2] , 0. , -p[0]*q[3]+p[3]*q[0] , p[0]*q[2]-p[2]*q[0] ,
-p[1]*q[3]+p[3]*q[1] , p[0]*q[3]-p[3]*q[0] , 0. , -p[0]*q[1]+p[1]*q[0] ,
p[1]*q[2]-p[2]*q[1] , -p[0]*q[2]+p[2]*q[0] , p[0]*q[1]-p[1]*q[0] , 0. );
}
template<typename Scal1, typename Scal2, typename Scal3>
Vec4<PROMOTE(Scal1,PROMOTE(Scal2,Scal3))>
ATOOLS::EpsilonTensorContraction(const Vec4<Scal1>& p1, const Vec4<Scal2>& p2,
const Vec4<Scal3>& p3) {
// v^{mu} = eps^{mu,alpha,beta,gamma}g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
// p1^{alpha'} p2^{beta'} p3^{gamma'} and permutations thereof
// maybe uncomment, but doubtful if any computational improvement
// if (IsEqual(p1,p2) || IsEqual(p1,p3) || IsEqual(p2,p3))
// return Vec4<PROMOTE(Scal1,PROMOTE(Scal2,Scal3))>(0.,0.,0.,0.);
return EpsilonTensorContraction(BuildTensor(p1,p2,p3),2,3,4);
}
template<typename Scal1, typename Scal2, typename Scal3, typename Scal4>
PROMOTE(Scal1,PROMOTE(Scal2,PROMOTE(Scal3,Scal4)))
ATOOLS::EpsilonTensorContraction(const Vec4<Scal1>& p1, const Vec4<Scal2>& p2,
const Vec4<Scal3>& p3, const Vec4<Scal4>& p4) {
// c = eps^{alpha,beta,gamma,delta} p1^alpha p2^beta p3^gamma p4^delta
// gamma_{alpha,alpha'} gamma_{beta,beta'}
// gamma_{gamma,gamma'} gamma_{delta,delta'}
return p1[0]*p2[1]*p3[2]*p4[3]
-p1[0]*p2[1]*p3[3]*p4[2]
-p1[0]*p2[2]*p3[1]*p4[3]
+p1[0]*p2[2]*p3[3]*p4[1]
+p1[0]*p2[3]*p3[1]*p4[2]
-p1[0]*p2[3]*p3[2]*p4[1]
-p1[1]*p2[0]*p3[2]*p4[3]
+p1[1]*p2[0]*p3[3]*p4[2]
+p1[1]*p2[2]*p3[0]*p4[3]
-p1[1]*p2[2]*p3[3]*p4[0]
-p1[1]*p2[3]*p3[0]*p4[2]
+p1[1]*p2[3]*p3[2]*p4[0]
+p1[2]*p2[0]*p3[1]*p4[3]
-p1[2]*p2[0]*p3[3]*p4[1]
-p1[2]*p2[1]*p3[0]*p4[3]
+p1[2]*p2[1]*p3[3]*p4[0]
+p1[2]*p2[3]*p3[0]*p4[1]
-p1[2]*p2[3]*p3[1]*p4[0]
-p1[3]*p2[0]*p3[1]*p4[2]
+p1[3]*p2[0]*p3[2]*p4[1]
+p1[3]*p2[1]*p3[0]*p4[2]
-p1[3]*p2[1]*p3[2]*p4[0]
-p1[3]*p2[2]*p3[0]*p4[1]
+p1[3]*p2[2]*p3[1]*p4[0];
}
template<typename Scalar>
Lorentz_Ten2<Scalar>
ATOOLS::EpsilonTensorContraction(const Lorentz_Ten2<Scalar>& t, int i, int j) {
// returns the contraction eps^{mu,nu,alpha,beta}g_{alpha,alpha'}g_{beta,beta'}t^{alpha'beta'}
if ((i==3) && (j==4))
return Lorentz_Ten2<Scalar>(
0. , -t.at(2,3)+t.at(3,2) , t.at(1,3)-t.at(3,1) , -t.at(1,2)+t.at(2,1) ,
t.at(2,3)-t.at(3,2) , 0. , -t.at(0,3)+t.at(3,0) , t.at(0,2)-t.at(2,0) ,
-t.at(1,3)+t.at(3,1) , t.at(0,3)-t.at(3,0) , 0. , -t.at(0,1)+t.at(1,0) ,
t.at(1,2)-t.at(2,1) , -t.at(0,2)+t.at(2,0) , t.at(0,1)-t.at(1,0) , 0. );
else if ((j==3) && (i==4))
return -Lorentz_Ten2<Scalar>(
0. , -t.at(2,3)+t.at(3,2) , t.at(1,3)-t.at(3,1) , -t.at(1,2)+t.at(2,1) ,
t.at(2,3)-t.at(3,2) , 0. , -t.at(0,3)+t.at(3,0) , t.at(0,2)-t.at(2,0) ,
-t.at(1,3)+t.at(3,1) , t.at(0,3)-t.at(3,0) , 0. , -t.at(0,1)+t.at(1,0) ,
t.at(1,2)-t.at(2,1) , -t.at(0,2)+t.at(2,0) , t.at(0,1)-t.at(1,0) , 0. );
return Lorentz_Ten2<Scalar>();
}
template<typename Scal1, typename Scal2>
Vec4<PROMOTE(Scal1,Scal2)>
ATOOLS::EpsilonTensorContraction(const Lorentz_Ten2<Scal1>& t, int i, int j,
const Vec4<Scal2>& p, int k) {
// v^{mu} = eps^{mu,alpha,beta,gamma}g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
// t^{alpha',beta'} p^{gamma'} and permutations thereof
return EpsilonTensorContraction(BuildTensor(t,p),i,j,k);
}
template<typename Scalar>
Vec4<Scalar>
ATOOLS::EpsilonTensorContraction(const Lorentz_Ten3<Scalar>& ten, int i, int j, int k) {
Lorentz_Ten3<Scalar> t;
// v^{mu} = eps^{mu,alpha,beta,gamma} t^{alpha',beta',gamma'} g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
if ((i==2) && (j==3) && (k==4))
t = ten;
// v^{mu} = eps^{mu,alpha,beta,gamma} t^{alpha',gamma',beta'} g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
else if ((i==2) && (j==4) && (k==3))
t = ten.Transpose(2,3);
// v^{mu} = eps^{mu,alpha,beta,gamma} t^{beta',alpha',gamma'} g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
else if ((i==3) && (j==2) && (k==4))
t = ten.Transpose(1,2);
// v^{mu} = eps^{mu,alpha,beta,gamma} t^{beta',gamma',alpha'} g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
else if ((i==3) && (j==4) && (k==2))
{t = ten.Transpose(2,3); t = t.Transpose(1,2);}
// v^{mu} = eps^{mu,alpha,beta,gamma} t^{gamma',alpha',beta'} g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
else if ((i==4) && (j==2) && (k==3))
{t = ten.Transpose(1,2); t = t.Transpose(2,3);}
// v^{mu} = eps^{mu,alpha,beta,gamma} t^{gamma',beta',alpha'} g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
else if ((i==4) && (j==3) && (k==2))
t = ten.Transpose(1,3);
else
return Vec4<Scalar>(0.,0.,0.,0.);
return Vec4<Scalar>(
t.at(3,2,1)-t.at(2,3,1)+t.at(1,3,2)-t.at(3,1,2)+t.at(2,1,3)-t.at(1,2,3),
t.at(3,2,0)-t.at(2,3,0)+ t.at(0,3,2)-t.at(3,0,2)+t.at(2,0,3)-t.at(0,2,3),
t.at(1,3,0)-t.at(3,1,0)+t.at(3,0,1)-t.at(0,3,1)+ t.at(0,1,3)-t.at(1,0,3),
t.at(2,1,0)-t.at(1,2,0)+t.at(0,2,1)-t.at(2,1,0)+t.at(1,0,2)-t.at(0,1,2) );
}
template<typename Scal1, typename Scal2>
PROMOTE(Scal1,Scal2)
ATOOLS::EpsilonTensorContraction(const Lorentz_Ten2<Scal1>& ten, int i, int j,
const Lorentz_Ten2<Scal2>& sen, int k, int l) {
// c = eps^{alpha,beta,gamma,delta} t^{alpha',beta'} s^{gamma',delta'}
// gamma_{alpha,alpha'} gamma_{beta,beta'}
// gamma_{gamma,gamma'} gamma_{delta,delta'}
return EpsilonTensorContraction(BuildTensor(ten,sen),i,j,k,l);
}
template<typename Scal1, typename Scal2>
PROMOTE(Scal1,Scal2)
ATOOLS::EpsilonTensorContraction(const Lorentz_Ten3<Scal1>& ten, int i, int j, int k,
const Vec4<Scal2>& v, int l) {
// c = eps^{alpha,beta,gamma,delta} t^{alpha',beta',gamma'} v^delta'
// gamma_{alpha,alpha'} gamma_{beta,beta'}
// gamma_{gamma,gamma'} gamma_{delta,delta'}
return EpsilonTensorContraction(BuildTensor(ten,v),i,j,k,l);
}
template<typename Scalar>
Scalar
ATOOLS::EpsilonTensorContraction(const Lorentz_Ten4<Scalar>& ten, int i, int j, int k, int l) {
Lorentz_Ten4<Scalar> t;
// c = eps^{alpha,beta,gamma,delta} t^{alpha',beta',gamma',delta'}
// gamma_{alpha,alpha'} gamma_{beta,beta'}
// gamma_{gamma,gamma'} gamma_{delta,delta'}
if ((i == 1) && (j == 2) && (k == 3) && (l == 4))
{t = ten;}
else if ((i == 1) && (j == 2) && (k == 4) && (l == 3))
{t = ten.Transpose(3,4);}
else if ((i == 1) && (j == 3) && (k == 2) && (l == 4))
{t = ten.Transpose(2,3);}
else if ((i == 1) && (j == 3) && (k == 4) && (l == 2))
{t = ten.Transpose(3,4); t = t.Transpose(2,3);}
else if ((i == 1) && (j == 4) && (k == 2) && (l == 3))
{t = ten.Transpose(2,3); t = t.Transpose(3,4);}
else if ((i == 1) && (j == 4) && (k == 3) && (l == 2))
{t = ten.Transpose(2,4);}
else if ((i == 2) && (j == 1) && (k == 3) && (l == 4))
{t = ten.Transpose(1,2);}
else if ((i == 2) && (j == 1) && (k == 4) && (l == 3))
{t = ten.Transpose(1,2); t = t.Transpose(3,4);}
else if ((i == 2) && (j == 3) && (k == 1) && (l == 4))
{t = ten.Transpose(2,3); t = t.Transpose(3,4);}
else if ((i == 2) && (j == 3) && (k == 4) && (l == 1))
{t = ten.Transpose(3,4); t = t.Transpose(2,3); t = t.Transpose(1,2);}
else if ((i == 2) && (j == 4) && (k == 1) && (l == 3))
{t = ten.Transpose(1,3); t = t.Transpose(2,3); t = t.Transpose(3,4);}
else if ((i == 2) && (j == 4) && (k == 3) && (l == 1))
{t = ten.Transpose(2,4); t = t.Transpose(1,2);}
else if ((i == 3) && (j == 1) && (k == 2) && (l == 4))
{t = ten.Transpose(1,2); t = t.Transpose(2,3);}
else if ((i == 3) && (j == 1) && (k == 4) && (l == 2))
{t = ten.Transpose(1,2); t = t.Transpose(2,4); t = t.Transpose(3,4);}
else if ((i == 3) && (j == 2) && (k == 1) && (l == 4))
{t = ten.Transpose(1,3);}
else if ((i == 3) && (j == 2) && (k == 4) && (l == 1))
{t = ten.Transpose(1,4); t = t.Transpose(3,4);}
else if ((i == 3) && (j == 4) && (k == 1) && (l == 2))
{t = ten.Transpose(1,3); t = t.Transpose(2,4);}
else if ((i == 3) && (j == 4) && (k == 2) && (l == 1))
{t = ten.Transpose(1,4); t = t.Transpose(2,3); t = t.Transpose(3,4);}
else if ((i == 4) && (j == 1) && (k == 2) && (l == 3))
{t = ten.Transpose(1,2); t = t.Transpose(2,3); t = t.Transpose(3,4);}
else if ((i == 4) && (j == 1) && (k == 3) && (l == 2))
{t = ten.Transpose(1,2); t = t.Transpose(2,4);}
else if ((i == 4) && (j == 2) && (k == 1) && (l == 3))
{t = ten.Transpose(1,3); t = t.Transpose(3,4);}
else if ((i == 4) && (j == 2) && (k == 3) && (l == 1))
{t = ten.Transpose(1,4);}
else if ((i == 4) && (j == 3) && (k == 1) && (l == 2))
{t = ten.Transpose(1,3); t = t.Transpose(2,4); t = t.Transpose(3,4);}
else if ((i == 4) && (j == 3) && (k == 2) && (l == 1))
{t = ten.Transpose(1,4); t = t.Transpose(2,3);}
else
return 0.;
return t.at(0,1,2,3)
-t.at(0,1,3,2)
-t.at(0,2,1,3)
+t.at(0,2,3,1)
+t.at(0,3,1,2)
-t.at(0,3,2,1)
-t.at(1,0,2,3)
+t.at(1,0,3,2)
+t.at(1,2,0,3)
-t.at(1,2,3,0)
-t.at(1,3,0,2)
+t.at(1,3,2,0)
+t.at(2,0,1,3)
-t.at(2,0,3,1)
-t.at(2,1,0,3)
+t.at(2,1,3,0)
+t.at(2,3,0,1)
-t.at(2,3,1,0)
-t.at(3,0,1,2)
+t.at(3,0,2,1)
+t.at(3,1,0,2)
-t.at(3,1,2,0)
-t.at(3,2,0,1)
+t.at(3,2,1,0);
}
#endif