#ifndef ATOOLS_Math_Tensor_Contractions_H
#define ATOOLS_Math_Tensor_Contractions_H
/********************************************************************************************
****** ******
****** tensor with itself ******
****** ******
********************************************************************************************/
template<typename Scalar>
Scalar
ATOOLS::Contraction(const Lorentz_Ten2<Scalar>& t) {
// c = t^{alpha,beta} g_{alpha,beta}
return t.at(0,0)-t.at(1,1)-t.at(2,2)-t.at(3,3);
}
template<typename Scalar>
Vec4<Scalar>
ATOOLS::Contraction(const Lorentz_Ten3<Scalar>& ten, int i, int j) {
Lorentz_Ten3<Scalar> t;
// v^{mu} = t{mu,alpha,alpha'} g_{alpha,alpha'}
if (((i==2) && (j==3)) || ((j==2) && (i==3)))
t = ten;
// v^{mu} = t{alpha,mu,alpha'} g_{alpha,alpha'}
else if (((i==1) && (j==3)) || ((j==1) && (i==3)))
t = ten.Transpose(1,2);
// v^{mu} = t{alpha,alpha',mu} g_{alpha,alpha'}
else if (((i==1) && (j==2)) || ((j==1) && (i==2)))
t = ten.Transpose(1,3);
else
return Vec4<Scalar>(0.,0.,0.,0.);
return Vec4<Scalar>(t.at(0,0,0)-t.at(0,1,1)-t.at(0,2,2)-t.at(0,3,3),
t.at(1,0,0)-t.at(1,1,1)-t.at(1,2,2)-t.at(1,3,3),
t.at(2,0,0)-t.at(2,1,1)-t.at(2,2,2)-t.at(2,3,3),
t.at(3,0,0)-t.at(3,1,1)-t.at(3,2,2)-t.at(3,3,3));
}
template<typename Scalar>
Lorentz_Ten2<Scalar>
ATOOLS::Contraction(const Lorentz_Ten4<Scalar>& ten, int i, int j) {
Lorentz_Ten4<Scalar> t;
// T^{mu,nu} = t^{mu,nu,alpha,alpha'} g_{alpha,alpha'}
if (((i==3) && (j==4)) || ((i==4) && (j==3)))
t = ten;
// T^{mu,nu} = t^{mu,alpha,nu,alpha'} g_{alpha,alpha'}
else if (((i==2) && (j==4)) || ((i==4) && (j==2)))
t = ten.Transpose(2,3);
// T^{mu,nu} = t^{alpha,mu,nu,alpha'} g_{alpha,alpha'}
else if (((i==1) && (j==4)) || ((i==4) && (j==1)))
{t = ten.Transpose(1,2); t = t.Transpose(2,3);}
// T^{mu,nu} = t^{mu,alpha,alpha',nu} g_{alpha,alpha'}
else if (((i==2) && (j==3)) || ((i==3) && (j==2)))
t = ten.Transpose(2,4);
// T^{mu,nu} = t^{alpha,mu,alpha',nu} g_{alpha,alpha'}
else if (((i==1) && (j==3)) || ((i==3) && (j==1)))
{t = ten.Transpose(1,2); t = t.Transpose(2,4);}
// T^{mu,nu} = t^{alpha,alpha',mu,nu} g_{alpha,alpha'}
else if (((i==1) && (j==2)) || ((i==2) && (j==1)))
{t = ten.Transpose(1,3); t = t.Transpose(2,4);}
else
return Lorentz_Ten2<Scalar>();
return Lorentz_Ten2<Scalar>(t.at(0,0,0,0)-t.at(0,0,1,1)-t.at(0,0,2,2)-t.at(0,0,3,3),
t.at(1,0,0,0)-t.at(1,0,1,1)-t.at(1,0,2,2)-t.at(1,0,3,3),
t.at(2,0,0,0)-t.at(2,0,1,1)-t.at(2,0,2,2)-t.at(2,0,3,3),
t.at(3,0,0,0)-t.at(3,0,1,1)-t.at(3,0,2,2)-t.at(3,0,3,3),
t.at(0,1,0,0)-t.at(0,1,1,1)-t.at(0,1,2,2)-t.at(0,1,3,3),
t.at(1,1,0,0)-t.at(1,1,1,1)-t.at(1,1,2,2)-t.at(1,1,3,3),
t.at(2,1,0,0)-t.at(2,1,1,1)-t.at(2,1,2,2)-t.at(2,1,3,3),
t.at(3,1,0,0)-t.at(3,1,1,1)-t.at(3,1,2,2)-t.at(3,1,3,3),
t.at(0,2,0,0)-t.at(0,2,1,1)-t.at(0,2,2,2)-t.at(0,2,3,3),
t.at(1,2,0,0)-t.at(1,2,1,1)-t.at(1,2,2,2)-t.at(1,2,3,3),
t.at(2,2,0,0)-t.at(2,2,1,1)-t.at(2,2,2,2)-t.at(2,2,3,3),
t.at(3,2,0,0)-t.at(3,2,1,1)-t.at(3,2,2,2)-t.at(3,2,3,3),
t.at(0,3,0,0)-t.at(0,3,1,1)-t.at(0,3,2,2)-t.at(0,3,3,3),
t.at(1,3,0,0)-t.at(1,3,1,1)-t.at(1,3,2,2)-t.at(1,3,3,3),
t.at(2,3,0,0)-t.at(2,3,1,1)-t.at(2,3,2,2)-t.at(2,3,3,3),
t.at(3,3,0,0)-t.at(3,3,1,1)-t.at(3,3,2,2)-t.at(3,3,3,3));
}
template<typename Scalar>
Scalar
ATOOLS::Contraction(const Lorentz_Ten4<Scalar>& ten, int i, int j, int k, int l) {
Lorentz_Ten4<Scalar> t;
// c = t^{alpha,alpha',beta,beta'} g_{alpha,alpha'} g_{beta,beta'}
if (((((i==1) && (j==2)) || ((i==2) && (j==1))) &&
(((k==3) && (l==4)) || ((k==4) && (l==3)))) ||
((((i==3) && (j==4)) || ((i==4) && (j==3))) &&
(((k==1) && (l==2)) || ((k==2) && (l==1)))))
t = ten;
// c = t^{alpha,beta,alpha',beta'} g_{alpha,alpha'} g_{beta,beta'}
else if (((((i==1) && (j==3)) || ((i==3) && (j==1))) &&
(((k==2) && (l==4)) || ((k==4) && (l==2)))) ||
((((i==2) && (j==4)) || ((i==4) && (j==2))) &&
(((k==1) && (l==3)) || ((k==3) && (l==1)))))
t = ten.Transpose(2,3);
// c = t^{alpha,beta,beta',alpha'} g_{alpha,alpha'} g_{beta,beta'}
else if (((((i==1) && (j==4)) || ((i==4) && (j==1))) &&
(((k==2) && (l==3)) || ((k==3) && (l==2)))) ||
((((i==2) && (j==3)) || ((i==3) && (j==2))) &&
(((k==1) && (l==4)) || ((k==4) && (l==1)))))
t = ten.Transpose(2,4);
else
return Scalar(0.);
return Scalar( t.at(0,0,0,0)-t.at(0,0,1,1)-t.at(0,0,2,2)-t.at(0,0,3,3)
-t.at(1,1,0,0)+t.at(1,1,1,1)+t.at(1,1,2,2)+t.at(1,1,3,3)
-t.at(2,2,0,0)+t.at(2,2,1,1)+t.at(2,2,2,2)+t.at(2,2,3,3)
-t.at(3,3,0,0)+t.at(3,3,1,1)+t.at(3,3,2,2)+t.at(3,3,3,3));
}
/********************************************************************************************
****** ******
****** tensor with vector ******
****** ******
********************************************************************************************/
template<typename Scal1, typename Scal2>
Vec4<PROMOTE(Scal1,Scal2)>
ATOOLS::Contraction(const Lorentz_Ten2<Scal1>& ten, int i, const Vec4<Scal2>& p) {
Lorentz_Ten2<Scal1> t;
if ((i==1) || (i==2)) {
// v^{mu} = t^{alpha,mu} p^{beta} g_{alpha,beta}
if (i==1)
t = ten.Transpose();
// v^{mu} = t^{mu,alpha} p^{beta} g_{alpha,beta}
else if (i==2)
t = ten;
return Vec4<PROMOTE(Scal1,Scal2)>(
t.at(0,0)*p[0]-t.at(0,1)*p[1]-t.at(0,2)*p[2]-t.at(0,3)*p[3],
t.at(1,0)*p[0]-t.at(1,1)*p[1]-t.at(1,2)*p[2]-t.at(1,3)*p[3],
t.at(2,0)*p[0]-t.at(2,1)*p[1]-t.at(2,2)*p[2]-t.at(2,3)*p[3],
t.at(3,0)*p[0]-t.at(3,1)*p[1]-t.at(3,2)*p[2]-t.at(3,3)*p[3]);
}
return Vec4<PROMOTE(Scal1,Scal2)>();
}
template<typename Scal1, typename Scal2, typename Scal3>
PROMOTE(PROMOTE(Scal1,Scal2),Scal3)
ATOOLS::Contraction(const Lorentz_Ten2<Scal1>& t,
const Vec4<Scal2>& p, const Vec4<Scal3>& q) {
// c = t^{alpha,beta} p^{alpha'} q^{beta'} g_{alpha',alpha} g_{beta,beta'}
return q[0]*(t.at(0,0)*p[0]-t.at(1,0)*p[1]-t.at(2,0)*p[2]-t.at(3,0)*p[3])
-q[1]*(t.at(0,1)*p[0]-t.at(1,1)*p[1]-t.at(2,1)*p[2]-t.at(3,1)*p[3])
-q[2]*(t.at(0,2)*p[0]-t.at(1,2)*p[1]-t.at(2,2)*p[2]-t.at(3,2)*p[3])
-q[3]*(t.at(0,3)*p[0]-t.at(1,3)*p[1]-t.at(2,3)*p[2]-t.at(3,3)*p[3]);
}
template<typename Scal1, typename Scal2>
Lorentz_Ten2<PROMOTE(Scal1,Scal2)>
ATOOLS::Contraction(const Lorentz_Ten3<Scal1>& ten, int i, const Vec4<Scal2>& v) {
Lorentz_Ten3<Scal1> t;
// c^{mu,nu} = t^{alpha,mu,nu} g_{alpha,alpha'} v^{alpha'}
if (i==1)
{t = ten.Transpose(1,2); t = t.Transpose(2,3);}
// c^{mu,nu} = t^{mu,alpha,nu} g_{alpha,alpha'} v^{alpha'}
else if (i==2)
t = ten.Transpose(2,3);
// c^{mu,nu} = t^{mu,nu,alpha} g_{alpha,alpha'} v^{alpha'}
else if (i==3)
t = ten;
else
return Lorentz_Ten2<PROMOTE(Scal1,Scal2)>();
return Lorentz_Ten2<PROMOTE(Scal1,Scal2)>(
t.at(0,0,0)*v[0]-t.at(0,0,1)*v[1]-t.at(0,0,2)*v[2]-t.at(0,0,3)*v[3],
t.at(1,0,0)*v[0]-t.at(1,0,1)*v[1]-t.at(1,0,2)*v[2]-t.at(1,0,3)*v[3],
t.at(2,0,0)*v[0]-t.at(2,0,1)*v[1]-t.at(2,0,2)*v[2]-t.at(2,0,3)*v[3],
t.at(3,0,0)*v[0]-t.at(3,0,1)*v[1]-t.at(3,0,2)*v[2]-t.at(3,0,3)*v[3],
t.at(0,1,0)*v[0]-t.at(0,1,1)*v[1]-t.at(0,1,2)*v[2]-t.at(0,1,3)*v[3],
t.at(1,1,0)*v[0]-t.at(1,1,1)*v[1]-t.at(1,1,2)*v[2]-t.at(1,1,3)*v[3],
t.at(2,1,0)*v[0]-t.at(2,1,1)*v[1]-t.at(2,1,2)*v[2]-t.at(2,1,3)*v[3],
t.at(3,1,0)*v[0]-t.at(3,1,1)*v[1]-t.at(3,1,2)*v[2]-t.at(3,1,3)*v[3],
t.at(0,2,0)*v[0]-t.at(0,2,1)*v[1]-t.at(0,2,2)*v[2]-t.at(0,2,3)*v[3],
t.at(1,2,0)*v[0]-t.at(1,2,1)*v[1]-t.at(1,2,2)*v[2]-t.at(1,2,3)*v[3],
t.at(2,2,0)*v[0]-t.at(2,2,1)*v[1]-t.at(2,2,2)*v[2]-t.at(2,2,3)*v[3],
t.at(3,2,0)*v[0]-t.at(3,2,1)*v[1]-t.at(3,2,2)*v[2]-t.at(3,2,3)*v[3],
t.at(0,3,0)*v[0]-t.at(0,3,1)*v[1]-t.at(0,3,2)*v[2]-t.at(0,3,3)*v[3],
t.at(1,3,0)*v[0]-t.at(1,3,1)*v[1]-t.at(1,3,2)*v[2]-t.at(1,3,3)*v[3],
t.at(2,3,0)*v[0]-t.at(2,3,1)*v[1]-t.at(2,3,2)*v[2]-t.at(2,3,3)*v[3],
t.at(3,3,0)*v[0]-t.at(3,3,1)*v[1]-t.at(3,3,2)*v[2]-t.at(3,3,3)*v[3]);
}
template<typename Scal1, typename Scal2>
Lorentz_Ten3<PROMOTE(Scal1,Scal2)>
ATOOLS::Contraction(const Lorentz_Ten4<Scal1>& ten, int i, const Vec4<Scal2>& v) {
Lorentz_Ten4<Scal1> t;
// c^{mu,nu,rho} = t^{alpha,mu,nu,rho} g_{alpha,alpha'} v^{alpha'}
if (i==1)
{t = ten.Transpose(1,2); t = t.Transpose(2,3); t = t.Transpose(3,4);}
// c^{mu,nu,rho} = t^{mu,alpha,nu,rho} g_{alpha,alpha'} v^{alpha'}
else if (i==2)
{t = ten.Transpose(2,3); t = t.Transpose(3,4);}
// c^{mu,nu,rho} = t^{mu,nu,alpha,rho} g_{alpha,alpha'} v^{alpha'}
else if (i==3)
t = ten.Transpose(3,4);
// c^{mu,nu,rho} = t^{mu,nu,rho,alpha} g_{alpha,alpha'} v^{alpha'}
else if (i==4)
t = ten;
else
return Lorentz_Ten3<PROMOTE(Scal1,Scal2)>();
return Lorentz_Ten3<PROMOTE(Scal1,Scal2)>(
t.at(0,0,0,0)*v[0]-t.at(0,0,0,1)*v[1]-t.at(0,0,0,2)*v[2]-t.at(0,0,0,3)*v[3],
t.at(1,0,0,0)*v[0]-t.at(1,0,0,1)*v[1]-t.at(1,0,0,2)*v[2]-t.at(1,0,0,3)*v[3],
t.at(2,0,0,0)*v[0]-t.at(2,0,0,1)*v[1]-t.at(2,0,0,2)*v[2]-t.at(2,0,0,3)*v[3],
t.at(3,0,0,0)*v[0]-t.at(3,0,0,1)*v[1]-t.at(3,0,0,2)*v[2]-t.at(3,0,0,3)*v[3],
t.at(0,1,0,0)*v[0]-t.at(0,1,0,1)*v[1]-t.at(0,1,0,2)*v[2]-t.at(0,1,0,3)*v[3],
t.at(1,1,0,0)*v[0]-t.at(1,1,0,1)*v[1]-t.at(1,1,0,2)*v[2]-t.at(1,1,0,3)*v[3],
t.at(2,1,0,0)*v[0]-t.at(2,1,0,1)*v[1]-t.at(2,1,0,2)*v[2]-t.at(2,1,0,3)*v[3],
t.at(3,1,0,0)*v[0]-t.at(3,1,0,1)*v[1]-t.at(3,1,0,2)*v[2]-t.at(3,1,0,3)*v[3],
t.at(0,2,0,0)*v[0]-t.at(0,2,0,1)*v[1]-t.at(0,2,0,2)*v[2]-t.at(0,2,0,3)*v[3],
t.at(1,2,0,0)*v[0]-t.at(1,2,0,1)*v[1]-t.at(1,2,0,2)*v[2]-t.at(1,2,0,3)*v[3],
t.at(2,2,0,0)*v[0]-t.at(2,2,0,1)*v[1]-t.at(2,2,0,2)*v[2]-t.at(2,2,0,3)*v[3],
t.at(3,2,0,0)*v[0]-t.at(3,2,0,1)*v[1]-t.at(3,2,0,2)*v[2]-t.at(3,2,0,3)*v[3],
t.at(0,3,0,0)*v[0]-t.at(0,3,0,1)*v[1]-t.at(0,3,0,2)*v[2]-t.at(0,3,0,3)*v[3],
t.at(1,3,0,0)*v[0]-t.at(1,3,0,1)*v[1]-t.at(1,3,0,2)*v[2]-t.at(1,3,0,3)*v[3],
t.at(2,3,0,0)*v[0]-t.at(2,3,0,1)*v[1]-t.at(2,3,0,2)*v[2]-t.at(2,3,0,3)*v[3],
t.at(3,3,0,0)*v[0]-t.at(3,3,0,1)*v[1]-t.at(3,3,0,2)*v[2]-t.at(3,3,0,3)*v[3],
t.at(0,0,1,0)*v[0]-t.at(0,0,1,1)*v[1]-t.at(0,0,1,2)*v[2]-t.at(0,0,1,3)*v[3],
t.at(1,0,1,0)*v[0]-t.at(1,0,1,1)*v[1]-t.at(1,0,1,2)*v[2]-t.at(1,0,1,3)*v[3],
t.at(2,0,1,0)*v[0]-t.at(2,0,1,1)*v[1]-t.at(2,0,1,2)*v[2]-t.at(2,0,1,3)*v[3],
t.at(3,0,1,0)*v[0]-t.at(3,0,1,1)*v[1]-t.at(3,0,1,2)*v[2]-t.at(3,0,1,3)*v[3],
t.at(0,1,1,0)*v[0]-t.at(0,1,1,1)*v[1]-t.at(0,1,1,2)*v[2]-t.at(0,1,1,3)*v[3],
t.at(1,1,1,0)*v[0]-t.at(1,1,1,1)*v[1]-t.at(1,1,1,2)*v[2]-t.at(1,1,1,3)*v[3],
t.at(2,1,1,0)*v[0]-t.at(2,1,1,1)*v[1]-t.at(2,1,1,2)*v[2]-t.at(2,1,1,3)*v[3],
t.at(3,1,1,0)*v[0]-t.at(3,1,1,1)*v[1]-t.at(3,1,1,2)*v[2]-t.at(3,1,1,3)*v[3],
t.at(0,2,1,0)*v[0]-t.at(0,2,1,1)*v[1]-t.at(0,2,1,2)*v[2]-t.at(0,2,1,3)*v[3],
t.at(1,2,1,0)*v[0]-t.at(1,2,1,1)*v[1]-t.at(1,2,1,2)*v[2]-t.at(1,2,1,3)*v[3],
t.at(2,2,1,0)*v[0]-t.at(2,2,1,1)*v[1]-t.at(2,2,1,2)*v[2]-t.at(2,2,1,3)*v[3],
t.at(3,2,1,0)*v[0]-t.at(3,2,1,1)*v[1]-t.at(3,2,1,2)*v[2]-t.at(3,2,1,3)*v[3],
t.at(0,3,1,0)*v[0]-t.at(0,3,1,1)*v[1]-t.at(0,3,1,2)*v[2]-t.at(0,3,1,3)*v[3],
t.at(1,3,1,0)*v[0]-t.at(1,3,1,1)*v[1]-t.at(1,3,1,2)*v[2]-t.at(1,3,1,3)*v[3],
t.at(2,3,1,0)*v[0]-t.at(2,3,1,1)*v[1]-t.at(2,3,1,2)*v[2]-t.at(2,3,1,3)*v[3],
t.at(3,3,1,0)*v[0]-t.at(3,3,1,1)*v[1]-t.at(3,3,1,2)*v[2]-t.at(3,3,1,3)*v[3],
t.at(0,0,2,0)*v[0]-t.at(0,0,2,1)*v[1]-t.at(0,0,2,2)*v[2]-t.at(0,0,2,3)*v[3],
t.at(1,0,2,0)*v[0]-t.at(1,0,2,1)*v[1]-t.at(1,0,2,2)*v[2]-t.at(1,0,2,3)*v[3],
t.at(2,0,2,0)*v[0]-t.at(2,0,2,1)*v[1]-t.at(2,0,2,2)*v[2]-t.at(2,0,2,3)*v[3],
t.at(3,0,2,0)*v[0]-t.at(3,0,2,1)*v[1]-t.at(3,0,2,2)*v[2]-t.at(3,0,2,3)*v[3],
t.at(0,1,2,0)*v[0]-t.at(0,1,2,1)*v[1]-t.at(0,1,2,2)*v[2]-t.at(0,1,2,3)*v[3],
t.at(1,1,2,0)*v[0]-t.at(1,1,2,1)*v[1]-t.at(1,1,2,2)*v[2]-t.at(1,1,2,3)*v[3],
t.at(2,1,2,0)*v[0]-t.at(2,1,2,1)*v[1]-t.at(2,1,2,2)*v[2]-t.at(2,1,2,3)*v[3],
t.at(3,1,2,0)*v[0]-t.at(3,1,2,1)*v[1]-t.at(3,1,2,2)*v[2]-t.at(3,1,2,3)*v[3],
t.at(0,2,2,0)*v[0]-t.at(0,2,2,1)*v[1]-t.at(0,2,2,2)*v[2]-t.at(0,2,2,3)*v[3],
t.at(1,2,2,0)*v[0]-t.at(1,2,2,1)*v[1]-t.at(1,2,2,2)*v[2]-t.at(1,2,2,3)*v[3],
t.at(2,2,2,0)*v[0]-t.at(2,2,2,1)*v[1]-t.at(2,2,2,2)*v[2]-t.at(2,2,2,3)*v[3],
t.at(3,2,2,0)*v[0]-t.at(3,2,2,1)*v[1]-t.at(3,2,2,2)*v[2]-t.at(3,2,2,3)*v[3],
t.at(0,3,2,0)*v[0]-t.at(0,3,2,1)*v[1]-t.at(0,3,2,2)*v[2]-t.at(0,3,2,3)*v[3],
t.at(1,3,2,0)*v[0]-t.at(1,3,2,1)*v[1]-t.at(1,3,2,2)*v[2]-t.at(1,3,2,3)*v[3],
t.at(2,3,2,0)*v[0]-t.at(2,3,2,1)*v[1]-t.at(2,3,2,2)*v[2]-t.at(2,3,2,3)*v[3],
t.at(3,3,2,0)*v[0]-t.at(3,3,2,1)*v[1]-t.at(3,3,2,2)*v[2]-t.at(3,3,2,3)*v[3],
t.at(0,0,3,0)*v[0]-t.at(0,0,3,1)*v[1]-t.at(0,0,3,2)*v[2]-t.at(0,0,3,3)*v[3],
t.at(1,0,3,0)*v[0]-t.at(1,0,3,1)*v[1]-t.at(1,0,3,2)*v[2]-t.at(1,0,3,3)*v[3],
t.at(2,0,3,0)*v[0]-t.at(2,0,3,1)*v[1]-t.at(2,0,3,2)*v[2]-t.at(2,0,3,3)*v[3],
t.at(3,0,3,0)*v[0]-t.at(3,0,3,1)*v[1]-t.at(3,0,3,2)*v[2]-t.at(3,0,3,3)*v[3],
t.at(0,1,3,0)*v[0]-t.at(0,1,3,1)*v[1]-t.at(0,1,3,2)*v[2]-t.at(0,1,3,3)*v[3],
t.at(1,1,3,0)*v[0]-t.at(1,1,3,1)*v[1]-t.at(1,1,3,2)*v[2]-t.at(1,1,3,3)*v[3],
t.at(2,1,3,0)*v[0]-t.at(2,1,3,1)*v[1]-t.at(2,1,3,2)*v[2]-t.at(2,1,3,3)*v[3],
t.at(3,1,3,0)*v[0]-t.at(3,1,3,1)*v[1]-t.at(3,1,3,2)*v[2]-t.at(3,1,3,3)*v[3],
t.at(0,2,3,0)*v[0]-t.at(0,2,3,1)*v[1]-t.at(0,2,3,2)*v[2]-t.at(0,2,3,3)*v[3],
t.at(1,2,3,0)*v[0]-t.at(1,2,3,1)*v[1]-t.at(1,2,3,2)*v[2]-t.at(1,2,3,3)*v[3],
t.at(2,2,3,0)*v[0]-t.at(2,2,3,1)*v[1]-t.at(2,2,3,2)*v[2]-t.at(2,2,3,3)*v[3],
t.at(3,2,3,0)*v[0]-t.at(3,2,3,1)*v[1]-t.at(3,2,3,2)*v[2]-t.at(3,2,3,3)*v[3],
t.at(0,3,3,0)*v[0]-t.at(0,3,3,1)*v[1]-t.at(0,3,3,2)*v[2]-t.at(0,3,3,3)*v[3],
t.at(1,3,3,0)*v[0]-t.at(1,3,3,1)*v[1]-t.at(1,3,3,2)*v[2]-t.at(1,3,3,3)*v[3],
t.at(2,3,3,0)*v[0]-t.at(2,3,3,1)*v[1]-t.at(2,3,3,2)*v[2]-t.at(2,3,3,3)*v[3],
t.at(3,3,3,0)*v[0]-t.at(3,3,3,1)*v[1]-t.at(3,3,3,2)*v[2]-t.at(3,3,3,3)*v[3]);
}
/********************************************************************************************
****** ******
****** tensor with tensor ******
****** ******
********************************************************************************************/
template<typename Scal1, typename Scal2>
PROMOTE(Scal1,Scal2)
ATOOLS::Contraction(const Lorentz_Ten2<Scal1>& ten, int i, int j,
const Lorentz_Ten2<Scal2>& sen, int k, int l) {
Lorentz_Ten2<Scal1> t; Lorentz_Ten2<Scal2> s;
// c = t^{alpha,beta} s^{alpha',beta'} g_{alpha',alpha} g_{beta,beta'}
if (((i==1) && (j==2) && (k==1) && (l==2)) ||
((i==2) && (j==1) && (k==2) && (l==1)))
{t = ten.Transpose(); s = sen;}
// c = t^{alpha,beta} s^{beta',alpha'} g_{alpha',alpha} g_{beta,beta'}
else if (((i==1) && (j==2) && (k==2) && (l==1)) ||
((i==2) && (j==1) && (k==1) && (l==2)))
{t = ten; s = sen;}
else
return PROMOTE(Scal1,Scal2)();
return PROMOTE(Scal1,Scal2)
(t.at(0,0)*s.at(0,0)-t.at(1,0)*s.at(0,1)-t.at(2,0)*s.at(0,2)-t.at(3,0)*s.at(0,3))
-(t.at(0,1)*s.at(1,0)-t.at(1,1)*s.at(1,1)-t.at(2,1)*s.at(1,2)-t.at(3,1)*s.at(1,3))
-(t.at(0,2)*s.at(2,0)-t.at(1,2)*s.at(2,1)-t.at(2,2)*s.at(2,2)-t.at(3,2)*s.at(2,3))
-(t.at(0,3)*s.at(3,0)-t.at(1,3)*s.at(3,1)-t.at(2,3)*s.at(3,2)-t.at(3,3)*s.at(3,3));
}
template<typename Scal1, typename Scal2>
Lorentz_Ten2<PROMOTE(Scal1,Scal2)>
ATOOLS::Contraction(const Lorentz_Ten2<Scal1>& ten,int i,
const Lorentz_Ten2<Scal2>& sen, int j) {
Lorentz_Ten2<Scal1> t; Lorentz_Ten2<Scal2> s;
if (((i==1) || (i==2)) && ((j==1) || (j==2))) {
// T^{mu,nu} = t^{alpha,mu} s^{beta,nu} g_{alpha,beta}
if ((i==1) && (j==1))
{t = ten.Transpose(); s = sen;}
// T^{mu,nu} = t^{alpha,mu} s^{nu,beta} g_{alpha,beta}
else if ((i==1) && (j==2))
{t = ten.Transpose(); s = sen.Transpose();}
// T^{mu,nu} = t^{mu,alpha} s^{beta,nu} g_{alpha,beta}
else if ((i==2) && (j==1))
{t = ten; s = sen;}
// T^{mu,nu} = t^{mu,alpha} s^{nu,beta} g_{alpha,beta}
else if ((i==2) && (j==2))
{t = ten; s = sen.Transpose();}
return Lorentz_Ten2<PROMOTE(Scal1,Scal2)>(
s.at(0,0)*t.at(0,0)-s.at(1,0)*t.at(0,1)-s.at(2,0)*t.at(0,2)-s.at(3,0)*t.at(0,3),
s.at(0,0)*t.at(1,0)-s.at(1,0)*t.at(1,1)-s.at(2,0)*t.at(1,2)-s.at(3,0)*t.at(1,3),
s.at(0,0)*t.at(2,0)-s.at(1,0)*t.at(2,1)-s.at(2,0)*t.at(2,2)-s.at(3,0)*t.at(2,3),
s.at(0,0)*t.at(3,0)-s.at(1,0)*t.at(3,1)-s.at(2,0)*t.at(3,2)-s.at(3,0)*t.at(3,3),
s.at(0,1)*t.at(0,0)-s.at(1,1)*t.at(0,1)-s.at(2,1)*t.at(0,2)-s.at(3,1)*t.at(0,3),
s.at(0,1)*t.at(1,0)-s.at(1,1)*t.at(1,1)-s.at(2,1)*t.at(1,2)-s.at(3,1)*t.at(1,3),
s.at(0,1)*t.at(2,0)-s.at(1,1)*t.at(2,1)-s.at(2,1)*t.at(2,2)-s.at(3,1)*t.at(2,3),
s.at(0,1)*t.at(3,0)-s.at(1,1)*t.at(3,1)-s.at(2,1)*t.at(3,2)-s.at(3,1)*t.at(3,3),
s.at(0,2)*t.at(0,0)-s.at(1,2)*t.at(0,1)-s.at(2,2)*t.at(0,2)-s.at(3,2)*t.at(0,3),
s.at(0,2)*t.at(1,0)-s.at(1,2)*t.at(1,1)-s.at(2,2)*t.at(1,2)-s.at(3,2)*t.at(1,3),
s.at(0,2)*t.at(2,0)-s.at(1,2)*t.at(2,1)-s.at(2,2)*t.at(2,2)-s.at(3,2)*t.at(2,3),
s.at(0,2)*t.at(3,0)-s.at(1,2)*t.at(3,1)-s.at(2,2)*t.at(3,2)-s.at(3,2)*t.at(3,3),
s.at(0,3)*t.at(0,0)-s.at(1,3)*t.at(0,1)-s.at(2,3)*t.at(0,2)-s.at(3,3)*t.at(0,3),
s.at(0,3)*t.at(1,0)-s.at(1,3)*t.at(1,1)-s.at(2,3)*t.at(1,2)-s.at(3,3)*t.at(1,3),
s.at(0,3)*t.at(2,0)-s.at(1,3)*t.at(2,1)-s.at(2,3)*t.at(2,2)-s.at(3,3)*t.at(2,3),
s.at(0,3)*t.at(3,0)-s.at(1,3)*t.at(3,1)-s.at(2,3)*t.at(3,2)-s.at(3,3)*t.at(3,3));
}
return Lorentz_Ten2<PROMOTE(Scal1,Scal2)>();
}
template<typename Scal1, typename Scal2>
Vec4<PROMOTE(Scal1,Scal2)>
ATOOLS::Contraction(const Lorentz_Ten3<Scal1>& ten, int i, int j,
const Lorentz_Ten2<Scal2>& sen, int k, int l) {
Lorentz_Ten3<Scal1> t; Lorentz_Ten2<Scal2> s;
// v^{mu} = t^{mu,alpha,beta} s^{alpha',beta'} g_{alpha,alpha'} g_{beta,beta'}
if (((i==2) && (j==3) && (k==1) && (l==2)) || ((i==3) && (j==2) && (k==2) && (l==1)))
// v^{mu} = t^{alpha,mu,beta} s^{alpha',beta'} g_{alpha,alpha'} g_{beta,beta'}
{t = ten; s = sen;}
else if (((i==1) && (j==3) && (k==1) && (l==2)) || ((i==3) && (j==1) && (k==2) && (l==1)))
{t = ten.Transpose(1,2); s = sen;}
// v^{mu} = t^{alpha,beta,mu} s^{alpha',beta'} g_{alpha,alpha'} g_{beta,beta'}
else if (((i==1) && (j==2) && (k==1) && (l==2)) || ((i==2) && (j==1) && (k==2) && (l==1)))
{t = ten.Transpose(2,3); t = t.Transpose(1,2); s = sen;}
// v^{mu} = t^{mu,alpha,beta} s^{beta',alpha'} g_{alpha,alpha'} g_{beta,beta'}
else if (((i==2) && (j==3) && (k==2) && (l==1)) || ((i==3) && (j==2) && (k==1) && (l==2)))
{t = ten; s = sen.Transpose();}
// v^{mu} = t^{alpha,mu,beta} s^{beta',alpha'} g_{alpha,alpha'} g_{beta,beta'}
else if (((i==1) && (j==3) && (k==2) && (l==1)) || ((i==3) && (j==1) && (k==1) && (l==2)))
{t = ten.Transpose(1,2); s = sen.Transpose();}
// v^{mu} = t^{alpha,beta,mu} s^{beta',alpha'} g_{alpha,alpha'} g_{beta,beta'}
else if (((i==1) && (j==2) && (k==2) && (l==1)) || ((i==2) && (j==1) && (k==1) && (l==2)))
{t = ten.Transpose(2,3); t = t.Transpose(1,2); s = sen.Transpose();}
else
return Vec4<PROMOTE(Scal1,Scal2)>(0.,0.,0.,0.);
return Vec4<PROMOTE(Scal1,Scal2)>(
t.at(0,0,0)*s.at(0,0)-t.at(0,1,0)*s.at(1,0)-t.at(0,2,0)*s.at(2,0)-t.at(0,3,0)*s.at(3,0)
-t.at(0,0,1)*s.at(0,1)+t.at(0,1,1)*s.at(1,1)+t.at(0,2,1)*s.at(2,1)+t.at(0,3,1)*s.at(3,1)
-t.at(0,0,2)*s.at(0,2)+t.at(0,1,2)*s.at(1,2)+t.at(0,2,2)*s.at(2,2)+t.at(0,3,2)*s.at(3,2)
-t.at(0,0,3)*s.at(0,3)+t.at(0,1,3)*s.at(1,3)+t.at(0,2,3)*s.at(2,3)+t.at(0,3,3)*s.at(3,3),
t.at(1,0,0)*s.at(0,0)-t.at(1,1,0)*s.at(1,0)-t.at(1,2,0)*s.at(2,0)-t.at(1,3,0)*s.at(3,0)
-t.at(1,0,1)*s.at(0,1)+t.at(1,1,1)*s.at(1,1)+t.at(1,2,1)*s.at(2,1)+t.at(1,3,1)*s.at(3,1)
-t.at(1,0,2)*s.at(0,2)+t.at(1,1,2)*s.at(1,2)+t.at(1,2,2)*s.at(2,2)+t.at(1,3,2)*s.at(3,2)
-t.at(1,0,3)*s.at(0,3)+t.at(1,1,3)*s.at(1,3)+t.at(1,2,3)*s.at(2,3)+t.at(1,3,3)*s.at(3,3),
t.at(2,0,0)*s.at(0,0)-t.at(2,1,0)*s.at(1,0)-t.at(2,2,0)*s.at(2,0)-t.at(2,3,0)*s.at(3,0)
-t.at(2,0,1)*s.at(0,1)+t.at(2,1,1)*s.at(1,1)+t.at(2,2,1)*s.at(2,1)+t.at(2,3,1)*s.at(3,1)
-t.at(2,0,2)*s.at(0,2)+t.at(2,1,2)*s.at(1,2)+t.at(2,2,2)*s.at(2,2)+t.at(2,3,2)*s.at(3,2)
-t.at(2,0,3)*s.at(0,3)+t.at(2,1,3)*s.at(1,3)+t.at(2,2,3)*s.at(2,3)+t.at(2,3,3)*s.at(3,3),
t.at(3,0,0)*s.at(0,0)-t.at(3,1,0)*s.at(1,0)-t.at(3,2,0)*s.at(2,0)-t.at(3,3,0)*s.at(3,0)
-t.at(3,0,1)*s.at(0,1)+t.at(3,1,1)*s.at(1,1)+t.at(3,2,1)*s.at(2,1)+t.at(3,3,1)*s.at(3,1)
-t.at(3,0,2)*s.at(0,2)+t.at(3,1,2)*s.at(1,2)+t.at(3,2,2)*s.at(2,2)+t.at(3,3,2)*s.at(3,2)
-t.at(3,0,3)*s.at(0,3)+t.at(3,1,3)*s.at(1,3)+t.at(3,2,3)*s.at(2,3)+t.at(3,3,3)*s.at(3,3));
}
template<typename Scal1, typename Scal2>
Lorentz_Ten3<PROMOTE(Scal1,Scal2)>
ATOOLS::Contraction(const Lorentz_Ten3<Scal1>& ten, int i,
const Lorentz_Ten2<Scal2>& sen, int j) {
Lorentz_Ten3<Scal1> t; Lorentz_Ten2<Scal2> s;
// T^{mu,nu,rho} = t^{mu,nu,alpha} s^{rho,alpha'} g_{alpha,alpha'}
if ((i==3) && (j==2))
{t = ten; s = sen;}
else if ((i==2) && (j==2))
{t = ten.Transpose(2,3); s = sen;}
else if ((i==1) && (j==2))
{t = ten.Transpose(1,2); t = t.Transpose(2,3); s = sen;}
else if ((i==3) && (j==1))
{t = ten; s = sen.Transpose();}
else if ((i==2) && (j==1))
{t = ten.Transpose(2,3); s = sen.Transpose();}
else if ((i==1) && (j==1))
{t = ten.Transpose(1,2); t = t.Transpose(2,3); s = sen.Transpose();}
else
return Lorentz_Ten3<PROMOTE(Scal1,Scal2)>();
return Lorentz_Ten3<PROMOTE(Scal1,Scal2)>(
t.at(0,0,0)*s.at(0,0)-t.at(0,0,1)*s.at(0,1)
-t.at(0,0,2)*s.at(0,2)-t.at(0,0,3)*s.at(0,3),
t.at(1,0,0)*s.at(0,0)-t.at(1,0,1)*s.at(0,1)
-t.at(1,0,2)*s.at(0,2)-t.at(1,0,3)*s.at(0,3),
t.at(2,0,0)*s.at(0,0)-t.at(2,0,1)*s.at(0,1)
-t.at(2,0,2)*s.at(0,2)-t.at(2,0,3)*s.at(0,3),
t.at(3,0,0)*s.at(0,0)-t.at(3,0,1)*s.at(0,1)
-t.at(3,0,2)*s.at(0,2)-t.at(3,0,3)*s.at(0,3),
t.at(0,1,0)*s.at(0,0)-t.at(0,1,1)*s.at(0,1)
-t.at(0,1,2)*s.at(0,2)-t.at(0,1,3)*s.at(0,3),
t.at(1,1,0)*s.at(0,0)-t.at(1,1,1)*s.at(0,1)
-t.at(1,1,2)*s.at(0,2)-t.at(1,1,3)*s.at(0,3),
t.at(2,1,0)*s.at(0,0)-t.at(2,1,1)*s.at(0,1)
-t.at(2,1,2)*s.at(0,2)-t.at(2,1,3)*s.at(0,3),
t.at(3,1,0)*s.at(0,0)-t.at(3,1,1)*s.at(0,1)
-t.at(3,1,2)*s.at(0,2)-t.at(3,1,3)*s.at(0,3),
t.at(0,2,0)*s.at(0,0)-t.at(0,2,1)*s.at(0,1)
-t.at(0,2,2)*s.at(0,2)-t.at(0,2,3)*s.at(0,3),
t.at(1,2,0)*s.at(0,0)-t.at(1,2,1)*s.at(0,1)
-t.at(1,2,2)*s.at(0,2)-t.at(1,2,3)*s.at(0,3),
t.at(2,2,0)*s.at(0,0)-t.at(2,2,1)*s.at(0,1)
-t.at(2,2,2)*s.at(0,2)-t.at(2,2,3)*s.at(0,3),
t.at(3,2,0)*s.at(0,0)-t.at(3,2,1)*s.at(0,1)
-t.at(3,2,2)*s.at(0,2)-t.at(3,2,3)*s.at(0,3),
t.at(0,3,0)*s.at(0,0)-t.at(0,3,1)*s.at(0,1)
-t.at(0,3,2)*s.at(0,2)-t.at(0,3,3)*s.at(0,3),
t.at(1,3,0)*s.at(0,0)-t.at(1,3,1)*s.at(0,1)
-t.at(1,3,2)*s.at(0,2)-t.at(1,3,3)*s.at(0,3),
t.at(2,3,0)*s.at(0,0)-t.at(2,3,1)*s.at(0,1)
-t.at(2,3,2)*s.at(0,2)-t.at(2,3,3)*s.at(0,3),
t.at(3,3,0)*s.at(0,0)-t.at(3,3,1)*s.at(0,1)
-t.at(3,3,2)*s.at(0,2)-t.at(3,3,3)*s.at(0,3),
t.at(0,0,0)*s.at(1,0)-t.at(0,0,1)*s.at(1,1)
-t.at(0,0,2)*s.at(1,2)-t.at(0,0,3)*s.at(1,3),
t.at(1,0,0)*s.at(1,0)-t.at(1,0,1)*s.at(1,1)
-t.at(1,0,2)*s.at(1,2)-t.at(1,0,3)*s.at(1,3),
t.at(2,0,0)*s.at(1,0)-t.at(2,0,1)*s.at(1,1)
-t.at(2,0,2)*s.at(1,2)-t.at(2,0,3)*s.at(1,3),
t.at(3,0,0)*s.at(1,0)-t.at(3,0,1)*s.at(1,1)
-t.at(3,0,2)*s.at(1,2)-t.at(3,0,3)*s.at(1,3),
t.at(0,1,0)*s.at(1,0)-t.at(0,1,1)*s.at(1,1)
-t.at(0,1,2)*s.at(1,2)-t.at(0,1,3)*s.at(1,3),
t.at(1,1,0)*s.at(1,0)-t.at(1,1,1)*s.at(1,1)
-t.at(1,1,2)*s.at(1,2)-t.at(1,1,3)*s.at(1,3),
t.at(2,1,0)*s.at(1,0)-t.at(2,1,1)*s.at(1,1)
-t.at(2,1,2)*s.at(1,2)-t.at(2,1,3)*s.at(1,3),
t.at(3,1,0)*s.at(1,0)-t.at(3,1,1)*s.at(1,1)
-t.at(3,1,2)*s.at(1,2)-t.at(3,1,3)*s.at(1,3),
t.at(0,2,0)*s.at(1,0)-t.at(0,2,1)*s.at(1,1)
-t.at(0,2,2)*s.at(1,2)-t.at(0,2,3)*s.at(1,3),
t.at(1,2,0)*s.at(1,0)-t.at(1,2,1)*s.at(1,1)
-t.at(1,2,2)*s.at(1,2)-t.at(1,2,3)*s.at(1,3),
t.at(2,2,0)*s.at(1,0)-t.at(2,2,1)*s.at(1,1)
-t.at(2,2,2)*s.at(1,2)-t.at(2,2,3)*s.at(1,3),
t.at(3,2,0)*s.at(1,0)-t.at(3,2,1)*s.at(1,1)
-t.at(3,2,2)*s.at(1,2)-t.at(3,2,3)*s.at(1,3),
t.at(0,3,0)*s.at(1,0)-t.at(0,3,1)*s.at(1,1)
-t.at(0,3,2)*s.at(1,2)-t.at(0,3,3)*s.at(1,3),
t.at(1,3,0)*s.at(1,0)-t.at(1,3,1)*s.at(1,1)
-t.at(1,3,2)*s.at(1,2)-t.at(1,3,3)*s.at(1,3),
t.at(2,3,0)*s.at(1,0)-t.at(2,3,1)*s.at(1,1)
-t.at(2,3,2)*s.at(1,2)-t.at(2,3,3)*s.at(1,3),
t.at(3,3,0)*s.at(1,0)-t.at(3,3,1)*s.at(1,1)
-t.at(3,3,2)*s.at(1,2)-t.at(3,3,3)*s.at(1,3),
t.at(0,0,0)*s.at(2,0)-t.at(0,0,1)*s.at(2,1)
-t.at(0,0,2)*s.at(2,2)-t.at(0,0,3)*s.at(2,3),
t.at(1,0,0)*s.at(2,0)-t.at(1,0,1)*s.at(2,1)
-t.at(1,0,2)*s.at(2,2)-t.at(1,0,3)*s.at(2,3),
t.at(2,0,0)*s.at(2,0)-t.at(2,0,1)*s.at(2,1)
-t.at(2,0,2)*s.at(2,2)-t.at(2,0,3)*s.at(2,3),
t.at(3,0,0)*s.at(2,0)-t.at(3,0,1)*s.at(2,1)
-t.at(3,0,2)*s.at(2,2)-t.at(3,0,3)*s.at(2,3),
t.at(0,1,0)*s.at(2,0)-t.at(0,1,1)*s.at(2,1)
-t.at(0,1,2)*s.at(2,2)-t.at(0,1,3)*s.at(2,3),
t.at(1,1,0)*s.at(2,0)-t.at(1,1,1)*s.at(2,1)
-t.at(1,1,2)*s.at(2,2)-t.at(1,1,3)*s.at(2,3),
t.at(2,1,0)*s.at(2,0)-t.at(2,1,1)*s.at(2,1)
-t.at(2,1,2)*s.at(2,2)-t.at(2,1,3)*s.at(2,3),
t.at(3,1,0)*s.at(2,0)-t.at(3,1,1)*s.at(2,1)
-t.at(3,1,2)*s.at(2,2)-t.at(3,1,3)*s.at(2,3),
t.at(0,2,0)*s.at(2,0)-t.at(0,2,1)*s.at(2,1)
-t.at(0,2,2)*s.at(2,2)-t.at(0,2,3)*s.at(2,3),
t.at(1,2,0)*s.at(2,0)-t.at(1,2,1)*s.at(2,1)
-t.at(1,2,2)*s.at(2,2)-t.at(1,2,3)*s.at(2,3),
t.at(2,2,0)*s.at(2,0)-t.at(2,2,1)*s.at(2,1)
-t.at(2,2,2)*s.at(2,2)-t.at(2,2,3)*s.at(2,3),
t.at(3,2,0)*s.at(2,0)-t.at(3,2,1)*s.at(2,1)
-t.at(3,2,2)*s.at(2,2)-t.at(3,2,3)*s.at(2,3),
t.at(0,3,0)*s.at(2,0)-t.at(0,3,1)*s.at(2,1)
-t.at(0,3,2)*s.at(2,2)-t.at(0,3,3)*s.at(2,3),
t.at(1,3,0)*s.at(2,0)-t.at(1,3,1)*s.at(2,1)
-t.at(1,3,2)*s.at(2,2)-t.at(1,3,3)*s.at(2,3),
t.at(2,3,0)*s.at(2,0)-t.at(2,3,1)*s.at(2,1)
-t.at(2,3,2)*s.at(2,2)-t.at(2,3,3)*s.at(2,3),
t.at(3,3,0)*s.at(2,0)-t.at(3,3,1)*s.at(2,1)
-t.at(3,3,2)*s.at(2,2)-t.at(3,3,3)*s.at(2,3),
t.at(0,0,0)*s.at(3,0)-t.at(0,0,1)*s.at(3,1)
-t.at(0,0,2)*s.at(3,2)-t.at(0,0,3)*s.at(3,3),
t.at(1,0,0)*s.at(3,0)-t.at(1,0,1)*s.at(3,1)
-t.at(1,0,2)*s.at(3,2)-t.at(1,0,3)*s.at(3,3),
t.at(2,0,0)*s.at(3,0)-t.at(2,0,1)*s.at(3,1)
-t.at(2,0,2)*s.at(3,2)-t.at(2,0,3)*s.at(3,3),
t.at(3,0,0)*s.at(3,0)-t.at(3,0,1)*s.at(3,1)
-t.at(3,0,2)*s.at(3,2)-t.at(3,0,3)*s.at(3,3),
t.at(0,1,0)*s.at(3,0)-t.at(0,1,1)*s.at(3,1)
-t.at(0,1,2)*s.at(3,2)-t.at(0,1,3)*s.at(3,3),
t.at(1,1,0)*s.at(3,0)-t.at(1,1,1)*s.at(3,1)
-t.at(1,1,2)*s.at(3,2)-t.at(1,1,3)*s.at(3,3),
t.at(2,1,0)*s.at(3,0)-t.at(2,1,1)*s.at(3,1)
-t.at(2,1,2)*s.at(3,2)-t.at(2,1,3)*s.at(3,3),
t.at(3,1,0)*s.at(3,0)-t.at(3,1,1)*s.at(3,1)
-t.at(3,1,2)*s.at(3,2)-t.at(3,1,3)*s.at(3,3),
t.at(0,2,0)*s.at(3,0)-t.at(0,2,1)*s.at(3,1)
-t.at(0,2,2)*s.at(3,2)-t.at(0,2,3)*s.at(3,3),
t.at(1,2,0)*s.at(3,0)-t.at(1,2,1)*s.at(3,1)
-t.at(1,2,2)*s.at(3,2)-t.at(1,2,3)*s.at(3,3),
t.at(2,2,0)*s.at(3,0)-t.at(2,2,1)*s.at(3,1)
-t.at(2,2,2)*s.at(3,2)-t.at(2,2,3)*s.at(3,3),
t.at(3,2,0)*s.at(3,0)-t.at(3,2,1)*s.at(3,1)
-t.at(3,2,2)*s.at(3,2)-t.at(3,2,3)*s.at(3,3),
t.at(0,3,0)*s.at(3,0)-t.at(0,3,1)*s.at(3,1)
-t.at(0,3,2)*s.at(3,2)-t.at(0,3,3)*s.at(3,3),
t.at(1,3,0)*s.at(3,0)-t.at(1,3,1)*s.at(3,1)
-t.at(1,3,2)*s.at(3,2)-t.at(1,3,3)*s.at(3,3),
t.at(2,3,0)*s.at(3,0)-t.at(2,3,1)*s.at(3,1)
-t.at(2,3,2)*s.at(3,2)-t.at(2,3,3)*s.at(3,3),
t.at(3,3,0)*s.at(3,0)-t.at(3,3,1)*s.at(3,1)
-t.at(3,3,2)*s.at(3,2)-t.at(3,3,3)*s.at(3,3));
}
template<typename Scal1, typename Scal2>
Vec4<PROMOTE(Scal1,Scal2)>
ATOOLS::Contraction(const Lorentz_Ten4<Scal1>& ten, int i, int j, int k,
const Lorentz_Ten3<Scal2>& sen, int l, int m, int n) {
// v^{mu} = t^{mu,alpha,beta,gamma} s^{alpha',beta',gamma'}
// g_{alpha,alpha'} g_{beta,beta'} g_{gamma,gamma'}
if ((i!=2) && (j!=3) && (k!=4) && (l!=2) && (m!=3) && (n!=4)) {
THROW(fatal_error,"tensor contraction not specified for these indizes ...");
}
Vec4<PROMOTE(Scal1,Scal2)> end;
for (unsigned short int I=0; I<4; ++I) {
PROMOTE(Scal1,Scal2) res(ten.at(0,0,0,0)-ten.at(0,0,0,0));
for (unsigned short int a=0; a<4; ++a)
for (unsigned short int b=0; b<4; ++b)
for (unsigned short int c=0; c<4; ++c) {
double ga(-1.),gb(-1.),gc(-1.);
if (a==0) ga=1.; if (b==0) gb=1.; if (c==0) gc=1.;
res += ga*gb*gc*ten.at(I,a,b,c)*sen.at(a,b,c);
}
end[I] = res;
}
return end;
}
template<typename Scal1, typename Scal2>
Lorentz_Ten3<PROMOTE(Scal1,Scal2)>
ATOOLS::Contraction(const Lorentz_Ten4<Scal1>& ten, int i, int j,
const Lorentz_Ten3<Scal2>& sen, int k, int l) {
// T^{mu,nu,rho} = t^{mu,nu,alpha,beta} s^{rho,alpha',beta'}
// g_{alpha,alpha'} g_{beta,beta'}
if ((i!=3) && (j!=4) && (k!=3) && (l!=4)) {
THROW(fatal_error,"tensor contraction not specified for these indizes ...");
}
Vec4<PROMOTE(Scal1,Scal2)> end;
for (unsigned short int I=0; I<4; ++I)
for (unsigned short int J=0; J<4; ++J)
for (unsigned short int K=0; K<4; ++K) {
PROMOTE(Scal1,Scal2) res(ten.at(0,0,0,0)-ten.at(0,0,0,0));
for (unsigned short int a=0; a<4; ++a)
for (unsigned short int b=0; b<4; ++b) {
double ga(-1.),gb(-1.);
if (a==0) ga=1.; if (b==0) gb=1.;
res += ga*gb*ten.at(I,J,a,b)*sen.at(K,a,b);
}
end.at(I,J,K) = res;
}
return end;
}
template<typename Scal1, typename Scal2>
PROMOTE(Scal1,Scal2)
ATOOLS::Contraction(const Lorentz_Ten4<Scal1>& ten, int i, int j, int k, int l,
const Lorentz_Ten4<Scal2>& sen, int m, int n, int o, int p) {
Lorentz_Ten4<Scal1> t; Lorentz_Ten4<Scal2> s;
// c = t^{alpha,beta,gamma,delta} s^{alpha',beta',gamma',delta'}
// g_{alpha,alpha'} g_{beta,beta'} g_{gamma.gamma'} g_{delta,delta'}
if ((i==m) && (j==n) && (k==o) && (l==p))
{t=ten; s=sen;}
else if ((i==m) && (j==n) && (k==p) && (l==o))
{t=ten; s=sen.Transpose(3,4);}
else if ((i==m) && (j==o) && (k==n) && (l==p))
{t=ten; s=sen.Transpose(2,3);}
else if ((i==m) && (j==o) && (k==p) && (l==n))
{t=ten; s=sen.Transpose(3,4); s=s.Transpose(2,3);}
else if ((i==n) && (j==m) && (k==o) && (l==p))
{t=ten; s=sen.Transpose(1,2);}
else if ((i==n) && (j==m) && (k==p) && (l==o))
{t=ten; s=sen.Transpose(1,2); s=s.Transpose(3,4);}
else if ((i==n) && (j==o) && (k==m) && (l==p))
{t=ten; s=sen.Transpose(1,3); s=s.Transpose(2,3);}
else if ((i==n) && (j==o) && (k==p) && (l==m))
{t=ten; s=sen.Transpose(3,4); s=s.Transpose(2,3); s=s.Transpose(1,2);}
else if ((i==n) && (j==p) && (k==m) && (l==o))
{t=ten; s=sen.Transpose(1,3); s=s.Transpose(2,4); s=s.Transpose(2,3);}
else if ((i==n) && (j==p) && (k==o) && (l==m))
{t=ten; s=sen.Transpose(2,4); s=s.Transpose(1,2);}
else if ((i==o) && (j==m) && (k==n) && (l==p))
{t=ten; s=sen.Transpose(1,2); s=s.Transpose(2,3);}
else if ((i==o) && (j==m) && (k==p) && (l==n))
{t=ten; s=sen.Transpose(1,2); s=s.Transpose(2,3); s=s.Transpose(2,4);}
else if ((i==o) && (j==n) && (k==m) && (l==p))
{t=ten; s=sen.Transpose(1,3);}
else if ((i==o) && (j==n) && (k==p) && (l==m))
{t=ten; s=sen.Transpose(1,4); s=s.Transpose(3,4);}
else if ((i==o) && (j==p) && (k==m) && (l==n))
{t=ten; s=sen.Transpose(1,3); s=s.Transpose(2,4);}
else if ((i==o) && (j==p) && (k==n) && (l==m))
{t=ten; s=sen.Transpose(1,4); s=s.Transpose(2,3); s=s.Transpose(3,4);}
else if ((i==p) && (j==m) && (k==n) && (l==o))
{t=ten; s=sen.Transpose(1,2); s=s.Transpose(2,3); s=s.Transpose(3,4);}
else if ((i==p) && (j==m) && (k==o) && (l==n))
{t=ten; s=sen.Transpose(1,2); s=s.Transpose(2,4);}
else if ((i==p) && (j==n) && (k==m) && (l==o))
{t=ten; s=sen.Transpose(1,3); s=s.Transpose(3,4);}
else if ((i==p) && (j==n) && (k==o) && (l==m))
{t=ten; s=sen.Transpose(1,4);}
else if ((i==p) && (j==o) && (k==m) && (l==n))
{t=ten; s=sen.Transpose(1,3); s=s.Transpose(2,4); s=s.Transpose(3,4);}
else if ((i==p) && (j==o) && (k==n) && (l==m))
{t=ten; s=sen.Transpose(1,4); s=s.Transpose(2,3);}
else
return PROMOTE(Scal1,Scal2)();
PROMOTE(Scal1,Scal2) res(ten.at(0,0,0,0)-ten.at(0,0,0,0));
for (unsigned short int a=0; a<4; ++a)
for (unsigned short int b=0; b<4; ++b)
for (unsigned short int c=0; c<4; ++c)
for (unsigned short int d=0; d<4; ++d) {
double ga(-1.),gb(-1.),gc(-1.),gd(-1.);
if (a==0) ga=1.; if (b==0) gb=1.; if (c==0) gc=1.; if (d==0) gd=1.;
res += ga*gb*gc*gd*t.at(a,b,c,d)*s.at(a,b,c,d);
}
return res;
}
template<typename Scal1, typename Scal2>
Lorentz_Ten2<PROMOTE(Scal1,Scal2)>
ATOOLS::Contraction(const Lorentz_Ten4<Scal1>& ten, int i, int j, int k,
const Lorentz_Ten4<Scal2>& sen, int l, int m, int n) {
Lorentz_Ten4<Scal1> t; Lorentz_Ten4<Scal2> s;
// T^{mu,nu} = t^{mu,alpha,beta,gamma} s^{nu,alpha',beta',gamma'}
// g_{alpha,alpha'} g_{beta,beta'} g_{gamma,gamma'}
if (((i==l) && (j==m) && (k==n)) && ((i==2) && (j==3) && (k==4)))
{t=ten; s=sen;}
else if (((i==l) && (j==m) && (k==n)) && ((i==1) && (j==3) && (k==4)))
{t=ten.Transpose(1,2); s=sen.Transpose(1,2);}
else if (((i==l) && (j==m) && (k==n)) && ((i==1) && (j==2) && (k==4)))
{t=ten.Transpose(1,3); s=sen.Transpose(1,3);}
else if (((i==l) && (j==m) && (k==n)) && ((i==1) && (j==2) && (k==3)))
{t=ten.Transpose(1,4); s=sen.Transpose(1,4);}
else if (((i==l) && (j==n) && (k==m)) && ((i==2) && (j==3) && (k==4)))
{t=ten; s=sen.Transpose(3,4);}
else if (((i==l) && (j==n) && (k==m)) && ((i==1) && (j==3) && (k==4)))
{t=ten.Transpose(1,2); s=sen.Transpose(1,2); s=s.Transpose(3,4);}
else if (((i==l) && (j==n) && (k==m)) && ((i==1) && (j==2) && (k==4)))
{t=ten.Transpose(1,3); s=sen.Transpose(1,3); s=s.Transpose(2,4);}
else if (((i==l) && (j==n) && (k==m)) && ((i==1) && (j==2) && (k==3)))
{t=ten.Transpose(1,4); s=sen.Transpose(1,4); s=s.Transpose(2,3);}
else if (((i==m) && (j==l) && (k==n)) && ((i==2) && (j==3) && (k==4)))
{t=ten; s=sen.Transpose(2,3);}
else if (((i==m) && (j==l) && (k==n)) && ((i==1) && (j==3) && (k==4)))
{t=ten.Transpose(1,2); s=sen.Transpose(1,2); s=s.Transpose(2,3);}
else if (((i==m) && (j==l) && (k==n)) && ((i==1) && (j==2) && (k==4)))
{t=ten.Transpose(1,3); s=sen.Transpose(1,3); s=s.Transpose(2,3);}
else if (((i==m) && (j==l) && (k==n)) && ((i==1) && (j==2) && (k==3)))
{t=ten.Transpose(1,4); s=sen.Transpose(1,4); s=s.Transpose(2,4);}
else if (((i==m) && (j==n) && (k==l)) && ((i==2) && (j==3) && (k==4)))
{t=ten; s=sen.Transpose(2,4); s=s.Transpose(3,4);}
else if (((i==m) && (j==n) && (k==l)) && ((i==1) && (j==3) && (k==4)))
{t=ten.Transpose(1,2); s=sen.Transpose(1,2); s=s.Transpose(2,4); s=s.Transpose(3,4);}
else if (((i==m) && (j==n) && (k==l)) && ((i==1) && (j==2) && (k==4)))
{t=ten.Transpose(1,3); s=sen.Transpose(1,3); s=s.Transpose(2,4); s=s.Transpose(2,3);}
else if (((i==m) && (j==n) && (k==l)) && ((i==1) && (j==2) && (k==3)))
{t=ten.Transpose(1,4); s=sen.Transpose(1,4); s=s.Transpose(2,4); s=s.Transpose(3,4);}
else if (((i==n) && (j==l) && (k==m)) && ((i==2) && (j==3) && (k==4)))
{t=ten; s=sen.Transpose(2,3); s=s.Transpose(3,4);}
else if (((i==n) && (j==l) && (k==m)) && ((i==1) && (j==3) && (k==4)))
{t=ten.Transpose(1,2); s=sen.Transpose(1,2); s=s.Transpose(2,3); s=s.Transpose(3,4);}
else if (((i==n) && (j==l) && (k==m)) && ((i==1) && (j==2) && (k==4)))
{t=ten.Transpose(1,3); s=sen.Transpose(1,3); s=s.Transpose(3,4); s=s.Transpose(2,3);}
else if (((i==n) && (j==l) && (k==m)) && ((i==1) && (j==2) && (k==3)))
{t=ten.Transpose(1,4); s=sen.Transpose(1,4); s=s.Transpose(2,3); s=s.Transpose(3,4);}
else if (((i==n) && (j==m) && (k==l)) && ((i==2) && (j==3) && (k==4)))
{t=ten; s=sen.Transpose(2,4);}
else if (((i==n) && (j==m) && (k==l)) && ((i==1) && (j==3) && (k==4)))
{t=ten.Transpose(1,2); s=sen.Transpose(1,2); s=s.Transpose(2,4);}
else if (((i==n) && (j==m) && (k==l)) && ((i==1) && (j==2) && (k==4)))
{t=ten.Transpose(1,3); s=sen.Transpose(1,3); s=s.Transpose(3,4);}
else if (((i==n) && (j==m) && (k==l)) && ((i==1) && (j==2) && (k==3)))
{t=ten.Transpose(1,4); s=sen.Transpose(1,4); s=s.Transpose(3,4);}
else {
THROW(fatal_error,"tensor contraction not specified for these indizes ...");
return Lorentz_Ten2<PROMOTE(Scal1,Scal2)>();
}
Lorentz_Ten2<PROMOTE(Scal1,Scal2)> end;
for (unsigned short int I=0; I<4; ++I)
for (unsigned short int J=0; J<4; ++J) {
PROMOTE(Scal1,Scal2) res(ten.at(0,0,0,0)-ten.at(0,0,0,0));
for (unsigned short int a=0; a<4; ++a)
for (unsigned short int b=0; b<4; ++b)
for (unsigned short int c=0; c<4; ++c) {
double ga(-1.),gb(-1.),gc(-1.);
if (a==0) ga=1.; if (b==0) gb=1.; if (c==0) gc=1.;
res += ga*gb*gc*t.at(I,a,b,c)*s.at(J,a,b,c);
}
end.at(I,J) = res;
}
return end;
}
template<typename Scal1, typename Scal2>
Lorentz_Ten4<PROMOTE(Scal1,Scal2)>
ATOOLS::Contraction(const Lorentz_Ten4<Scal1>& ten, int i, int j,
const Lorentz_Ten4<Scal2>& sen, int k, int l) {
// T^{mu,nu,rho,sigma} = t^{mu,nu,alpha,beta} s^{rho,sigma,alpha',beta'}
// g_{alpha,alpha'} g_{beta,beta'}
if ((i!=3) && (j!=4) && (k!=3) && (l!=4)) {
THROW(fatal_error,"tensor contraction not specified for these indizes ...");
return Lorentz_Ten4<PROMOTE(Scal1,Scal2)>();
}
Lorentz_Ten4<PROMOTE(Scal1,Scal2)> end;
for (unsigned short int I=0; I<4; ++I)
for (unsigned short int J=0; J<4; ++J)
for (unsigned short int K=0; K<4; ++K)
for (unsigned short int L=0; L<4; ++L) {
PROMOTE(Scal1,Scal2) res(ten.at(0,0,0,0)-ten.at(0,0,0,0));
for (unsigned short int a=0; a<4; ++a)
for (unsigned short int b=0; b<4; ++b) {
double ga(-1.),gb(-1.);
if (a==0) ga=1.; if (b==0) gb=1.;
res += ga*gb*ten.at(I,J,a,b)*sen.at(K,L,a,b);
}
end.at(I,J) = res;
}
return end;
}
#endif