Back to home page

EIC code displayed by LXR

 
 

    


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

0001 #ifndef ATOOLS_Math_Tensor_H
0002 #define ATOOLS_Math_Tensor_H
0003 
0004 #include <iostream>
0005 #include "ATOOLS/Math/Lorentz_Ten2.H"
0006 #include "ATOOLS/Math/Lorentz_Ten3.H"
0007 #include "ATOOLS/Math/Lorentz_Ten4.H"
0008 #include "ATOOLS/Math/Vector.H"
0009 #include "ATOOLS/Math/MyComplex.H"
0010 #include "ATOOLS/Math/MathTools.H"
0011 
0012 #include "ATOOLS/Org/Exception.H"
0013 #include "ATOOLS/Org/Message.H"
0014 
0015 
0016 namespace ATOOLS {
0017 
0018   /*
0019    * declarations and specialisations for 4 dimensional 2nd, 3rd and 4th rank tensors of doubles
0020    */
0021   typedef Lorentz_Ten2<double>  Lorentz_Ten2D;
0022   typedef Lorentz_Ten3<double>  Lorentz_Ten3D;
0023   typedef Lorentz_Ten4<double>  Lorentz_Ten4D;
0024 
0025   bool IsEqual(const Lorentz_Ten2D& t1, const Lorentz_Ten2D& t2, const double crit=1.0e-12);
0026   bool IsEqual(const Lorentz_Ten3D& t1, const Lorentz_Ten3D& t2, const double crit=1.0e-12);
0027   bool IsEqual(const Lorentz_Ten4D& t1, const Lorentz_Ten4D& t2, const double crit=1.0e-12);
0028 
0029   inline bool operator==(const Lorentz_Ten2D& t1, const Lorentz_Ten2D& t2) {
0030     return IsEqual(t1,t2);
0031   }
0032   inline bool operator==(const Lorentz_Ten3D& t1, const Lorentz_Ten3D& t2) {
0033     return IsEqual(t1,t2);
0034   }
0035   inline bool operator==(const Lorentz_Ten4D& t1, const Lorentz_Ten4D& t2) {
0036     return IsEqual(t1,t2);
0037   }
0038   inline bool operator!=(const Lorentz_Ten2D& t1, const Lorentz_Ten2D& t2) {
0039     return !IsEqual(t1,t2);
0040   }
0041   inline bool operator!=(const Lorentz_Ten3D& t1, const Lorentz_Ten3D& t2) {
0042     return !IsEqual(t1,t2);
0043   }
0044   inline bool operator!=(const Lorentz_Ten4D& t1, const Lorentz_Ten4D& t2) {
0045     return !IsEqual(t1,t2);
0046   }
0047 
0048   // **
0049   // declarations and specialisations for 4 dimensional 2nd, 3rd and 4th rank tensors of Complex
0050   // **
0051   typedef Lorentz_Ten2<Complex> Lorentz_Ten2C;
0052   typedef Lorentz_Ten3<Complex> Lorentz_Ten3C;
0053   typedef Lorentz_Ten4<Complex> Lorentz_Ten4C;
0054 
0055   bool IsEqual(const Lorentz_Ten2C& t1, const Lorentz_Ten2C& t2, const double crit=1.0e-12);
0056   bool IsEqual(const Lorentz_Ten3C& t1, const Lorentz_Ten3C& t2, const double crit=1.0e-12);
0057   bool IsEqual(const Lorentz_Ten4C& t1, const Lorentz_Ten4C& t2, const double crit=1.0e-12);
0058 
0059   inline bool operator==(const Lorentz_Ten2C& t1, const Lorentz_Ten2C& t2) {
0060     return IsEqual(t1,t2);
0061   }
0062   inline bool operator==(const Lorentz_Ten3C& t1, const Lorentz_Ten3C& t2) {
0063     return IsEqual(t1,t2);
0064   }
0065   inline bool operator==(const Lorentz_Ten4C& t1, const Lorentz_Ten4C& t2) {
0066     return IsEqual(t1,t2);
0067   }
0068   inline bool operator!=(const Lorentz_Ten2C& t1, const Lorentz_Ten2C& t2) {
0069     return !IsEqual(t1,t2);
0070   }
0071   inline bool operator!=(const Lorentz_Ten3C& t1, const Lorentz_Ten3C& t2) {
0072     return !IsEqual(t1,t2);
0073   }
0074   inline bool operator!=(const Lorentz_Ten4C& t1, const Lorentz_Ten4C& t2) {
0075     return !IsEqual(t1,t2);
0076   }
0077 
0078   inline Lorentz_Ten2C conj(const Lorentz_Ten2C& t) {
0079       Complex x[4][4];
0080       for (unsigned short int i=0; i<4; ++i)
0081         for (unsigned short int j=0; j<4; ++j)
0082           x[i][j] = std::conj(t.at(i,j));
0083       return Lorentz_Ten2C(x);
0084   }
0085 
0086   inline Lorentz_Ten3C conj(const Lorentz_Ten3C& t) {
0087       Complex x[4][4][4];
0088       Complex z = 0.;
0089       for (unsigned short int i=0; i<4; ++i)
0090         for (unsigned short int j=0; j<4; ++j)
0091           for (unsigned short int k=0; k<4; ++k)
0092             x[i][j][k] = std::conj(t.at(i,j,k));
0093       return Lorentz_Ten3C(x);
0094   }
0095 
0096   inline Lorentz_Ten4C conj(const Lorentz_Ten4C& t) {
0097       Complex x[4][4][4][4];
0098       for (unsigned short int i=0; i<4; ++i)
0099         for (unsigned short int j=0; j<4; ++j)
0100           for (unsigned short int k=0; k<4; ++k)
0101             for (unsigned short int l=0; l<4; ++l)
0102             x[i][j][k][l] = std::conj(t.at(i,j,k,l));
0103       return Lorentz_Ten4C(x);
0104   }
0105 
0106 
0107   // meaningful only for the transpostion of a pair of indizes (i<->j)
0108   inline Lorentz_Ten2C hermconj(const Lorentz_Ten2C& t) {
0109     return conj(t).Transpose();
0110   }
0111 
0112   inline Lorentz_Ten3C hermconj(const Lorentz_Ten3C& t, unsigned short int i, unsigned short int j) {
0113     return conj(t).Transpose(i,j);
0114   }
0115 
0116   inline Lorentz_Ten4C hermconj(const Lorentz_Ten4C& t, unsigned short int i, unsigned short int j) {
0117     return conj(t).Transpose(i,j);
0118   }
0119 
0120 
0121 
0122 
0123 
0124   /************************************************************************************/
0125   // define special tensors
0126   // all tensor defined with upper indices
0127 
0128   // metric tensor g^{mu,nu} = g_{mu,nu}, componentwise
0129 
0130   inline Lorentz_Ten2D MetricTensor() {
0131     return Lorentz_Ten2D( 1., 0., 0., 0.,
0132                           0.,-1., 0., 0.,
0133                           0., 0.,-1., 0.,
0134                           0., 0., 0.,-1.);
0135   }
0136 
0137   inline Lorentz_Ten4D EpsilonTensor() {
0138     return Lorentz_Ten4D(
0139         0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0140         0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,-1., 0.,
0141         0., 0., 0., 0., 0., 0., 0.,-1., 0., 0., 0., 0., 0., 1., 0., 0.,
0142         0., 0., 0., 0., 0., 0., 1., 0., 0.,-1., 0., 0., 0., 0., 0., 0.,
0143 
0144         0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,-1., 0., 0., 1., 0.,
0145         0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0146         0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.,-1., 0., 0., 0.,
0147         0., 0.,-1., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,
0148 
0149         0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.,-1., 0., 0.,
0150         0., 0., 0.,-1., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.,
0151         0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0152         0., 1., 0., 0.,-1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0153 
0154         0., 0., 0., 0., 0., 0.,-1., 0., 0., 1., 0., 0., 0., 0., 0., 0.,
0155         0., 0., 1., 0., 0., 0., 0., 0.,-1., 0., 0., 0., 0., 0., 0., 0.,
0156         0.,-1., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0157         0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.);
0158   }
0159 
0160   /***************************************************************************/
0161   // build tensor from vectors and other lower rank tensors
0162 
0163   template<typename Scal1, typename Scal2>                                // ok
0164   Lorentz_Ten2<PROMOTE(Scal1,Scal2)> 
0165   BuildTensor(const Vec4<Scal1>&, const Vec4<Scal2>&);
0166 
0167   template<typename Scal1, typename Scal2, typename Scal3>                // ok
0168   Lorentz_Ten3<PROMOTE(Scal1,PROMOTE(Scal2,Scal3))> 
0169   BuildTensor(const Vec4<Scal1>&, const Vec4<Scal2>&, const Vec4<Scal3>&);
0170 
0171   template<typename Scal1, typename Scal2>                                // ok
0172   Lorentz_Ten3<PROMOTE(Scal1,Scal2)> 
0173   BuildTensor(const Lorentz_Ten2<Scal1>&, const Vec4<Scal2>&);
0174 
0175   template<typename Scal1, typename Scal2, typename Scal3, typename Scal4>// ok
0176   Lorentz_Ten4<PROMOTE(Scal1,PROMOTE(Scal2,PROMOTE(Scal3,Scal4)))> 
0177   BuildTensor(const Vec4<Scal1>&, const Vec4<Scal2>&,
0178               const Vec4<Scal3>&, const Vec4<Scal4>&);
0179 
0180   template<typename Scal1, typename Scal2>                                // ok
0181   Lorentz_Ten4<PROMOTE(Scal1,Scal2)> 
0182   BuildTensor(const Lorentz_Ten2<Scal1>&, const Lorentz_Ten2<Scal2>&);
0183 
0184   template<typename Scal1, typename Scal2>                                // ok
0185   Lorentz_Ten4<PROMOTE(Scal1,Scal2)> 
0186   BuildTensor(const Lorentz_Ten3<Scal1>&, const Vec4<Scal2>&);
0187 
0188 
0189   /***************************************************************************/
0190   // tensor contractions
0191 
0192   // epsilon tensor contractions  
0193 
0194   template<typename Scalar>                                               // ok
0195   Lorentz_Ten3<Scalar>
0196   EpsilonTensorContraction(const Vec4<Scalar>&);
0197 
0198   template<typename Scal1, typename Scal2>                                // ok
0199   Lorentz_Ten2<PROMOTE(Scal1,Scal2)>
0200   EpsilonTensorContraction(const Vec4<Scal1>&, const Vec4<Scal2>&);
0201 
0202   template<typename Scal1, typename Scal2, typename Scal3>                // ok
0203   Vec4<PROMOTE(Scal1,PROMOTE(Scal2,Scal3))>
0204   EpsilonTensorContraction(const Vec4<Scal1>&, const Vec4<Scal2>&,
0205                            const Vec4<Scal3>&);
0206 
0207   template<typename Scal1, typename Scal2, typename Scal3, typename Scal4>// ok
0208   PROMOTE(Scal1,PROMOTE(Scal2,PROMOTE(Scal3,Scal4)))
0209   EpsilonTensorContraction(const Vec4<Scal1>&, const Vec4<Scal2>&,
0210                            const Vec4<Scal3>&, const Vec4<Scal4>&);
0211 
0212   template<typename Scalar>                                               // ok
0213   Lorentz_Ten2<Scalar>
0214   EpsilonTensorContraction(const Lorentz_Ten2<Scalar>&, int, int);
0215 
0216   template<typename Scal1, typename Scal2>                                // ok
0217   Vec4<PROMOTE(Scal1,Scal2)>
0218   EpsilonTensorContraction(const Lorentz_Ten2<Scal1>&, int, int,
0219                            const Vec4<Scal2>&, int);
0220 
0221   template<typename Scal1, typename Scal2>                                // ok
0222   PROMOTE(Scal1,Scal2)
0223   EpsilonTensorContraction(const Lorentz_Ten2<Scal1>&, int, int,
0224                            const Lorentz_Ten2<Scal2>&, int, int);
0225 
0226   template<typename Scalar>                                               // ok
0227   Vec4<Scalar>
0228   EpsilonTensorContraction(const Lorentz_Ten3<Scalar>&, int, int, int);
0229 
0230   template<typename Scal1, typename Scal2>                                // ok
0231   PROMOTE(Scal1,Scal2)
0232   EpsilonTensorContraction(const Lorentz_Ten3<Scal1>&, int, int, int,
0233                            const Vec4<Scal2>&, int);
0234 
0235   template<typename Scalar>                                               // ok
0236   Scalar
0237   EpsilonTensorContraction(const Lorentz_Ten4<Scalar>&, int, int, int, int);
0238 
0239 
0240   // others
0241 
0242   // tensor with itself
0243   template<typename Scalar>                                               // ok
0244   Scalar
0245   Contraction(const Lorentz_Ten2<Scalar>&);
0246 
0247   template<typename Scalar>                                               // ok
0248   Vec4<Scalar>
0249   Contraction(const Lorentz_Ten3<Scalar>&, int, int);
0250 
0251   template<typename Scalar>                                               // ok
0252   Lorentz_Ten2<Scalar>
0253   Contraction(const Lorentz_Ten4<Scalar>&, int, int);
0254 
0255   template<typename Scalar>                                               // ok
0256   Scalar
0257   Contraction(const Lorentz_Ten4<Scalar>&, int, int, int, int);
0258 
0259 
0260   // tensor with vector
0261   template<typename Scal1, typename Scal2>                                // ok
0262   Vec4<PROMOTE(Scal1,Scal2)>
0263   Contraction(const Lorentz_Ten2<Scal1>&, int, const Vec4<Scal2>&);
0264 
0265   template<typename Scal1, typename Scal2, typename Scal3>                // ok
0266   PROMOTE(PROMOTE(Scal1,Scal2),Scal3)
0267   Contraction(const Lorentz_Ten2<Scal1>&, const Vec4<Scal2>&,
0268               const Vec4<Scal3>&);
0269 
0270   template<typename Scal1, typename Scal2>                                // ok
0271   Lorentz_Ten2<PROMOTE(Scal1,Scal2)>
0272   Contraction(const Lorentz_Ten3<Scal1>&, int, const Vec4<Scal2>&);
0273 
0274   template<typename Scal1, typename Scal2>                                // ok
0275   Lorentz_Ten3<PROMOTE(Scal1,Scal2)>
0276   Contraction(const Lorentz_Ten4<Scal1>&, int, const Vec4<Scal2>&);
0277 
0278 
0279   // tensor with other tensor
0280   template<typename Scal1, typename Scal2>                                // ok
0281   PROMOTE(Scal1,Scal2)
0282   Contraction(const Lorentz_Ten2<Scal1>&, int, int,
0283               const Lorentz_Ten2<Scal2>&, int, int);
0284 
0285   template<typename Scal1, typename Scal2>                                // ok
0286   Lorentz_Ten2<PROMOTE(Scal1,Scal2)>
0287   Contraction(const Lorentz_Ten2<Scal1>&, int, const Lorentz_Ten2<Scal2>&, int);
0288 
0289   template<typename Scal1, typename Scal2>                                // ok
0290   Vec4<PROMOTE(Scal1,Scal2)>
0291   Contraction(const Lorentz_Ten3<Scal1>&, int, int,
0292               const Lorentz_Ten2<Scal2>&, int, int);
0293 
0294   template<typename Scal1, typename Scal2>                                // ok
0295   Lorentz_Ten3<PROMOTE(Scal1,Scal2)>
0296   Contraction(const Lorentz_Ten3<Scal1>&, int,
0297               const Lorentz_Ten2<Scal2>&, int);
0298 
0299   template<typename Scal1, typename Scal2>
0300   PROMOTE(Scal1,Scal2)
0301   Contraction(const Lorentz_Ten3<Scal1>&, int, int, int,
0302               const Lorentz_Ten3<Scal2>&, int, int, int);
0303 
0304   template<typename Scal1, typename Scal2>
0305   Lorentz_Ten2<PROMOTE(Scal1,Scal2)>
0306   Contraction(const Lorentz_Ten3<Scal1>&, int, int,
0307               const Lorentz_Ten3<Scal2>&, int, int);
0308 
0309   template<typename Scal1, typename Scal2>
0310   Lorentz_Ten4<PROMOTE(Scal1,Scal2)>
0311   Contraction(const Lorentz_Ten3<Scal1>&, int,
0312               const Lorentz_Ten3<Scal2>&, int);
0313 
0314   template<typename Scal1, typename Scal2>
0315   Lorentz_Ten2<PROMOTE(Scal1,Scal2)>
0316   Contraction(const Lorentz_Ten4<Scal1>&, int, int,
0317               const Lorentz_Ten2<Scal2>&, int, int);
0318 
0319   template<typename Scal1, typename Scal2>
0320   Lorentz_Ten4<PROMOTE(Scal1,Scal2)>
0321   Contraction(const Lorentz_Ten4<Scal1>&, int,
0322               const Lorentz_Ten2<Scal2>&, int);
0323 
0324   template<typename Scal1, typename Scal2>                                //(tbc)
0325   Vec4<PROMOTE(Scal1,Scal2)>
0326   Contraction(const Lorentz_Ten4<Scal1>&, int, int, int,
0327               const Lorentz_Ten3<Scal2>&, int, int, int);
0328 
0329   template<typename Scal1, typename Scal2>                                //(tbc)
0330   Lorentz_Ten3<PROMOTE(Scal1,Scal2)>
0331   Contraction(const Lorentz_Ten4<Scal1>&, int, int,
0332               const Lorentz_Ten3<Scal2>&, int, int);
0333 
0334   template<typename Scal1, typename Scal2>                                // tbc
0335   PROMOTE(Scal1,Scal2)
0336   Contraction(const Lorentz_Ten4<Scal1>&, int, int, int, int,
0337               const Lorentz_Ten4<Scal2>&, int, int, int, int);
0338 
0339   template<typename Scal1, typename Scal2>                                // tbc
0340   Lorentz_Ten2<PROMOTE(Scal1,Scal2)>
0341   Contraction(const Lorentz_Ten4<Scal1>&, int, int, int,
0342               const Lorentz_Ten4<Scal2>&, int, int, int);
0343 
0344   template<typename Scal1, typename Scal2>                                //(tbc)
0345   Lorentz_Ten4<PROMOTE(Scal1,Scal2)>
0346   Contraction(const Lorentz_Ten4<Scal1>&, int, int,
0347               const Lorentz_Ten4<Scal2>&, int, int);
0348 
0349 }
0350 
0351 //*********************************************************************************
0352 //*********************************************************************************
0353 // tensor algebra implementation
0354 //*********************************************************************************
0355 
0356 using namespace std;
0357 using namespace ATOOLS;
0358 
0359 // build tensors
0360 
0361 #include "ATOOLS/Math/Tensor_Build.H"
0362 
0363 // tensor contractions
0364 
0365 #include "ATOOLS/Math/Tensor_Contractions_Epsilon.H"
0366 
0367 // others
0368 
0369 #include "ATOOLS/Math/Tensor_Contractions.H"
0370 
0371 #endif