ATOOLS/Math/Tensor_Contractions_Epsilon.H

0001 #ifndef ATOOLS_Math_Tensor_Contractions_Epsilon_H
0002 #define ATOOLS_Math_Tensor_Contractions_Epsilon_H
0003
0004 /********************************************************************************************
0005 ******                                                                                 ******
0006 ******       epsilon tensor contractions explicitely                                   ******
0007 ******                                                                                 ******
0008 ********************************************************************************************/
0009
0010 template<typename Scalar>
0011 Lorentz_Ten3<Scalar>
0012 ATOOLS::EpsilonTensorContraction(const Vec4<Scalar>& v) {
0013   // T^{mu,nu,rho} = eps^{mu,nu,rho,alpha} g_{alpha,alpha'} v^{alpha}
0014   return Lorentz_Ten3<Scalar>(
0015              0.  , 0.  , 0.  , 0.  ,
0016              0.  , 0.  ,-v[3], v[2],
0017              0.  , v[3], 0.  ,-v[1],
0018              0.  ,-v[2], v[1], 0.  ,
0019              0.  , 0.  , v[3],-v[2],
0020              0.  , 0.  , 0.  , 0.  ,
0021             -v[3], 0.  , 0.  , v[0],
0022              v[2], 0.  ,-v[0], 0.  ,
0023              0.  ,-v[3], 0.  , v[1],
0024              v[3], 0.  , 0.  ,-v[0],
0025              0.  , 0.  , 0.  , 0.  ,
0026             -v[1], v[0], 0.  , 0.  ,
0027              0.  , v[2],-v[1], 0.  ,
0028             -v[2], 0.  , v[0], 0.  ,
0029              v[1],-v[0], 0.  , 0.  ,
0030              0.  , 0.  , 0.  , 0.  );
0031 }
0032
0033 template<typename Scal1, typename Scal2>
0034 Lorentz_Ten2<PROMOTE(Scal1,Scal2)>
0035 ATOOLS::EpsilonTensorContraction(const Vec4<Scal1>& p, const Vec4<Scal2>& q) {
0036   // returns epsilon tensor contractions
0037   //  T^{mu,nu} = eps^{mu,nu,alpha,beta}g_{alpha,alpha'}g_{beta,beta'}p^alpha' q^beta'
0038   return Lorentz_Ten2<PROMOTE(Scal1,Scal2)>(
0039                                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] ,
0040                 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] ,
0041                -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] ,
0042                 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.   );
0043 }
0044
0045 template<typename Scal1, typename Scal2, typename Scal3>
0046 Vec4<PROMOTE(Scal1,PROMOTE(Scal2,Scal3))>
0047 ATOOLS::EpsilonTensorContraction(const Vec4<Scal1>& p1, const Vec4<Scal2>& p2,
0048                                  const Vec4<Scal3>& p3) {
0049   // v^{mu} = eps^{mu,alpha,beta,gamma}g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
0050   //              p1^{alpha'} p2^{beta'} p3^{gamma'} and permutations thereof
0051 // maybe uncomment, but doubtful if any computational improvement
0052 //   if (IsEqual(p1,p2) || IsEqual(p1,p3)  || IsEqual(p2,p3))
0053 //     return Vec4<PROMOTE(Scal1,PROMOTE(Scal2,Scal3))>(0.,0.,0.,0.);
0054   return EpsilonTensorContraction(BuildTensor(p1,p2,p3),2,3,4);
0055 }
0056
0057 template<typename Scal1, typename Scal2, typename Scal3, typename Scal4>
0058 PROMOTE(Scal1,PROMOTE(Scal2,PROMOTE(Scal3,Scal4)))
0059 ATOOLS::EpsilonTensorContraction(const Vec4<Scal1>& p1, const Vec4<Scal2>& p2,
0060                                  const Vec4<Scal3>& p3, const Vec4<Scal4>& p4) {
0061   // c = eps^{alpha,beta,gamma,delta} p1^alpha p2^beta p3^gamma p4^delta
0062   //                          gamma_{alpha,alpha'} gamma_{beta,beta'}
0063   //                          gamma_{gamma,gamma'} gamma_{delta,delta'}
0064   return   p1[0]*p2[1]*p3[2]*p4[3]
0065           -p1[0]*p2[1]*p3[3]*p4[2]
0066           -p1[0]*p2[2]*p3[1]*p4[3]
0067           +p1[0]*p2[2]*p3[3]*p4[1]
0068           +p1[0]*p2[3]*p3[1]*p4[2]
0069           -p1[0]*p2[3]*p3[2]*p4[1]
0070
0071           -p1[1]*p2[0]*p3[2]*p4[3]
0072           +p1[1]*p2[0]*p3[3]*p4[2]
0073           +p1[1]*p2[2]*p3[0]*p4[3]
0074           -p1[1]*p2[2]*p3[3]*p4[0]
0075           -p1[1]*p2[3]*p3[0]*p4[2]
0076           +p1[1]*p2[3]*p3[2]*p4[0]
0077
0078           +p1[2]*p2[0]*p3[1]*p4[3]
0079           -p1[2]*p2[0]*p3[3]*p4[1]
0080           -p1[2]*p2[1]*p3[0]*p4[3]
0081           +p1[2]*p2[1]*p3[3]*p4[0]
0082           +p1[2]*p2[3]*p3[0]*p4[1]
0083           -p1[2]*p2[3]*p3[1]*p4[0]
0084
0085           -p1[3]*p2[0]*p3[1]*p4[2]
0086           +p1[3]*p2[0]*p3[2]*p4[1]
0087           +p1[3]*p2[1]*p3[0]*p4[2]
0088           -p1[3]*p2[1]*p3[2]*p4[0]
0089           -p1[3]*p2[2]*p3[0]*p4[1]
0090           +p1[3]*p2[2]*p3[1]*p4[0];
0091 }
0092
0093 template<typename Scalar>
0094 Lorentz_Ten2<Scalar>
0095 ATOOLS::EpsilonTensorContraction(const Lorentz_Ten2<Scalar>& t, int i, int j) {
0096 // returns the contraction eps^{mu,nu,alpha,beta}g_{alpha,alpha'}g_{beta,beta'}t^{alpha'beta'}
0097   if ((i==3) && (j==4))
0098     return Lorentz_Ten2<Scalar>(
0099       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) ,
0100       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) ,
0101      -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) ,
0102       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. );
0103   else if ((j==3) && (i==4))
0104     return -Lorentz_Ten2<Scalar>(
0105       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) ,
0106       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) ,
0107      -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) ,
0108       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. );
0109   return Lorentz_Ten2<Scalar>();
0110 }
0111
0112 template<typename Scal1, typename Scal2>
0113 Vec4<PROMOTE(Scal1,Scal2)>
0114 ATOOLS::EpsilonTensorContraction(const Lorentz_Ten2<Scal1>& t, int i, int j,
0115                                  const Vec4<Scal2>& p, int k) {
0116   // v^{mu} = eps^{mu,alpha,beta,gamma}g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
0117   //              t^{alpha',beta'} p^{gamma'} and permutations thereof
0118   return EpsilonTensorContraction(BuildTensor(t,p),i,j,k);
0119 }
0120
0121 template<typename Scalar>
0122 Vec4<Scalar>
0123 ATOOLS::EpsilonTensorContraction(const Lorentz_Ten3<Scalar>& ten, int i, int j, int k) {
0124   Lorentz_Ten3<Scalar> t;
0125   // v^{mu} = eps^{mu,alpha,beta,gamma} t^{alpha',beta',gamma'} g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
0126   if      ((i==2) && (j==3) && (k==4))
0127      t = ten;
0128   // v^{mu} = eps^{mu,alpha,beta,gamma} t^{alpha',gamma',beta'} g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
0129   else if ((i==2) && (j==4) && (k==3))
0130      t = ten.Transpose(2,3);
0131   // v^{mu} = eps^{mu,alpha,beta,gamma} t^{beta',alpha',gamma'} g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
0132   else if ((i==3) && (j==2) && (k==4))
0133      t = ten.Transpose(1,2);
0134   // v^{mu} = eps^{mu,alpha,beta,gamma} t^{beta',gamma',alpha'} g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
0135   else if ((i==3) && (j==4) && (k==2))
0136     {t = ten.Transpose(2,3); t = t.Transpose(1,2);}
0137   // v^{mu} = eps^{mu,alpha,beta,gamma} t^{gamma',alpha',beta'} g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
0138   else if ((i==4) && (j==2) && (k==3))
0139     {t = ten.Transpose(1,2); t = t.Transpose(2,3);}
0140   // v^{mu} = eps^{mu,alpha,beta,gamma} t^{gamma',beta',alpha'} g_{alpha,alpha'}g_{beta,beta'}g_{gamma,gamma'}
0141   else if ((i==4) && (j==3) && (k==2))
0142      t = ten.Transpose(1,3);
0143   else
0144     return Vec4<Scalar>(0.,0.,0.,0.);
0145   return Vec4<Scalar>(
0146                             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),
0147     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),
0148     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),
0149     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)                        );
0150 }
0151
0152 template<typename Scal1, typename Scal2>
0153 PROMOTE(Scal1,Scal2)
0154 ATOOLS::EpsilonTensorContraction(const Lorentz_Ten2<Scal1>& ten, int i, int j,
0155                                  const Lorentz_Ten2<Scal2>& sen, int k, int l) {
0156   // c = eps^{alpha,beta,gamma,delta} t^{alpha',beta'} s^{gamma',delta'}
0157   //                                  gamma_{alpha,alpha'} gamma_{beta,beta'}
0158   //                                  gamma_{gamma,gamma'} gamma_{delta,delta'}
0159   return EpsilonTensorContraction(BuildTensor(ten,sen),i,j,k,l);
0160 }
0161
0162 template<typename Scal1, typename Scal2>
0163 PROMOTE(Scal1,Scal2)
0164 ATOOLS::EpsilonTensorContraction(const Lorentz_Ten3<Scal1>& ten, int i, int j, int k,
0165                                  const Vec4<Scal2>& v, int l) {
0166   // c = eps^{alpha,beta,gamma,delta} t^{alpha',beta',gamma'} v^delta'
0167   //                                  gamma_{alpha,alpha'} gamma_{beta,beta'}
0168   //                                  gamma_{gamma,gamma'} gamma_{delta,delta'}
0169   return EpsilonTensorContraction(BuildTensor(ten,v),i,j,k,l);
0170 }
0171
0172 template<typename Scalar>
0173 Scalar
0174 ATOOLS::EpsilonTensorContraction(const Lorentz_Ten4<Scalar>& ten, int i, int j, int k, int l) {
0175   Lorentz_Ten4<Scalar> t;
0176   // c = eps^{alpha,beta,gamma,delta} t^{alpha',beta',gamma',delta'}
0177   //                                  gamma_{alpha,alpha'} gamma_{beta,beta'}
0178   //                                  gamma_{gamma,gamma'} gamma_{delta,delta'}
0179   if      ((i == 1) && (j == 2) && (k == 3) && (l == 4))
0180     {t = ten;}
0181   else if ((i == 1) && (j == 2) && (k == 4) && (l == 3))
0182     {t = ten.Transpose(3,4);}
0183   else if ((i == 1) && (j == 3) && (k == 2) && (l == 4))
0184     {t = ten.Transpose(2,3);}
0185   else if ((i == 1) && (j == 3) && (k == 4) && (l == 2))
0186     {t = ten.Transpose(3,4); t = t.Transpose(2,3);}
0187   else if ((i == 1) && (j == 4) && (k == 2) && (l == 3))
0188     {t = ten.Transpose(2,3); t = t.Transpose(3,4);}
0189   else if ((i == 1) && (j == 4) && (k == 3) && (l == 2))
0190     {t = ten.Transpose(2,4);}
0191   else if ((i == 2) && (j == 1) && (k == 3) && (l == 4))
0192     {t = ten.Transpose(1,2);}
0193   else if ((i == 2) && (j == 1) && (k == 4) && (l == 3))
0194     {t = ten.Transpose(1,2); t = t.Transpose(3,4);}
0195   else if ((i == 2) && (j == 3) && (k == 1) && (l == 4))
0196     {t = ten.Transpose(2,3); t = t.Transpose(3,4);}
0197   else if ((i == 2) && (j == 3) && (k == 4) && (l == 1))
0198     {t = ten.Transpose(3,4); t = t.Transpose(2,3); t = t.Transpose(1,2);}
0199   else if ((i == 2) && (j == 4) && (k == 1) && (l == 3))
0200     {t = ten.Transpose(1,3); t = t.Transpose(2,3); t = t.Transpose(3,4);}
0201   else if ((i == 2) && (j == 4) && (k == 3) && (l == 1))
0202     {t = ten.Transpose(2,4); t = t.Transpose(1,2);}
0203   else if ((i == 3) && (j == 1) && (k == 2) && (l == 4))
0204     {t = ten.Transpose(1,2); t = t.Transpose(2,3);}
0205   else if ((i == 3) && (j == 1) && (k == 4) && (l == 2))
0206     {t = ten.Transpose(1,2); t = t.Transpose(2,4); t = t.Transpose(3,4);}
0207   else if ((i == 3) && (j == 2) && (k == 1) && (l == 4))
0208     {t = ten.Transpose(1,3);}
0209   else if ((i == 3) && (j == 2) && (k == 4) && (l == 1))
0210     {t = ten.Transpose(1,4); t = t.Transpose(3,4);}
0211   else if ((i == 3) && (j == 4) && (k == 1) && (l == 2))
0212     {t = ten.Transpose(1,3); t = t.Transpose(2,4);}
0213   else if ((i == 3) && (j == 4) && (k == 2) && (l == 1))
0214     {t = ten.Transpose(1,4); t = t.Transpose(2,3); t = t.Transpose(3,4);}
0215   else if ((i == 4) && (j == 1) && (k == 2) && (l == 3))
0216     {t = ten.Transpose(1,2); t = t.Transpose(2,3); t = t.Transpose(3,4);}
0217   else if ((i == 4) && (j == 1) && (k == 3) && (l == 2))
0218     {t = ten.Transpose(1,2); t = t.Transpose(2,4);}
0219   else if ((i == 4) && (j == 2) && (k == 1) && (l == 3))
0220     {t = ten.Transpose(1,3); t = t.Transpose(3,4);}
0221   else if ((i == 4) && (j == 2) && (k == 3) && (l == 1))
0222     {t = ten.Transpose(1,4);}
0223   else if ((i == 4) && (j == 3) && (k == 1) && (l == 2))
0224     {t = ten.Transpose(1,3); t = t.Transpose(2,4); t = t.Transpose(3,4);}
0225   else if ((i == 4) && (j == 3) && (k == 2) && (l == 1))
0226     {t = ten.Transpose(1,4); t = t.Transpose(2,3);}
0227   else
0228     return 0.;
0229   return   t.at(0,1,2,3)
0230           -t.at(0,1,3,2)
0231           -t.at(0,2,1,3)
0232           +t.at(0,2,3,1)
0233           +t.at(0,3,1,2)
0234           -t.at(0,3,2,1)
0235
0236           -t.at(1,0,2,3)
0237           +t.at(1,0,3,2)
0238           +t.at(1,2,0,3)
0239           -t.at(1,2,3,0)
0240           -t.at(1,3,0,2)
0241           +t.at(1,3,2,0)
0242
0243           +t.at(2,0,1,3)
0244           -t.at(2,0,3,1)
0245           -t.at(2,1,0,3)
0246           +t.at(2,1,3,0)
0247           +t.at(2,3,0,1)
0248           -t.at(2,3,1,0)
0249
0250           -t.at(3,0,1,2)
0251           +t.at(3,0,2,1)
0252           +t.at(3,1,0,2)
0253           -t.at(3,1,2,0)
0254           -t.at(3,2,0,1)
0255           +t.at(3,2,1,0);
0256 }
0257
0258 #endif