Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-04-19 09:09:51 

0001 #ifndef ATOOLS_Math_Tensor_Contractions_H
0002 #define ATOOLS_Math_Tensor_Contractions_H
0003 
0004 /********************************************************************************************
0005 ******                                                                                 ******
0006 ******       tensor with itself                                                        ******
0007 ******                                                                                 ******
0008 ********************************************************************************************/
0009 
0010 template<typename Scalar>
0011 Scalar
0012 ATOOLS::Contraction(const Lorentz_Ten2<Scalar>& t) {
0013   // c = t^{alpha,beta} g_{alpha,beta}
0014   return t.at(0,0)-t.at(1,1)-t.at(2,2)-t.at(3,3);
0015 }
0016 
0017 template<typename Scalar>
0018 Vec4<Scalar>
0019 ATOOLS::Contraction(const Lorentz_Ten3<Scalar>& ten, int i, int j) {
0020   Lorentz_Ten3<Scalar> t;
0021   // v^{mu} = t{mu,alpha,alpha'} g_{alpha,alpha'}
0022   if      (((i==2) && (j==3)) || ((j==2) && (i==3)))
0023     t = ten;
0024   // v^{mu} = t{alpha,mu,alpha'} g_{alpha,alpha'}
0025   else if (((i==1) && (j==3)) || ((j==1) && (i==3)))
0026     t = ten.Transpose(1,2);
0027   // v^{mu} = t{alpha,alpha',mu} g_{alpha,alpha'}
0028   else if (((i==1) && (j==2)) || ((j==1) && (i==2)))
0029     t = ten.Transpose(1,3);
0030   else
0031     return Vec4<Scalar>(0.,0.,0.,0.);
0032   return Vec4<Scalar>(t.at(0,0,0)-t.at(0,1,1)-t.at(0,2,2)-t.at(0,3,3),
0033                       t.at(1,0,0)-t.at(1,1,1)-t.at(1,2,2)-t.at(1,3,3),
0034                       t.at(2,0,0)-t.at(2,1,1)-t.at(2,2,2)-t.at(2,3,3),
0035                       t.at(3,0,0)-t.at(3,1,1)-t.at(3,2,2)-t.at(3,3,3));
0036 }
0037 
0038 template<typename Scalar>
0039 Lorentz_Ten2<Scalar>
0040 ATOOLS::Contraction(const Lorentz_Ten4<Scalar>& ten, int i, int j) {
0041   Lorentz_Ten4<Scalar> t;
0042   // T^{mu,nu} = t^{mu,nu,alpha,alpha'} g_{alpha,alpha'}
0043   if      (((i==3) && (j==4)) || ((i==4) && (j==3)))
0044      t = ten;
0045   // T^{mu,nu} = t^{mu,alpha,nu,alpha'} g_{alpha,alpha'}
0046   else if (((i==2) && (j==4)) || ((i==4) && (j==2)))
0047      t = ten.Transpose(2,3);
0048   // T^{mu,nu} = t^{alpha,mu,nu,alpha'} g_{alpha,alpha'}
0049   else if (((i==1) && (j==4)) || ((i==4) && (j==1)))
0050     {t = ten.Transpose(1,2); t = t.Transpose(2,3);}
0051   // T^{mu,nu} = t^{mu,alpha,alpha',nu} g_{alpha,alpha'}
0052   else if (((i==2) && (j==3)) || ((i==3) && (j==2)))
0053      t = ten.Transpose(2,4);
0054   // T^{mu,nu} = t^{alpha,mu,alpha',nu} g_{alpha,alpha'}
0055   else if (((i==1) && (j==3)) || ((i==3) && (j==1)))
0056     {t = ten.Transpose(1,2); t = t.Transpose(2,4);}
0057   // T^{mu,nu} = t^{alpha,alpha',mu,nu} g_{alpha,alpha'}
0058   else if (((i==1) && (j==2)) || ((i==2) && (j==1)))
0059     {t = ten.Transpose(1,3); t = t.Transpose(2,4);}
0060   else
0061     return Lorentz_Ten2<Scalar>();
0062   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),
0063                               t.at(1,0,0,0)-t.at(1,0,1,1)-t.at(1,0,2,2)-t.at(1,0,3,3),
0064                               t.at(2,0,0,0)-t.at(2,0,1,1)-t.at(2,0,2,2)-t.at(2,0,3,3),
0065                               t.at(3,0,0,0)-t.at(3,0,1,1)-t.at(3,0,2,2)-t.at(3,0,3,3),
0066                               t.at(0,1,0,0)-t.at(0,1,1,1)-t.at(0,1,2,2)-t.at(0,1,3,3),
0067                               t.at(1,1,0,0)-t.at(1,1,1,1)-t.at(1,1,2,2)-t.at(1,1,3,3),
0068                               t.at(2,1,0,0)-t.at(2,1,1,1)-t.at(2,1,2,2)-t.at(2,1,3,3),
0069                               t.at(3,1,0,0)-t.at(3,1,1,1)-t.at(3,1,2,2)-t.at(3,1,3,3),
0070                               t.at(0,2,0,0)-t.at(0,2,1,1)-t.at(0,2,2,2)-t.at(0,2,3,3),
0071                               t.at(1,2,0,0)-t.at(1,2,1,1)-t.at(1,2,2,2)-t.at(1,2,3,3),
0072                               t.at(2,2,0,0)-t.at(2,2,1,1)-t.at(2,2,2,2)-t.at(2,2,3,3),
0073                               t.at(3,2,0,0)-t.at(3,2,1,1)-t.at(3,2,2,2)-t.at(3,2,3,3),
0074                               t.at(0,3,0,0)-t.at(0,3,1,1)-t.at(0,3,2,2)-t.at(0,3,3,3),
0075                               t.at(1,3,0,0)-t.at(1,3,1,1)-t.at(1,3,2,2)-t.at(1,3,3,3),
0076                               t.at(2,3,0,0)-t.at(2,3,1,1)-t.at(2,3,2,2)-t.at(2,3,3,3),
0077                               t.at(3,3,0,0)-t.at(3,3,1,1)-t.at(3,3,2,2)-t.at(3,3,3,3));
0078 }
0079 
0080 template<typename Scalar>
0081 Scalar
0082 ATOOLS::Contraction(const Lorentz_Ten4<Scalar>& ten, int i, int j, int k, int l) {
0083   Lorentz_Ten4<Scalar> t;
0084   // c = t^{alpha,alpha',beta,beta'} g_{alpha,alpha'} g_{beta,beta'}
0085   if      (((((i==1) && (j==2)) || ((i==2) && (j==1)))  && 
0086             (((k==3) && (l==4)) || ((k==4) && (l==3)))) ||
0087            ((((i==3) && (j==4)) || ((i==4) && (j==3)))  && 
0088             (((k==1) && (l==2)) || ((k==2) && (l==1)))))
0089     t = ten;
0090   // c = t^{alpha,beta,alpha',beta'} g_{alpha,alpha'} g_{beta,beta'}
0091   else if (((((i==1) && (j==3)) || ((i==3) && (j==1)))  && 
0092             (((k==2) && (l==4)) || ((k==4) && (l==2)))) ||
0093            ((((i==2) && (j==4)) || ((i==4) && (j==2)))  && 
0094             (((k==1) && (l==3)) || ((k==3) && (l==1)))))
0095     t = ten.Transpose(2,3);
0096   // c = t^{alpha,beta,beta',alpha'} g_{alpha,alpha'} g_{beta,beta'}
0097   else if (((((i==1) && (j==4)) || ((i==4) && (j==1)))  && 
0098             (((k==2) && (l==3)) || ((k==3) && (l==2)))) ||
0099            ((((i==2) && (j==3)) || ((i==3) && (j==2)))  && 
0100             (((k==1) && (l==4)) || ((k==4) && (l==1)))))
0101     t = ten.Transpose(2,4);
0102   else
0103     return Scalar(0.);
0104   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)
0105                 -t.at(1,1,0,0)+t.at(1,1,1,1)+t.at(1,1,2,2)+t.at(1,1,3,3)
0106                 -t.at(2,2,0,0)+t.at(2,2,1,1)+t.at(2,2,2,2)+t.at(2,2,3,3)
0107                 -t.at(3,3,0,0)+t.at(3,3,1,1)+t.at(3,3,2,2)+t.at(3,3,3,3));
0108 }
0109 
0110 /********************************************************************************************
0111 ******                                                                                 ******
0112 ******       tensor with vector                                                        ******
0113 ******                                                                                 ******
0114 ********************************************************************************************/
0115 
0116 template<typename Scal1, typename Scal2> 
0117 Vec4<PROMOTE(Scal1,Scal2)> 
0118 ATOOLS::Contraction(const Lorentz_Ten2<Scal1>& ten, int i, const Vec4<Scal2>& p) {
0119   Lorentz_Ten2<Scal1> t;
0120   if ((i==1) || (i==2)) {
0121     // v^{mu} = t^{alpha,mu} p^{beta} g_{alpha,beta}
0122     if      (i==1)
0123       t = ten.Transpose();
0124     // v^{mu} = t^{mu,alpha} p^{beta} g_{alpha,beta}
0125     else if (i==2)
0126       t = ten;
0127     return Vec4<PROMOTE(Scal1,Scal2)>(
0128                         t.at(0,0)*p[0]-t.at(0,1)*p[1]-t.at(0,2)*p[2]-t.at(0,3)*p[3],
0129                         t.at(1,0)*p[0]-t.at(1,1)*p[1]-t.at(1,2)*p[2]-t.at(1,3)*p[3],
0130                         t.at(2,0)*p[0]-t.at(2,1)*p[1]-t.at(2,2)*p[2]-t.at(2,3)*p[3],
0131                         t.at(3,0)*p[0]-t.at(3,1)*p[1]-t.at(3,2)*p[2]-t.at(3,3)*p[3]);
0132   }
0133   return Vec4<PROMOTE(Scal1,Scal2)>();
0134 }
0135 
0136 template<typename Scal1, typename Scal2, typename Scal3> 
0137 PROMOTE(PROMOTE(Scal1,Scal2),Scal3)
0138 ATOOLS::Contraction(const Lorentz_Ten2<Scal1>& t,
0139                     const Vec4<Scal2>& p, const Vec4<Scal3>& q) {
0140   // c = t^{alpha,beta} p^{alpha'} q^{beta'} g_{alpha',alpha} g_{beta,beta'}
0141   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])
0142          -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])
0143          -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])
0144          -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]);
0145 }
0146 
0147 template<typename Scal1, typename Scal2>
0148 Lorentz_Ten2<PROMOTE(Scal1,Scal2)>
0149 ATOOLS::Contraction(const Lorentz_Ten3<Scal1>& ten, int i, const Vec4<Scal2>& v) {
0150   Lorentz_Ten3<Scal1> t;
0151   // c^{mu,nu} = t^{alpha,mu,nu} g_{alpha,alpha'} v^{alpha'}
0152   if      (i==1)
0153     {t = ten.Transpose(1,2); t = t.Transpose(2,3);}
0154   // c^{mu,nu} = t^{mu,alpha,nu} g_{alpha,alpha'} v^{alpha'}
0155   else if (i==2)
0156      t = ten.Transpose(2,3);
0157   // c^{mu,nu} = t^{mu,nu,alpha} g_{alpha,alpha'} v^{alpha'}
0158   else if (i==3)
0159      t = ten;
0160   else
0161     return Lorentz_Ten2<PROMOTE(Scal1,Scal2)>();
0162   return Lorentz_Ten2<PROMOTE(Scal1,Scal2)>(
0163             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],
0164             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],
0165             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],
0166             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],
0167             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],
0168             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],
0169             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],
0170             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],
0171             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],
0172             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],
0173             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],
0174             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],
0175             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],
0176             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],
0177             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],
0178             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]);
0179 }
0180 
0181 template<typename Scal1, typename Scal2>
0182 Lorentz_Ten3<PROMOTE(Scal1,Scal2)>
0183 ATOOLS::Contraction(const Lorentz_Ten4<Scal1>& ten, int i, const Vec4<Scal2>& v) {
0184   Lorentz_Ten4<Scal1> t;
0185   // c^{mu,nu,rho} = t^{alpha,mu,nu,rho} g_{alpha,alpha'} v^{alpha'}
0186   if      (i==1)
0187     {t = ten.Transpose(1,2); t = t.Transpose(2,3); t = t.Transpose(3,4);}
0188   // c^{mu,nu,rho} = t^{mu,alpha,nu,rho} g_{alpha,alpha'} v^{alpha'}
0189   else if (i==2)
0190     {t = ten.Transpose(2,3); t = t.Transpose(3,4);}
0191   // c^{mu,nu,rho} = t^{mu,nu,alpha,rho} g_{alpha,alpha'} v^{alpha'}
0192   else if (i==3)
0193      t = ten.Transpose(3,4);
0194   // c^{mu,nu,rho} = t^{mu,nu,rho,alpha} g_{alpha,alpha'} v^{alpha'}
0195   else if (i==4)
0196      t = ten;
0197   else
0198     return Lorentz_Ten3<PROMOTE(Scal1,Scal2)>();
0199   return Lorentz_Ten3<PROMOTE(Scal1,Scal2)>(
0200     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],
0201     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],
0202     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],
0203     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],
0204     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],
0205     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],
0206     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],
0207     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],
0208     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],
0209     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],
0210     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],
0211     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],
0212     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],
0213     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],
0214     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],
0215     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],
0216 
0217     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],
0218     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],
0219     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],
0220     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],
0221     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],
0222     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],
0223     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],
0224     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],
0225     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],
0226     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],
0227     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],
0228     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],
0229     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],
0230     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],
0231     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],
0232     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],
0233 
0234     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],
0235     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],
0236     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],
0237     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],
0238     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],
0239     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],
0240     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],
0241     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],
0242     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],
0243     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],
0244     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],
0245     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],
0246     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],
0247     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],
0248     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],
0249     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],
0250 
0251     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],
0252     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],
0253     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],
0254     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],
0255     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],
0256     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],
0257     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],
0258     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],
0259     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],
0260     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],
0261     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],
0262     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],
0263     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],
0264     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],
0265     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],
0266     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]);
0267 }
0268 
0269 /********************************************************************************************
0270 ******                                                                                 ******
0271 ******       tensor with tensor                                                        ******
0272 ******                                                                                 ******
0273 ********************************************************************************************/
0274 
0275 template<typename Scal1, typename Scal2> 
0276 PROMOTE(Scal1,Scal2) 
0277 ATOOLS::Contraction(const Lorentz_Ten2<Scal1>& ten, int i, int j,
0278                     const Lorentz_Ten2<Scal2>& sen, int k, int l) {
0279   Lorentz_Ten2<Scal1> t; Lorentz_Ten2<Scal2> s;
0280   // c = t^{alpha,beta} s^{alpha',beta'} g_{alpha',alpha} g_{beta,beta'}
0281   if      (((i==1) && (j==2) && (k==1) && (l==2)) ||
0282            ((i==2) && (j==1) && (k==2) && (l==1)))
0283     {t = ten.Transpose(); s = sen;}
0284   // c = t^{alpha,beta} s^{beta',alpha'} g_{alpha',alpha} g_{beta,beta'}
0285   else if (((i==1) && (j==2) && (k==2) && (l==1)) ||
0286            ((i==2) && (j==1) && (k==1) && (l==2)))
0287     {t = ten;             s = sen;}
0288   else 
0289     return PROMOTE(Scal1,Scal2)();
0290   return PROMOTE(Scal1,Scal2)
0291     (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))
0292    -(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))
0293    -(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))
0294    -(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));
0295 }
0296 
0297 template<typename Scal1, typename Scal2> 
0298 Lorentz_Ten2<PROMOTE(Scal1,Scal2)> 
0299 ATOOLS::Contraction(const Lorentz_Ten2<Scal1>& ten,int i,
0300                     const Lorentz_Ten2<Scal2>& sen, int j) {
0301   Lorentz_Ten2<Scal1> t; Lorentz_Ten2<Scal2> s;
0302   if (((i==1) || (i==2)) && ((j==1) || (j==2))) {
0303     // T^{mu,nu} = t^{alpha,mu} s^{beta,nu} g_{alpha,beta}
0304     if      ((i==1) && (j==1))
0305       {t = ten.Transpose(); s = sen;}
0306     // T^{mu,nu} = t^{alpha,mu} s^{nu,beta} g_{alpha,beta}
0307     else if ((i==1) && (j==2))
0308       {t = ten.Transpose(); s = sen.Transpose();}
0309     // T^{mu,nu} = t^{mu,alpha} s^{beta,nu} g_{alpha,beta}
0310     else if ((i==2) && (j==1))
0311       {t = ten; s = sen;}
0312     // T^{mu,nu} = t^{mu,alpha} s^{nu,beta} g_{alpha,beta}
0313     else if ((i==2) && (j==2))
0314       {t = ten; s = sen.Transpose();}
0315     return Lorentz_Ten2<PROMOTE(Scal1,Scal2)>(
0316       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),
0317       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),
0318       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),
0319       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),
0320       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),
0321       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),
0322       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),
0323       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),
0324       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),
0325       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),
0326       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),
0327       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),
0328       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),
0329       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),
0330       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),
0331       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));
0332   }
0333   return Lorentz_Ten2<PROMOTE(Scal1,Scal2)>();
0334 }
0335 
0336 template<typename Scal1, typename Scal2> 
0337 Vec4<PROMOTE(Scal1,Scal2)> 
0338 ATOOLS::Contraction(const Lorentz_Ten3<Scal1>& ten, int i, int j,
0339                     const Lorentz_Ten2<Scal2>& sen, int k, int l) {
0340   Lorentz_Ten3<Scal1> t; Lorentz_Ten2<Scal2> s;
0341   // v^{mu} = t^{mu,alpha,beta} s^{alpha',beta'} g_{alpha,alpha'} g_{beta,beta'}
0342   if      (((i==2) && (j==3) && (k==1) && (l==2)) || ((i==3) && (j==2) && (k==2) && (l==1)))
0343   // v^{mu} = t^{alpha,mu,beta} s^{alpha',beta'} g_{alpha,alpha'} g_{beta,beta'}
0344     {t = ten; s = sen;}
0345   else if (((i==1) && (j==3) && (k==1) && (l==2)) || ((i==3) && (j==1) && (k==2) && (l==1)))
0346     {t = ten.Transpose(1,2); s = sen;}
0347   // v^{mu} = t^{alpha,beta,mu} s^{alpha',beta'} g_{alpha,alpha'} g_{beta,beta'}
0348   else if (((i==1) && (j==2) && (k==1) && (l==2)) || ((i==2) && (j==1) && (k==2) && (l==1)))
0349     {t = ten.Transpose(2,3); t = t.Transpose(1,2); s = sen;}
0350   // v^{mu} = t^{mu,alpha,beta} s^{beta',alpha'} g_{alpha,alpha'} g_{beta,beta'}
0351   else if (((i==2) && (j==3) && (k==2) && (l==1)) || ((i==3) && (j==2) && (k==1) && (l==2)))
0352     {t = ten; s = sen.Transpose();}
0353   // v^{mu} = t^{alpha,mu,beta} s^{beta',alpha'} g_{alpha,alpha'} g_{beta,beta'}
0354   else if (((i==1) && (j==3) && (k==2) && (l==1)) || ((i==3) && (j==1) && (k==1) && (l==2)))
0355     {t = ten.Transpose(1,2); s = sen.Transpose();}
0356   // v^{mu} = t^{alpha,beta,mu} s^{beta',alpha'} g_{alpha,alpha'} g_{beta,beta'}
0357   else if (((i==1) && (j==2) && (k==2) && (l==1)) || ((i==2) && (j==1) && (k==1) && (l==2)))
0358     {t = ten.Transpose(2,3); t = t.Transpose(1,2); s = sen.Transpose();}
0359   else
0360     return Vec4<PROMOTE(Scal1,Scal2)>(0.,0.,0.,0.);
0361   return Vec4<PROMOTE(Scal1,Scal2)>(
0362     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)
0363    -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)
0364    -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)
0365    -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),
0366     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)
0367    -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)
0368    -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)
0369    -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),
0370     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)
0371    -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)
0372    -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)
0373    -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),
0374     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)
0375    -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)
0376    -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)
0377    -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));
0378 }
0379 
0380 template<typename Scal1, typename Scal2> 
0381 Lorentz_Ten3<PROMOTE(Scal1,Scal2)> 
0382 ATOOLS::Contraction(const Lorentz_Ten3<Scal1>& ten, int i,
0383                     const Lorentz_Ten2<Scal2>& sen, int j) {
0384   Lorentz_Ten3<Scal1> t; Lorentz_Ten2<Scal2> s;
0385   // T^{mu,nu,rho} = t^{mu,nu,alpha} s^{rho,alpha'} g_{alpha,alpha'}
0386   if      ((i==3) && (j==2))
0387     {t = ten; s = sen;}
0388   else if ((i==2) && (j==2))
0389     {t = ten.Transpose(2,3); s = sen;}
0390   else if ((i==1) && (j==2))
0391     {t = ten.Transpose(1,2); t = t.Transpose(2,3); s = sen;}
0392   else if ((i==3) && (j==1))
0393     {t = ten; s = sen.Transpose();}
0394   else if ((i==2) && (j==1))
0395     {t = ten.Transpose(2,3); s = sen.Transpose();}
0396   else if ((i==1) && (j==1))
0397     {t = ten.Transpose(1,2); t = t.Transpose(2,3); s = sen.Transpose();}
0398   else
0399     return Lorentz_Ten3<PROMOTE(Scal1,Scal2)>();
0400   return Lorentz_Ten3<PROMOTE(Scal1,Scal2)>( 
0401               t.at(0,0,0)*s.at(0,0)-t.at(0,0,1)*s.at(0,1)
0402                             -t.at(0,0,2)*s.at(0,2)-t.at(0,0,3)*s.at(0,3),
0403               t.at(1,0,0)*s.at(0,0)-t.at(1,0,1)*s.at(0,1)
0404                             -t.at(1,0,2)*s.at(0,2)-t.at(1,0,3)*s.at(0,3),
0405               t.at(2,0,0)*s.at(0,0)-t.at(2,0,1)*s.at(0,1)
0406                             -t.at(2,0,2)*s.at(0,2)-t.at(2,0,3)*s.at(0,3),
0407               t.at(3,0,0)*s.at(0,0)-t.at(3,0,1)*s.at(0,1)
0408                             -t.at(3,0,2)*s.at(0,2)-t.at(3,0,3)*s.at(0,3),
0409               t.at(0,1,0)*s.at(0,0)-t.at(0,1,1)*s.at(0,1)
0410                             -t.at(0,1,2)*s.at(0,2)-t.at(0,1,3)*s.at(0,3),
0411               t.at(1,1,0)*s.at(0,0)-t.at(1,1,1)*s.at(0,1)
0412                             -t.at(1,1,2)*s.at(0,2)-t.at(1,1,3)*s.at(0,3),
0413               t.at(2,1,0)*s.at(0,0)-t.at(2,1,1)*s.at(0,1)
0414                             -t.at(2,1,2)*s.at(0,2)-t.at(2,1,3)*s.at(0,3),
0415               t.at(3,1,0)*s.at(0,0)-t.at(3,1,1)*s.at(0,1)
0416                             -t.at(3,1,2)*s.at(0,2)-t.at(3,1,3)*s.at(0,3),
0417               t.at(0,2,0)*s.at(0,0)-t.at(0,2,1)*s.at(0,1)
0418                             -t.at(0,2,2)*s.at(0,2)-t.at(0,2,3)*s.at(0,3),
0419               t.at(1,2,0)*s.at(0,0)-t.at(1,2,1)*s.at(0,1)
0420                             -t.at(1,2,2)*s.at(0,2)-t.at(1,2,3)*s.at(0,3),
0421               t.at(2,2,0)*s.at(0,0)-t.at(2,2,1)*s.at(0,1)
0422                             -t.at(2,2,2)*s.at(0,2)-t.at(2,2,3)*s.at(0,3),
0423               t.at(3,2,0)*s.at(0,0)-t.at(3,2,1)*s.at(0,1)
0424                             -t.at(3,2,2)*s.at(0,2)-t.at(3,2,3)*s.at(0,3),
0425               t.at(0,3,0)*s.at(0,0)-t.at(0,3,1)*s.at(0,1)
0426                             -t.at(0,3,2)*s.at(0,2)-t.at(0,3,3)*s.at(0,3),
0427               t.at(1,3,0)*s.at(0,0)-t.at(1,3,1)*s.at(0,1)
0428                             -t.at(1,3,2)*s.at(0,2)-t.at(1,3,3)*s.at(0,3),
0429               t.at(2,3,0)*s.at(0,0)-t.at(2,3,1)*s.at(0,1)
0430                             -t.at(2,3,2)*s.at(0,2)-t.at(2,3,3)*s.at(0,3),
0431               t.at(3,3,0)*s.at(0,0)-t.at(3,3,1)*s.at(0,1)
0432                             -t.at(3,3,2)*s.at(0,2)-t.at(3,3,3)*s.at(0,3),
0433 
0434               t.at(0,0,0)*s.at(1,0)-t.at(0,0,1)*s.at(1,1)
0435                             -t.at(0,0,2)*s.at(1,2)-t.at(0,0,3)*s.at(1,3),
0436               t.at(1,0,0)*s.at(1,0)-t.at(1,0,1)*s.at(1,1)
0437                             -t.at(1,0,2)*s.at(1,2)-t.at(1,0,3)*s.at(1,3),
0438               t.at(2,0,0)*s.at(1,0)-t.at(2,0,1)*s.at(1,1)
0439                             -t.at(2,0,2)*s.at(1,2)-t.at(2,0,3)*s.at(1,3),
0440               t.at(3,0,0)*s.at(1,0)-t.at(3,0,1)*s.at(1,1)
0441                             -t.at(3,0,2)*s.at(1,2)-t.at(3,0,3)*s.at(1,3),
0442               t.at(0,1,0)*s.at(1,0)-t.at(0,1,1)*s.at(1,1)
0443                             -t.at(0,1,2)*s.at(1,2)-t.at(0,1,3)*s.at(1,3),
0444               t.at(1,1,0)*s.at(1,0)-t.at(1,1,1)*s.at(1,1)
0445                             -t.at(1,1,2)*s.at(1,2)-t.at(1,1,3)*s.at(1,3),
0446               t.at(2,1,0)*s.at(1,0)-t.at(2,1,1)*s.at(1,1)
0447                             -t.at(2,1,2)*s.at(1,2)-t.at(2,1,3)*s.at(1,3),
0448               t.at(3,1,0)*s.at(1,0)-t.at(3,1,1)*s.at(1,1)
0449                             -t.at(3,1,2)*s.at(1,2)-t.at(3,1,3)*s.at(1,3),
0450               t.at(0,2,0)*s.at(1,0)-t.at(0,2,1)*s.at(1,1)
0451                             -t.at(0,2,2)*s.at(1,2)-t.at(0,2,3)*s.at(1,3),
0452               t.at(1,2,0)*s.at(1,0)-t.at(1,2,1)*s.at(1,1)
0453                             -t.at(1,2,2)*s.at(1,2)-t.at(1,2,3)*s.at(1,3),
0454               t.at(2,2,0)*s.at(1,0)-t.at(2,2,1)*s.at(1,1)
0455                             -t.at(2,2,2)*s.at(1,2)-t.at(2,2,3)*s.at(1,3),
0456               t.at(3,2,0)*s.at(1,0)-t.at(3,2,1)*s.at(1,1)
0457                             -t.at(3,2,2)*s.at(1,2)-t.at(3,2,3)*s.at(1,3),
0458               t.at(0,3,0)*s.at(1,0)-t.at(0,3,1)*s.at(1,1)
0459                             -t.at(0,3,2)*s.at(1,2)-t.at(0,3,3)*s.at(1,3),
0460               t.at(1,3,0)*s.at(1,0)-t.at(1,3,1)*s.at(1,1)
0461                             -t.at(1,3,2)*s.at(1,2)-t.at(1,3,3)*s.at(1,3),
0462               t.at(2,3,0)*s.at(1,0)-t.at(2,3,1)*s.at(1,1)
0463                             -t.at(2,3,2)*s.at(1,2)-t.at(2,3,3)*s.at(1,3),
0464               t.at(3,3,0)*s.at(1,0)-t.at(3,3,1)*s.at(1,1)
0465                             -t.at(3,3,2)*s.at(1,2)-t.at(3,3,3)*s.at(1,3),
0466 
0467               t.at(0,0,0)*s.at(2,0)-t.at(0,0,1)*s.at(2,1)
0468                             -t.at(0,0,2)*s.at(2,2)-t.at(0,0,3)*s.at(2,3),
0469               t.at(1,0,0)*s.at(2,0)-t.at(1,0,1)*s.at(2,1)
0470                             -t.at(1,0,2)*s.at(2,2)-t.at(1,0,3)*s.at(2,3),
0471               t.at(2,0,0)*s.at(2,0)-t.at(2,0,1)*s.at(2,1)
0472                             -t.at(2,0,2)*s.at(2,2)-t.at(2,0,3)*s.at(2,3),
0473               t.at(3,0,0)*s.at(2,0)-t.at(3,0,1)*s.at(2,1)
0474                             -t.at(3,0,2)*s.at(2,2)-t.at(3,0,3)*s.at(2,3),
0475               t.at(0,1,0)*s.at(2,0)-t.at(0,1,1)*s.at(2,1)
0476                             -t.at(0,1,2)*s.at(2,2)-t.at(0,1,3)*s.at(2,3),
0477               t.at(1,1,0)*s.at(2,0)-t.at(1,1,1)*s.at(2,1)
0478                             -t.at(1,1,2)*s.at(2,2)-t.at(1,1,3)*s.at(2,3),
0479               t.at(2,1,0)*s.at(2,0)-t.at(2,1,1)*s.at(2,1)
0480                             -t.at(2,1,2)*s.at(2,2)-t.at(2,1,3)*s.at(2,3),
0481               t.at(3,1,0)*s.at(2,0)-t.at(3,1,1)*s.at(2,1)
0482                             -t.at(3,1,2)*s.at(2,2)-t.at(3,1,3)*s.at(2,3),
0483               t.at(0,2,0)*s.at(2,0)-t.at(0,2,1)*s.at(2,1)
0484                             -t.at(0,2,2)*s.at(2,2)-t.at(0,2,3)*s.at(2,3),
0485               t.at(1,2,0)*s.at(2,0)-t.at(1,2,1)*s.at(2,1)
0486                             -t.at(1,2,2)*s.at(2,2)-t.at(1,2,3)*s.at(2,3),
0487               t.at(2,2,0)*s.at(2,0)-t.at(2,2,1)*s.at(2,1)
0488                             -t.at(2,2,2)*s.at(2,2)-t.at(2,2,3)*s.at(2,3),
0489               t.at(3,2,0)*s.at(2,0)-t.at(3,2,1)*s.at(2,1)
0490                             -t.at(3,2,2)*s.at(2,2)-t.at(3,2,3)*s.at(2,3),
0491               t.at(0,3,0)*s.at(2,0)-t.at(0,3,1)*s.at(2,1)
0492                             -t.at(0,3,2)*s.at(2,2)-t.at(0,3,3)*s.at(2,3),
0493               t.at(1,3,0)*s.at(2,0)-t.at(1,3,1)*s.at(2,1)
0494                             -t.at(1,3,2)*s.at(2,2)-t.at(1,3,3)*s.at(2,3),
0495               t.at(2,3,0)*s.at(2,0)-t.at(2,3,1)*s.at(2,1)
0496                             -t.at(2,3,2)*s.at(2,2)-t.at(2,3,3)*s.at(2,3),
0497               t.at(3,3,0)*s.at(2,0)-t.at(3,3,1)*s.at(2,1)
0498                             -t.at(3,3,2)*s.at(2,2)-t.at(3,3,3)*s.at(2,3),
0499 
0500               t.at(0,0,0)*s.at(3,0)-t.at(0,0,1)*s.at(3,1)
0501                             -t.at(0,0,2)*s.at(3,2)-t.at(0,0,3)*s.at(3,3),
0502               t.at(1,0,0)*s.at(3,0)-t.at(1,0,1)*s.at(3,1)
0503                             -t.at(1,0,2)*s.at(3,2)-t.at(1,0,3)*s.at(3,3),
0504               t.at(2,0,0)*s.at(3,0)-t.at(2,0,1)*s.at(3,1)
0505                             -t.at(2,0,2)*s.at(3,2)-t.at(2,0,3)*s.at(3,3),
0506               t.at(3,0,0)*s.at(3,0)-t.at(3,0,1)*s.at(3,1)
0507                             -t.at(3,0,2)*s.at(3,2)-t.at(3,0,3)*s.at(3,3),
0508               t.at(0,1,0)*s.at(3,0)-t.at(0,1,1)*s.at(3,1)
0509                             -t.at(0,1,2)*s.at(3,2)-t.at(0,1,3)*s.at(3,3),
0510               t.at(1,1,0)*s.at(3,0)-t.at(1,1,1)*s.at(3,1)
0511                             -t.at(1,1,2)*s.at(3,2)-t.at(1,1,3)*s.at(3,3),
0512               t.at(2,1,0)*s.at(3,0)-t.at(2,1,1)*s.at(3,1)
0513                             -t.at(2,1,2)*s.at(3,2)-t.at(2,1,3)*s.at(3,3),
0514               t.at(3,1,0)*s.at(3,0)-t.at(3,1,1)*s.at(3,1)
0515                             -t.at(3,1,2)*s.at(3,2)-t.at(3,1,3)*s.at(3,3),
0516               t.at(0,2,0)*s.at(3,0)-t.at(0,2,1)*s.at(3,1)
0517                             -t.at(0,2,2)*s.at(3,2)-t.at(0,2,3)*s.at(3,3),
0518               t.at(1,2,0)*s.at(3,0)-t.at(1,2,1)*s.at(3,1)
0519                             -t.at(1,2,2)*s.at(3,2)-t.at(1,2,3)*s.at(3,3),
0520               t.at(2,2,0)*s.at(3,0)-t.at(2,2,1)*s.at(3,1)
0521                             -t.at(2,2,2)*s.at(3,2)-t.at(2,2,3)*s.at(3,3),
0522               t.at(3,2,0)*s.at(3,0)-t.at(3,2,1)*s.at(3,1)
0523                             -t.at(3,2,2)*s.at(3,2)-t.at(3,2,3)*s.at(3,3),
0524               t.at(0,3,0)*s.at(3,0)-t.at(0,3,1)*s.at(3,1)
0525                             -t.at(0,3,2)*s.at(3,2)-t.at(0,3,3)*s.at(3,3),
0526               t.at(1,3,0)*s.at(3,0)-t.at(1,3,1)*s.at(3,1)
0527                             -t.at(1,3,2)*s.at(3,2)-t.at(1,3,3)*s.at(3,3),
0528               t.at(2,3,0)*s.at(3,0)-t.at(2,3,1)*s.at(3,1)
0529                             -t.at(2,3,2)*s.at(3,2)-t.at(2,3,3)*s.at(3,3),
0530               t.at(3,3,0)*s.at(3,0)-t.at(3,3,1)*s.at(3,1)
0531                             -t.at(3,3,2)*s.at(3,2)-t.at(3,3,3)*s.at(3,3));
0532 
0533 }
0534 
0535 template<typename Scal1, typename Scal2>
0536 Vec4<PROMOTE(Scal1,Scal2)>
0537 ATOOLS::Contraction(const Lorentz_Ten4<Scal1>& ten, int i, int j, int k,
0538                     const Lorentz_Ten3<Scal2>& sen, int l, int m, int n) {
0539   // v^{mu} = t^{mu,alpha,beta,gamma} s^{alpha',beta',gamma'}
0540   //                              g_{alpha,alpha'} g_{beta,beta'} g_{gamma,gamma'}
0541   if ((i!=2) && (j!=3) && (k!=4) && (l!=2) && (m!=3) && (n!=4)) {
0542     THROW(fatal_error,"tensor contraction not specified for these indizes ...");
0543   }
0544   Vec4<PROMOTE(Scal1,Scal2)> end;
0545   for (unsigned short int I=0; I<4; ++I) {
0546     PROMOTE(Scal1,Scal2) res(ten.at(0,0,0,0)-ten.at(0,0,0,0));
0547     for (unsigned short int a=0; a<4; ++a)
0548       for (unsigned short int b=0; b<4; ++b)
0549         for (unsigned short int c=0; c<4; ++c) {
0550           double ga(-1.),gb(-1.),gc(-1.);
0551           if (a==0) ga=1.; if (b==0) gb=1.; if (c==0) gc=1.;
0552           res += ga*gb*gc*ten.at(I,a,b,c)*sen.at(a,b,c);
0553         }
0554     end[I] = res;
0555   }
0556   return end;
0557 }
0558 
0559 template<typename Scal1, typename Scal2>
0560 Lorentz_Ten3<PROMOTE(Scal1,Scal2)>
0561 ATOOLS::Contraction(const Lorentz_Ten4<Scal1>& ten, int i, int j,
0562                     const Lorentz_Ten3<Scal2>& sen, int k, int l) {
0563   // T^{mu,nu,rho} = t^{mu,nu,alpha,beta} s^{rho,alpha',beta'}
0564   //                                g_{alpha,alpha'} g_{beta,beta'}
0565   if ((i!=3) && (j!=4) && (k!=3) && (l!=4)) {
0566     THROW(fatal_error,"tensor contraction not specified for these indizes ...");
0567   }
0568   Vec4<PROMOTE(Scal1,Scal2)> end;
0569   for (unsigned short int I=0; I<4; ++I)
0570     for (unsigned short int J=0; J<4; ++J)
0571       for (unsigned short int K=0; K<4; ++K) {
0572         PROMOTE(Scal1,Scal2) res(ten.at(0,0,0,0)-ten.at(0,0,0,0));
0573         for (unsigned short int a=0; a<4; ++a)
0574           for (unsigned short int b=0; b<4; ++b) {
0575             double ga(-1.),gb(-1.);
0576             if (a==0) ga=1.; if (b==0) gb=1.;
0577             res += ga*gb*ten.at(I,J,a,b)*sen.at(K,a,b);
0578           }
0579         end.at(I,J,K) = res;
0580       }
0581   return end;
0582 }
0583 
0584 template<typename Scal1, typename Scal2>
0585 PROMOTE(Scal1,Scal2)
0586 ATOOLS::Contraction(const Lorentz_Ten4<Scal1>& ten, int i, int j, int k, int l,
0587                     const Lorentz_Ten4<Scal2>& sen, int m, int n, int o, int p) {
0588   Lorentz_Ten4<Scal1> t; Lorentz_Ten4<Scal2> s;
0589   // c = t^{alpha,beta,gamma,delta} s^{alpha',beta',gamma',delta'}
0590   //            g_{alpha,alpha'} g_{beta,beta'} g_{gamma.gamma'} g_{delta,delta'}
0591   if       ((i==m) && (j==n) && (k==o) && (l==p))
0592     {t=ten; s=sen;}
0593   else if  ((i==m) && (j==n) && (k==p) && (l==o))
0594     {t=ten; s=sen.Transpose(3,4);}
0595   else if  ((i==m) && (j==o) && (k==n) && (l==p))
0596     {t=ten; s=sen.Transpose(2,3);}
0597   else if  ((i==m) && (j==o) && (k==p) && (l==n))
0598     {t=ten; s=sen.Transpose(3,4); s=s.Transpose(2,3);}
0599   else if  ((i==n) && (j==m) && (k==o) && (l==p))
0600     {t=ten; s=sen.Transpose(1,2);}
0601   else if  ((i==n) && (j==m) && (k==p) && (l==o))
0602     {t=ten; s=sen.Transpose(1,2); s=s.Transpose(3,4);}
0603   else if  ((i==n) && (j==o) && (k==m) && (l==p))
0604     {t=ten; s=sen.Transpose(1,3); s=s.Transpose(2,3);}
0605   else if  ((i==n) && (j==o) && (k==p) && (l==m))
0606     {t=ten; s=sen.Transpose(3,4); s=s.Transpose(2,3); s=s.Transpose(1,2);}
0607   else if  ((i==n) && (j==p) && (k==m) && (l==o))
0608     {t=ten; s=sen.Transpose(1,3); s=s.Transpose(2,4); s=s.Transpose(2,3);}
0609   else if  ((i==n) && (j==p) && (k==o) && (l==m))
0610     {t=ten; s=sen.Transpose(2,4); s=s.Transpose(1,2);}
0611   else if  ((i==o) && (j==m) && (k==n) && (l==p))
0612     {t=ten; s=sen.Transpose(1,2); s=s.Transpose(2,3);}
0613   else if  ((i==o) && (j==m) && (k==p) && (l==n))
0614     {t=ten; s=sen.Transpose(1,2); s=s.Transpose(2,3); s=s.Transpose(2,4);}
0615   else if  ((i==o) && (j==n) && (k==m) && (l==p))
0616     {t=ten; s=sen.Transpose(1,3);}
0617   else if  ((i==o) && (j==n) && (k==p) && (l==m))
0618     {t=ten; s=sen.Transpose(1,4); s=s.Transpose(3,4);}
0619   else if  ((i==o) && (j==p) && (k==m) && (l==n))
0620     {t=ten; s=sen.Transpose(1,3); s=s.Transpose(2,4);}
0621   else if  ((i==o) && (j==p) && (k==n) && (l==m))
0622     {t=ten; s=sen.Transpose(1,4); s=s.Transpose(2,3); s=s.Transpose(3,4);}
0623   else if  ((i==p) && (j==m) && (k==n) && (l==o))
0624     {t=ten; s=sen.Transpose(1,2); s=s.Transpose(2,3); s=s.Transpose(3,4);}
0625   else if  ((i==p) && (j==m) && (k==o) && (l==n))
0626     {t=ten; s=sen.Transpose(1,2); s=s.Transpose(2,4);}
0627   else if  ((i==p) && (j==n) && (k==m) && (l==o))
0628     {t=ten; s=sen.Transpose(1,3); s=s.Transpose(3,4);}
0629   else if  ((i==p) && (j==n) && (k==o) && (l==m))
0630     {t=ten; s=sen.Transpose(1,4);}
0631   else if  ((i==p) && (j==o) && (k==m) && (l==n))
0632     {t=ten; s=sen.Transpose(1,3); s=s.Transpose(2,4); s=s.Transpose(3,4);}
0633   else if  ((i==p) && (j==o) && (k==n) && (l==m))
0634     {t=ten; s=sen.Transpose(1,4); s=s.Transpose(2,3);}
0635   else
0636     return PROMOTE(Scal1,Scal2)();
0637   PROMOTE(Scal1,Scal2) res(ten.at(0,0,0,0)-ten.at(0,0,0,0));
0638   for (unsigned short int a=0; a<4; ++a)
0639     for (unsigned short int b=0; b<4; ++b)
0640       for (unsigned short int c=0; c<4; ++c)
0641         for (unsigned short int d=0; d<4; ++d) {
0642           double ga(-1.),gb(-1.),gc(-1.),gd(-1.);
0643           if (a==0) ga=1.; if (b==0) gb=1.; if (c==0) gc=1.; if (d==0) gd=1.;
0644           res += ga*gb*gc*gd*t.at(a,b,c,d)*s.at(a,b,c,d);
0645         }
0646   return res;
0647 }
0648 
0649 template<typename Scal1, typename Scal2>
0650 Lorentz_Ten2<PROMOTE(Scal1,Scal2)>
0651 ATOOLS::Contraction(const Lorentz_Ten4<Scal1>& ten, int i, int j, int k,
0652                     const Lorentz_Ten4<Scal2>& sen, int l, int m, int n) {
0653   Lorentz_Ten4<Scal1> t; Lorentz_Ten4<Scal2> s;
0654   // T^{mu,nu} = t^{mu,alpha,beta,gamma} s^{nu,alpha',beta',gamma'}
0655   //                    g_{alpha,alpha'} g_{beta,beta'} g_{gamma,gamma'}
0656   if      (((i==l) && (j==m) && (k==n)) && ((i==2) && (j==3) && (k==4)))
0657     {t=ten; s=sen;}
0658   else if (((i==l) && (j==m) && (k==n)) && ((i==1) && (j==3) && (k==4)))
0659     {t=ten.Transpose(1,2); s=sen.Transpose(1,2);}
0660   else if (((i==l) && (j==m) && (k==n)) && ((i==1) && (j==2) && (k==4)))
0661     {t=ten.Transpose(1,3); s=sen.Transpose(1,3);}
0662   else if (((i==l) && (j==m) && (k==n)) && ((i==1) && (j==2) && (k==3)))
0663     {t=ten.Transpose(1,4); s=sen.Transpose(1,4);}
0664 
0665   else if (((i==l) && (j==n) && (k==m)) && ((i==2) && (j==3) && (k==4)))
0666     {t=ten; s=sen.Transpose(3,4);}
0667   else if (((i==l) && (j==n) && (k==m)) && ((i==1) && (j==3) && (k==4)))
0668     {t=ten.Transpose(1,2); s=sen.Transpose(1,2); s=s.Transpose(3,4);}
0669   else if (((i==l) && (j==n) && (k==m)) && ((i==1) && (j==2) && (k==4)))
0670     {t=ten.Transpose(1,3); s=sen.Transpose(1,3); s=s.Transpose(2,4);}
0671   else if (((i==l) && (j==n) && (k==m)) && ((i==1) && (j==2) && (k==3)))
0672     {t=ten.Transpose(1,4); s=sen.Transpose(1,4); s=s.Transpose(2,3);}
0673 
0674   else if (((i==m) && (j==l) && (k==n)) && ((i==2) && (j==3) && (k==4)))
0675     {t=ten; s=sen.Transpose(2,3);}
0676   else if (((i==m) && (j==l) && (k==n)) && ((i==1) && (j==3) && (k==4)))
0677     {t=ten.Transpose(1,2); s=sen.Transpose(1,2); s=s.Transpose(2,3);}
0678   else if (((i==m) && (j==l) && (k==n)) && ((i==1) && (j==2) && (k==4)))
0679     {t=ten.Transpose(1,3); s=sen.Transpose(1,3); s=s.Transpose(2,3);}
0680   else if (((i==m) && (j==l) && (k==n)) && ((i==1) && (j==2) && (k==3)))
0681     {t=ten.Transpose(1,4); s=sen.Transpose(1,4); s=s.Transpose(2,4);}
0682 
0683   else if (((i==m) && (j==n) && (k==l)) && ((i==2) && (j==3) && (k==4)))
0684     {t=ten; s=sen.Transpose(2,4); s=s.Transpose(3,4);}
0685   else if (((i==m) && (j==n) && (k==l)) && ((i==1) && (j==3) && (k==4)))
0686     {t=ten.Transpose(1,2); s=sen.Transpose(1,2); s=s.Transpose(2,4); s=s.Transpose(3,4);}
0687   else if (((i==m) && (j==n) && (k==l)) && ((i==1) && (j==2) && (k==4)))
0688     {t=ten.Transpose(1,3); s=sen.Transpose(1,3); s=s.Transpose(2,4); s=s.Transpose(2,3);}
0689   else if (((i==m) && (j==n) && (k==l)) && ((i==1) && (j==2) && (k==3)))
0690     {t=ten.Transpose(1,4); s=sen.Transpose(1,4); s=s.Transpose(2,4); s=s.Transpose(3,4);}
0691 
0692   else if (((i==n) && (j==l) && (k==m)) && ((i==2) && (j==3) && (k==4)))
0693     {t=ten; s=sen.Transpose(2,3); s=s.Transpose(3,4);}
0694   else if (((i==n) && (j==l) && (k==m)) && ((i==1) && (j==3) && (k==4)))
0695     {t=ten.Transpose(1,2); s=sen.Transpose(1,2); s=s.Transpose(2,3); s=s.Transpose(3,4);}
0696   else if (((i==n) && (j==l) && (k==m)) && ((i==1) && (j==2) && (k==4)))
0697     {t=ten.Transpose(1,3); s=sen.Transpose(1,3); s=s.Transpose(3,4); s=s.Transpose(2,3);}
0698   else if (((i==n) && (j==l) && (k==m)) && ((i==1) && (j==2) && (k==3)))
0699     {t=ten.Transpose(1,4); s=sen.Transpose(1,4); s=s.Transpose(2,3); s=s.Transpose(3,4);}
0700 
0701   else if (((i==n) && (j==m) && (k==l)) && ((i==2) && (j==3) && (k==4)))
0702     {t=ten; s=sen.Transpose(2,4);}
0703   else if (((i==n) && (j==m) && (k==l)) && ((i==1) && (j==3) && (k==4)))
0704     {t=ten.Transpose(1,2); s=sen.Transpose(1,2); s=s.Transpose(2,4);}
0705   else if (((i==n) && (j==m) && (k==l)) && ((i==1) && (j==2) && (k==4)))
0706     {t=ten.Transpose(1,3); s=sen.Transpose(1,3); s=s.Transpose(3,4);}
0707   else if (((i==n) && (j==m) && (k==l)) && ((i==1) && (j==2) && (k==3)))
0708     {t=ten.Transpose(1,4); s=sen.Transpose(1,4); s=s.Transpose(3,4);}
0709   else {
0710     THROW(fatal_error,"tensor contraction not specified for these indizes ...");
0711     return Lorentz_Ten2<PROMOTE(Scal1,Scal2)>();
0712   }
0713   Lorentz_Ten2<PROMOTE(Scal1,Scal2)> end;
0714   for (unsigned short int I=0; I<4; ++I)
0715     for (unsigned short int J=0; J<4; ++J) {
0716       PROMOTE(Scal1,Scal2) res(ten.at(0,0,0,0)-ten.at(0,0,0,0));
0717       for (unsigned short int a=0; a<4; ++a)
0718         for (unsigned short int b=0; b<4; ++b)
0719           for (unsigned short int c=0; c<4; ++c) {
0720             double ga(-1.),gb(-1.),gc(-1.);
0721             if (a==0) ga=1.; if (b==0) gb=1.; if (c==0) gc=1.;
0722             res += ga*gb*gc*t.at(I,a,b,c)*s.at(J,a,b,c);
0723           }
0724       end.at(I,J) = res;
0725     }
0726   return end;
0727 }
0728 
0729 template<typename Scal1, typename Scal2>
0730 Lorentz_Ten4<PROMOTE(Scal1,Scal2)>
0731 ATOOLS::Contraction(const Lorentz_Ten4<Scal1>& ten, int i, int j,
0732                     const Lorentz_Ten4<Scal2>& sen, int k, int l) {
0733   // T^{mu,nu,rho,sigma} = t^{mu,nu,alpha,beta} s^{rho,sigma,alpha',beta'}
0734   //                                        g_{alpha,alpha'} g_{beta,beta'}
0735   if      ((i!=3) && (j!=4) && (k!=3) && (l!=4)) {
0736     THROW(fatal_error,"tensor contraction not specified for these indizes ...");
0737     return Lorentz_Ten4<PROMOTE(Scal1,Scal2)>();
0738   }
0739   Lorentz_Ten4<PROMOTE(Scal1,Scal2)> end;
0740   for (unsigned short int I=0; I<4; ++I)
0741     for (unsigned short int J=0; J<4; ++J)
0742       for (unsigned short int K=0; K<4; ++K)
0743         for (unsigned short int L=0; L<4; ++L) {
0744           PROMOTE(Scal1,Scal2) res(ten.at(0,0,0,0)-ten.at(0,0,0,0));
0745           for (unsigned short int a=0; a<4; ++a)
0746             for (unsigned short int b=0; b<4; ++b) {
0747               double ga(-1.),gb(-1.);
0748               if (a==0) ga=1.; if (b==0) gb=1.;
0749               res += ga*gb*ten.at(I,J,a,b)*sen.at(K,L,a,b);
0750             }
0751           end.at(I,J) = res;
0752         }
0753   return end;
0754 }
0755 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!