(***********************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which typically
will have the same name as this file except ending in ".nb" instead of
".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent to
using the Evaluate Initialization Cells menu command in the front end.

DO NOT EDIT THIS FILE.  This entire file is regenerated automatically 
each time the parent Notebook file is saved in the Mathematica front end.
Any changes you make to this file will be overwritten.
***********************************************************************)

(*:Mathematica Version: 4.0  *)

(* :Copyright: Copyright 1989-1999, Wolfram Research, Inc.*)

(*:Context: LinearAlgebra`TensorManipulation`  *)

(*:Title: Basic Tensor Manipulation Functions  *)

(*:Authors:   Juan Ramon Ruiz Tolosa and Enrique Castillo *)

(*:Summary:
This package provides basic tensor manipulation routines
*)


(*:Limitations:

*)

(*:Sources:   
J. R. Ruiz-Tolosa and E. Castillo.
From vectors to Tensors
*)

(*:Keywords: Tensors *)

(*:Requirements: None.  *)

(*:Package Version: 1.0 *)

(*:History:   Created June 7, 2003 by J. R. Ruiz-Tolosa and E. Castillo
*)

BeginPackage["LinearAlgebra`TensorManipulation`"];

MatrixTensorProduct::usage=
"MatrixTensorProduct[A,B] gives a new matrix that is the tensor product of 
matrices A and B. The matrices A and B can be of any dimension."

ListMatrixTensorProduct::usage=
"ListMatrixTensorProduct[{A, B,...}] gives a new matrix that is the tensor product of 
matrices A, B, etc., i.e., of all matrices in the list in the given order. The matrices in the 
list can be of any dimension."

SchurProd::usage=
"SchurProd[A, B] gives the Schur product of two matrices. A and B must have the same dimensions."

P1CTP::usage=
"P1CTP[m,p] returns a matrix P1 such that the tensor product of any two matrices B and A of 
dimensions {p,q} and {m,n}, respectively is identical to P1 M P2, where M is the tensor product of 
A and B."

P2CTP::usage=
"P2CTP[m,p] returns a matrix P2 such that the tensor product of any two matrices B and A of 
dimensions {p,q} and {m,n}, respectively is identical to P1 M P2, where M is the tensor product of 
A and B."

CommuteTensorProductMatrices::usage=
"CommuteTensorProductMatrices[m,n,p,q] gives two matrices P1 and P2 such that the tensor product
 of any two matrices B and A of dimensions {p,q} and {m,n}, respectively is identical to P1 M P2, where
 M is the tensor product of A and B."

QTPSchurMatrix::usage=
"QTPSchurMatrix[n] returns a matrix Q of dimensions {m,m*m} such that the Schur produc of any pair
of matrices A and B of dimension {m,n} is the product of Q M P, where M is the tensor product of 
A and B, and P is a matrix of dimension {n*n,n}."

PTPSchurMatrix::usage=
"PTPSchurMatrix[n] returns a matrix P of dimensions {n*n,n} such that the Schur produc of any pair
of matrices A and B of dimension {m,n} is the product of Q M P, where M is the tensor product of
A and B, and Q is a matrix of dimension {m,m*m}."

QPTPSchurMatrices::usage=
"QPTPSchurMatrices[n,m] returns a list {Q,P} of two matrices of dimensions {m,m*m} and {n*n,n}, 
respectively, such that the Schur produc of any matrices A and B of dimension {m,n} is  the product 
of Q M P, where M is the tensor product of A and B."

TensorMatrixForm::usage=
"TensorMatrixForm[tensor,nrowindices] returns the tensor components in matrix form grouping the first 
nrowindices indices as rows."

ExtensionMatrix::usage=
"ExtensionMatrix[G] returns the extended (column) form of matrix G."

CondensationMatrix::usage=
"CondensationMatrix[A,m,n] returns the original form of a matrix with m rows and n columns from 
its extended form. The product mn must be identical to the column length."

VectorTensorProduct::usage=
"VectorTensorProduct[VList] returns a list with two lists: the extended list of the components of the tensor
 product of the vectors in VList, and the list of its ranges."

TensorProductComponents::usage=
"TensorProductComponents[VList] returns the components of the tensor product of the vectors in 
VList together with their scores and component indices."

HomogeneousTensorComponents::usage=
"HomogeneousTensorComponents[homogeneoustensorc] returns the components of the homogeneous tensor
together with their scores and component indices. An homogeneous tensor is given by a list of two
 lists: the list of all components and the list of its linear space dimension and its order."

RowTensorComponentIndex::usage=
"RowTensorComponentIndex[n,r] gives the tensor component indices corresponding to the n-th component
of a tensor which associated linear spaces dimensions are in the list r."

TensorComponentRowIndex::usage=
"TensorComponentRowIndex[r,in] gives the score of the component
of a tensor which associated linear spaces dimensions are in the list r."

TensorContractionMatrix::usage=
"TensorContractionMatrix[n,s,p,q] returns the matrix C that allows the contraction of indices 
p and q of an homogeneous tensor of order s defined on a linear space of dimension n, i.e., R=C T, where R and T are the extended forms of the contracted and initial 
tensors, respectively."

TensorContractionHomogeneous::usage=
"TensorContractionHomogeneous[H,p,q] returns the homogeneous tensor resulting from contracting
indices p and q of the homogeneous tensor H."

TensorContraction::usage=
"TensorContraction[H,p,q] returns the homogeneous tensor resulting from contracting
indices p and q of the tensor H."

TensorContractionListHomogeneous::usage=
"TensorContractionListHomogeneous[H,ContList] returns the homogeneous tensor resulting after
contracting the pairs of indices in ContList of the homogeneous tensor H."

TensorContractionList::usage=
"TensorContractionList[H,ContList] returns the homogeneous tensor resulting after
contracting the pairs of indices in ContList of the homogeneous tensor H."

CovaContra::usage=
"CovaContra[G,x] Returns the covariant components of a vector with contravariant components x. G is the Gram matrix."

ContraCova::usage=
"ContraCova[G,x] Returns the contravariant components of a vector with covariant components x. G is the Gram matrix."

CovaContraBaseChange::usage=
"CovaContraBaseChange[G,C,x] Returns the covariant components of a vector with contravariant components x after a 
basis change with associated matrix C. G is the Gram matrix."

CovaCovaBaseChange::usage=
"CovaCovaBaseChange[G,C,x] Returns the covariant components of a vector with covariant components x after a 
basis change with associated matrix C. G is the Gram matrix."

ContraCovaBaseChange::usage=
"ContraCovaBaseChange[G,C,x] Returns the contravariant components of a vector with covariant components x after a 
basis change with associated matrix C. G is the Gram matrix."

ContraContraBaseChange::usage=
"ContraContraBaseChange[C,x] Returns the contravariant components of a vector with contravariant components x after a 
basis change with associated matrix C."

GrammMatrixBaseChange::usage=
"GrammMatrixBaseChange[G,C] returns the new Gram matrix after a basis change of associated matrix C."

BaseTensorChange::usage=
"BaseTensorChange[CM,CovariancesList] returns the matrix that allows calculate the extended component matrix of a tensor
of covariaces given in CmatricesList from its extended component matrix after changing the basis in the initial linear 
spaces using the matrices in the list CM."

InverseBaseTensorChange::usage=
"InverseBaseTensorChange[CM,CovariancesList] returns the matrix associated with the inverse of the transformation
given by BaseTensorChange[CM,CovariancesList]."

ListEigenValuesVectors::usage=
"EigenValuesVectors[MList] returns a list with the lists of eigenvalues and eigenvectors of the matrices in
MList."

TensorProductEigenValuesVectors::usage=
"TensorProductEigenValuesVectors[MList] returns a list with the lists of eigenvalues and eigenvectors of the
tensor product of the matrices in MList."

InitCounter::usage=
"InitCounter[r] initializes a counter with r ones."

AddCounter::usage=
"AddCounter[in,r] Adds one unit to the counter in a multibase format."

TensorComponents::usage=
"TensorComponents[tensor] Prints the components of a tensor with their corresponding indices."

ListTensorProduct::usage=
"ListTensorProduct[TList] performs the tensor product of the tensors in the list."
    
TensorProduct::usage=
"TensorProduct[tensor1,tensor2]:=performs the tensor product of the tensors tensor1 and tensor2."

TensorListProductComponents::usage=
"TensorListProductComponents[VList] returns the tensor product of the tensors in the given list."

VectorListTensorProduct::usage=
"VectorListTensorProduct[VList] returns the tensor product of the tensors in the given list."

TensorProductInversion::usage=
"TensorProductInversion[A,B] obtains the tensor product of the inverse matrices A and B."

IsTensor::usage=
"IsTensor[tensor] tests if tensor is a tensor."

HomogeneousTensorComponents::usage=
"HomogeneousTensorComponents[homogeneoustensorc] returns the components of a homogeneous tensor."

ValidContraction::usage=
"ValidContraction[tensor,p,q] Tests if the contraction of indices p and q is valid for tensor."

NewDualBasis::usage=
"NewDualBasis[OldBasis,CMatrix] returns the new dual basis after a change of matrix CMatrix."

NewContComponents::usage=
"NewContComponents[OldContComponents,CMatrix] returns the new contravariant components after 
a change of matrix CMatrix."

NewCovComponents::usage=
"NewCovComponents[OldCovComponents,CMatrix] returns the new covariant components after 
a change of matrix CMatrix."

NewGramMatrix::usage=
"NewGramMatrix[CMatrix] returns the new gram matrix after a change of matrix CMatrix."

TemplateSymmetricTensor::usage=
"TemplateSymmetricTensor[order,n] returns a template of the indicated symmetric tensor."

AddSymmetricCounter::usage=
"AddSymmetricCounter[insim,r] adds one unit to the insim symmetric counter."

SymmetrizeTensorCompact::usage=
"SymmetrizeTensorCompact[tensor] returns the symmetrized compacted tensor of tensor."
 
SymmetrizeTensor::usage=
"SymmetrizeTensor[tensor] returns the symmetrized tensor of tensor."
  
SymmetricTensorToTensor::usage=
"SymmetricTensorToTensor[symmetrictensor] converts the symmetric tensor to a general tensor format."
  
IsSymmetricTensor::usage=
"IsSymmetricTensor[tensor] test if the tensor is symmetric."

SymmetricTensorStrictComponents::usage=
"SymmetricTensorStrictComponents[tensor] returns the strict components of the given symmetric tensor."

AddAntiSymmetricCounter::usage=
"AddAntiSymmetricCounter[insim,r] adds one unit to the insim antisymmetric counter."

InitAntiSymmetricCounter::usage=
"InitAntiSymmetricCounter[order,n] initializes an antisymmetric counter."

TemplateAntiSymmetricTensor::usage=
"TemplateAntiSymmetricTensor[order,n] returns a template of the indicated antisymmetric tensor."

AntiSymmetrizeTensorCompact::usage=
"AntiSymmetrizeTensorCompact[tensor] returns the antisymmetrized compacted tensor of tensor."

AntiSymmetrizeTensor::usage=
"AntiSymmetrizeTensor[tensor] returns the antisymmetrized tensor of tensor."

AntySymmetricTensorToTensor::usage=
"AntySymmetricTensorToTensor[antisymmetrictensor] converts the antisymmetric tensor to a general tensor format."

IsAntiSymmetricTensor::usage=
"IsAntiSymmetricTensor[tensor] test if the tensor is antisymmetric."

AntiSymmetricTensorStrictComponents::usage=
"AntiSymmetricTensorStrictComponents[tensor] returns the strict components of the given antisymmetric tensor."
 
DeltaK::usage=
"DeltaK[i,j] retuns the Kronecker delta value."

GeneralizedKroneckerDelta::usage=
"GeneralizedKroneckerDelta[upperindices,lowerindices] retuns the generalized Kronecker delta value."

LeviCivita::usage=
"LeviCivita[indices] retuns the Levi Civita value."

CreateLeviCivita::usage=
"CreateLeviCivita[order,n,var] generates a Levi Civita tensor."
 
AdjTensor::usage=
"AdjTensor[tensor] returns the adjoint tensor of the given tensor."
 
IndexExchange::usage=
"IndexExchange[tensor,changes,Gram] exchanges indices in a tensor."
 
InitStrictCounter::usage=
"InitStrictCounter[order] initializes an strict components counter."

AddStrictCounter::usage=
"AddStrictCounter[insim,r] add one unit to an strict components counter."

TemplateStrictTensor::usage=
"TemplateStrictTensor[order,n] returns the template of the given tensor strict components."
 
ExteriorProduct::usage=
"ExteriorProduct[vectors] performs the exterior product of the given vectors."
 
EuclideanTensorScalarProduct::usage=
"EuclideanTensorScalarProduct[tensor1,tensor2,Gram] performs the Euclidean tensor scalar product of tensor1 and tensor2 using the given Gram matrix."
 
EuclideanTensorScalarProductMatrix::usage=
"EuclideanTensorScalarProductMatrix[tensor1,tensor2,Gram]  performs the Euclidean tensor scalar product of tensor1 and tensor2 using the given Gram matrix."
 
TensorModulus::usage=
"TensorModulus[tensor1,Gram] returns the modulus of the given tensor."
 
TensorAngle::usage=
"TensorAngle[tensor1,tensor2,Gram] returns the angle between the two given tensors."
 
EuclideanAdjTensor::usage=
"EuclideanAdjTensor[tensor,Gram] returns the Euclidean adjoint of the given tensor."
 
ProjectionTensor::usage=
"ProjectionTensor[Direction] returns the proyection tensor associated with the given direction."

ProjectionComplementTensor::usage=
"ProjectionComplementTensor[Direction] returns the complement of the proyection tensor associated with the given direction."

MomentTensor::usage=
"MomentTensor[r]  returns the moment tensor associated with the given vector."

AxisMomentTensor::usage=
"AxisMomentTensor[r,axis] returns the moment tensor associated with the given vector and axis."

RotationTensor::usage=
"RotationTensor[angle,axis] returns the rotation tensor associated with the given angle and axis."

ReflexionTensor::usage=
"ReflexionTensor[axis] returns the reflection tensor associated with the given axis."

AffinityTensor::usage=
"AffinityTensor[A,L] returns the affinity tensor associated with the given A and L."

HomoteticTensor::usage=
"HomoteticTensor[center,k] returns the homotetic tensor associated with the given center and k."

TranslationTensor::usage=
"TranslationTensor[t] returns the translation tensor associated with the given translation t."

RotationTensor2D::usage=
"RotationTensor2D[angle] returns the rotation 2D tensor associated with the given angle."

CentralSymmetryTensor::usage=
"CentralSymmetryTensor[center] returns the central symmetry tensor associated with the given center."

AxialSymmetryTensor2D::usage=
"AxialSymmetryTensor2D[A,B,C1] returns the axial symmetry tensor associated with the given A, B and C1."

Permute::usage=
"Permute[in, perm] permutes the list according to permutation perm."

Permutation::usage=
"Permutation[tensor,perm] Permutes the tensor according to the permutation perm."

(**********************************************************************************************)

Begin["`Private`"]

Procedures
(* CHAPTER 1 *)
CovaContra[GrammMatrix_,x_]:=Inverse[GrammMatrix].x
  
ContraCova[GrammMatrix_,x_]:=GrammMatrix.x
  
CovaContraBaseChange[GrammMatrix_,C1_,x_]:=Inverse[GrammMatrix.C1].x
  
CovaCovaBaseChange[GrammMatrix_,C1_,x_]:=Transpose[C1].x
  
ContraCovaBaseChange[GrammMatrix_,C1_,x_]:=
  Inverse[C1.Transpose[C1].GrammMatrix.C1].x
  
ContraContraBaseChange[C1_,x_]:=Inverse[C1].x
  
GrammMatrixBaseChange[G_,C1_]:=Transpose[C1].G.C1
  
MatrixTensorProduct[A_,B_]:=
  Module[{H,i,j,k,l,m=Length[A],n=Length[A[[1]]],p=Length[B],
      q=Length[B[[1]]]},H=Table[0,{i,1,m*p},{j,1,n*q}];
    Do[ii=(i-1)*p;
      Do[jj=(j-1)*q;
        Do[H[[ii+k,jj+l]]=A[[i,j]]*B[[k,l]];,{k,1,p},{l,1,q}],{j,1,n}],{i,1,
        m}];
    Return[H]]

ListMatrixTensorProduct[MList_]:=
  Module[{i,n=Length[MList],aux},aux=IdentityMatrix[1];
    Do[aux=MatrixTensorProduct[aux,MList[[i]]],{i,1,n}];
    Return[aux]]

SchurProd::invalid="`1` has not the same dimensions than `2`.";

SchurProd[A_,B_]:=Module[{C1,i,j},C1=A;
      Do[C1[[i,j]]=A[[i,j]]*B[[i,j]],{i,1,Length[A]},{j,1,Length[A[[1]]]}];
      Return[C1]]/;SameDimensions[A,B]
    
SameDimensions[A_,B_]:=
  Module[{aux},aux=Dimensions[A]\[Equal]Dimensions[B];
    If[!aux,\
          Message[SchurProd::invalid,A,B]];Return[aux]]

TensorProductInversion[A_,B_]:=MatrixTensorProduct[Inverse[A],Inverse[B]]

P1CTP[m_,p_]:=Module[{P,i,k},P=Table[0,{i,1,m*p},{k,1,m*p}];
    Do[P[[(k-1)*m+i,(i-1)*p+k]]=1,{i,1,m},{k,1,p}];
    Return[P]]

P2CTP[n_,q_]:=Module[{P,j,l},P=Table[0,{j,1,n*q},{l,1,n*q}];
    Do[P[[(j-1)*q+l,(l-1)*n+j]]=1,{j,1,n},{l,1,q}];
    Return[P]]

(*P2CTP[n_,q_]:=Transpose[P1CTP[n,q]]*)

CommuteTensorProductMatrices[m_,n_,p_,q_]:={P1CTP[m,p],P2CTP[n,q]}

PTPSchurMatrix[n_]:=Module[{P1,i,j},P1=Table[0,{j,1,n^2},{i,1,n}];
    Do[P1[[n*(j-1)+j,j]]=1,{j,1,n}];
    Return[P1]]

QTPSchurMatrix[m_]:=Module[{Q1,i,j},Q1=Table[0,{i,1,m},{j,1,m^2}];
    Do[Q1[[i,m*(i-1)+i]]=1,{i,1,m}];
    Return[Q1]]

QPTPSchurMatrices[n_,m_]:={QTPSchurMatrix[m],PTPSchurMatrix[n]}

ExtensionMatrix[G_?MatrixQ]:=Module[{aux,i},aux=Flatten[G];
    Return[Table[{aux[[i]]},{i,1,Length[aux]}]]]

CondensationMatrix[A_,m_,n_]:=Module[{B,i,j,k=1},B=Table[0,{i,1,m},{j,1,n}];
      Do[B[[i,j]]=A[[k,1]];k+=1,{i,1,m},{j,1,n}];Return[B];]/;
    CondensationMatrixCheck[A,m,n]
    
CondensationMatrixCheck[A_,m_,n_]:=m*n\[Equal]Dimensions[A][[1]]
  
  
(*CHAPTER  2*)
IsTensor::invalidtensor="`1` is not a tensor.";

IsTensor[tensor_]:=Module[{aux=False,aux1,m,i},
    If[ListQ[tensor],
      If[Length[tensor]\[Equal]2,
        If[ListQ[tensor[[1]]],
          If[ListQ[tensor[[2]]],
            
            If[Length[tensor[[1]]]==
                Product[Abs[tensor[[2,i]]],{i,1,Length[tensor[[2]]]}],
              aux=True;
              If[aux,
                i=1;aux1=True;
                m=Length[tensor[[2]]];
                While[i<=m && aux1,
                  
                  If[Length[tensor[[2,i]]]>0 && IntegerQ[tensor[[2,i]]],
                    aux1=False];
                  i++;
                  ];
                aux=aux&&aux1;
                ]
              ]
            ]
          ]
        ]
      ];
    If[!aux,Message[IsTensor::invalidtensor,tensor]];
    Return[aux]
    ]

(*Initializes a counter to a list of r ones*)

InitCounter[r_]:=Table[1,{i,1,r}]

(*Adds one unit to the counter in a multibase format*)

AddCounter[in_List,r_List]:=Module[{Done=False,rr=Length[r],in1=in,i,j},i=rr;
    While[!Done&&i\[GreaterEqual]1,If[in1[[i]]<r[[i]],Done=True;in1[[i]]+=1;
          Do[in1[[j]]=1,{j,i+1,rr}],i-=1];];
    Return[in1]]

(*Calculates the dimension of the tensor product of the vectors in VList*)

CalculateSizeTensProd[VList_]:=
  Product[Length[VList[[i]]],{i,1,Length[VList]}]

(*Returns the list of dimensions of the linear spaces of the vectors in\
     
       VList*)

CalculateRangesTensProd[VList_]:=
  Table[Length[VList[[i]]],{i,1,Length[VList]}]

TensorComponents[tensor_]:=
  Module[{in,i,range,ranges,ranges1,Size,comp={},count=0},
      If[tensor[[2]]\[Equal]{},Return[{1,{0},tensor[[1]]}],
        range=Length[tensor[[2]]];
        in=InitCounter[range];
        ranges=Abs[tensor[[2]]];
        ranges1=tensor[[2]];
        Size=Product[ranges[[i]],{i,1,Length[ranges]}];
        Do[count=count+1;
          aa={in,in};
          Do[If[ranges1[[j]]>0,aa[[2,j]]="o",aa[[1,j]]="o"],{j,1,range}];
          AppendTo[comp,{count,MatrixForm[aa],tensor[[1,i]]}];
          in=AddCounter[in,ranges],{i,1,Size}];
        Return[MatrixForm[comp]]]]/;IsTensor[tensor]
    
VectorListTensorProduct[VList_]:=
  Module[{in,i,ranges,Size,comp={},count=0,aux},
    in=InitCounter[Length[VList]];
    ranges=CalculateRangesTensProd[VList];
    Size=CalculateSizeTensProd[VList];
    Do[count+=1;
      aux=Product[VList[[j,in[[j]]]],{j,1,Length[VList]}];
      AppendTo[comp,aux];
      in=AddCounter[in,ranges];,{i,1,Size}];
    Return[{comp,ranges}]]

TensorListProductComponents[VList_]:=
  Module[{in,i,ranges,Size,comp={},count=0,aux},
    in=InitCounter[Length[VList]];
    ranges=CalculateRangesTensProd[VList];
    Size=CalculateSizeTensProd[VList];
    Do[count+=1;
      aux=Product[VList[[j,in[[j]]]],{j,1,Length[VList]}];
      AppendTo[comp,{count,in,aux}];
      in=AddCounter[in,ranges];,{i,1,Size}];
    Return[comp]]

HomogeneousTensorComponents[homogeneoustensorc_]:=
  Module[{in,i,range,ranges,n,Size,comp={},count=0},
    range=homogeneoustensorc[[2,2]];
    in=InitCounter[range];
    n=homogeneoustensorc[[2,1]];
    ranges=Table[n,{i,1,range}];
    Size=n^range;
    Do[count=count+1;
      AppendTo[comp,{count,in,homogeneoustensorc[[1,i]]}];
      in=AddCounter[in,ranges];,{i,1,Size}];
    Return[MatrixForm[comp]]]

RowTensorComponentIndex[n_,r_List]:=Module[{m,i,q=n-1,in},m=Length[r];
    in=Range[m];
    Do[in[[i]]=Mod[q,r[[i]]]+1;q=Quotient[q,r[[i]]],{i,m,2,-1}];
    in[[1]]=q+1;
    Return[in]]

TensorComponentRowIndex[r_List,in_List]:=
  Module[{aux=in[[1]]-1,i,n=Length[in]},Do[aux=aux*r[[i]]+in[[i]]-1,{i,2,n}];
      Return[aux+1]]/;Length[r]\[Equal]Length[in]
(*CHAPTER  3*)
ValidContraction::invalidtensor="`1` is not a tensor.";
ValidContraction::invalidindex="`1` is not a valid index.";
ValidContraction::invalidorder="`1` is not a valid order for contraction.";
ValidContraction::invalidvalency=
    " the valencies of the indices to be contracted must be different.";

ValidContraction[tensor_,p_,q_]:=Module[{aux=False},
    If[IsTensor[tensor],
      If[IntegerQ[p],
        If[IntegerQ[q],
          If[Length[tensor[[2]]]>1,
            If[p<=Length[tensor[[2]]],
              If[q<=Length[tensor[[2]]],
                If[tensor[[2,p]]*tensor[[2,q]]<0,
                  aux=True
                  ,Message[ValidContraction::invalidvalency]]
                ,Message[ValidContraction::invalidindex,q]]
              ,Message[ValidContraction::invalidindex,p]]
            ,Message[ValidContraction::invalidorder,Length[tensor[[2]]]]]
          ,Message[ValidContraction::invalidindex,q]]
        ,Message[ValidContraction::invalidindex,p]]
      ,Message[ValidContraction::invalidtensor,tensor]];
    Return[aux]
    ]

(*Returns a counter of contraction indices list initialized to ones*)

InitContractCounter:=Module[{i},incre={};
    jump=0;
    Do[If[i\[NotEqual]p1&&i\[NotEqual]q1,AppendTo[incre,increments[[i]]],
        jump=jump+increments[[i]]],{i,1,Length[increments]}];
    in=Table[1,{i,1,Length[incre]}];
    Return[in]]

(*Adds one unit to the counter (in list) according to the bases in incre*)

AddContractCounter:=Module[{done=0,j,i},i=rr;
    While[(done\[NotEqual]1&&i\[GreaterEqual]1),
      If[in[[i]]<n,done=1;in[[i]]=in[[i]]+1;counter=counter+incre[[i]];
        Do[in[[j]]=1;counter=counter-(n-1)*incre[[j]],{j,i+1,rr}],i=i-1]];]

TensorContraction[tensor_,p_,q_]:=
  Module[{tensorc,i,j,r,s,order,in,Size,
      neworder},(*tensor has two components:
          the list of components and a list with the dimensions of the \
associated linear spaces*)If[Length[tensor]\[Equal]2,n=Abs[tensor[[2,1]]];
      order=Length[tensor[[2]]];
      If[order>1,If[tensor[[2,p]]*tensor[[2,q]]<0,neworder={};
          
          Do[If[i\[NotEqual]p&&i\[NotEqual]q,
              AppendTo[neworder,tensor[[2,i]]]],{i,1,order}];
          rr=order-2;
          Size=n^rr;
          tensorc=Table[0,{i,1,Size}];
          
          If[p\[LessEqual]order&&q\[LessEqual]order&&p\[GreaterEqual]1&&
              q\[GreaterEqual]1&&p\[NotEqual]q,
            increments=Table[n^(order-i),{i,1,order}];
            p1=p;
            q1=q;
            in=InitContractCounter;
            counter=1;
            
            Do[tensorc[[i]]=
                Sum[tensor[[1,r]],{r,counter,counter+(n-1)*jump,jump}];
              AddContractCounter,{i,1,Size}],
            Print["Erroneous contraction indices"];Return[]],
          Print["The contracted indices must have different valency"];
          Return[]],Print["The order of the tensor must be at least 2"];
        Return[]],
      Print["Tensors must have two components: The list of components and the \
ranges-valency"];Return[]];
    Return[{tensorc,neworder}]]

 TensorContractionMatrix[n1_,order_,p_,q_]:=
  Module[{ChangeMatrix,i,j,r,s,in,Size},rr=order-2;
    Size=n^rr;
    n=n1;
    ChangeMatrix=Table[0,{i,1,Size},{j,1,n^order}];
    If[p\[LessEqual]order&&q\[LessEqual]order&&p\[GreaterEqual]1&&
        q\[GreaterEqual]1&&p\[NotEqual]q,
      increments=Table[n^(order-i),{i,1,order}];
      p1=p;
      q1=q;
      in=InitContractCounter;
      counter=1;
      Do[Do[ChangeMatrix[[i,r]]=1,{r,counter,counter+(n-1)*jump,jump}];
        AddContractCounter,{i,1,Size}],
      Print["Erroneous contraction indices"];];
    Return[ChangeMatrix];] 
  
TensorContractionList[homogeneoustensor_,ContList_]:=
  Module[{aux,i,j,k,k1,cL,aa},cL=ContList;
    aux=homogeneoustensor;
    Do[aux=TensorContraction[aux,cL[[i,1]],cL[[i,2]]];
      Do[aa=Table[0,{k,1,2}];
        Do[If[cL[[i,k1]]<cL[[j,k]],aa[[k]]=aa[[k]]+1],{k,1,2},{k1,1,2}];
        cL[[j,1]]=cL[[j,1]]-aa[[1]];
        cL[[j,2]]=cL[[j,2]]-aa[[2]],{j,i+1,Length[cL]}];,{i,1,
        Length[ContList]}];
    Return[aux]] 

NewDualBasis[OldBasis_,CMatrix_]:=OldBasis.CMatrix ;

NewContComponents[OldContComponents_,CMatrix_]:=
    Inverse[CMatrix].OldContComponents;

NewCovComponents[OldCovComponents_,CMatrix_]:=
    Transpose[CMatrix].OldCovComponents ;

NewGramMatrix[CMatrix_]:=Transpose[CMatrix].CMatrix;

(*CHAPTER  4*)
BaseTensorChange[CmatricesList_,CovariancesList_]:=
  Module[{i,CC={}},n=Length[CmatricesList];
    If[n\[NotEqual]Length[CovariancesList],
      Print["Data matrices must have the same length"];Return[],
      Do[If[CovariancesList[[i]]\[Equal]1,
          AppendTo[CC,Inverse[CmatricesList[[i]]]],
          AppendTo[CC,Transpose[CmatricesList[[i]]]]],{i,1,n}];
      Return[ListMatrixTensorProduct[CC]]]]

(*InverseBaseTensorChange[CmatricesList_,CovariancesList_]:=
    Inverse[\
          BaseTensorChange[CmatricesList,CovariancesList]]*) 


InverseBaseTensorChange[CmatricesList_,CovariancesList_]:=
  Module[{i,CC={}},n=Length[CmatricesList];
    If[n\[NotEqual]Length[CovariancesList],
      Print["Data matrices must have the same length"];Return[],
      Do[If[CovariancesList[[i]]\[Equal]1,AppendTo[CC,CmatricesList[[i]]],
          AppendTo[CC,Transpose[Inverse[CmatricesList[[i]]]]]],{i,1,n}];
      Return[ListMatrixTensorProduct[CC]]]]

(*CHAPTER  5*)
 ListEigenValuesVectors[MList_]:=
  Module[{eigenval={},eigenvec={},i},
    Do[AppendTo[eigenval,N[Eigenvalues[MList[[i]]]]];
      AppendTo[eigenvec,N[Eigenvectors[MList[[i]]]]],{i,1,Length[MList]}];
    Return[Chop[{eigenval,eigenvec}]]] 

TensorProductEigenValuesVectors[MList_]:=
  Module[{eigenval,eigenvec},aux=ListEigenValuesVectors[MList];
    eigenval=VectorListTensorProduct[aux[[1]]][[1]];
    eigenvec=ListMatrixTensorProduct[aux[[2]]];
    Return[{eigenval,eigenvec}]]
(*CHAPTER  6*)
 TemplateSymmetricTensor[order_,n_]:=
  Module[{ranges,i},ranges=Table[n,{i,1,order}];
    insim=InitCounter[order];
    Size=Binomial[n+order-1,order];
    Do[Print[insim];
      insim=AddSymmetricCounter[insim,ranges];,{i,1,Size}];] 

AddSymmetricCounter[insim_List,r_List]:=
  Module[{Done=False,rr=Length[r],in1sim=insim,j},i=rr;
    While[!Done&&i\[GreaterEqual]1,
      If[in1sim[[i]]<r[[i]],Done=True;in1sim[[i]]+=1;
          Do[in1sim[[j]]=in1sim[[i]],{j,i+1,rr}],i-=1];];
    Return[in1sim]]

 SymmetrizeTensorCompact[tensor_]:=
  Module[{i,m,order=Length[tensor[[2]]],ranges=Abs[tensor[[2]]],
      n=Abs[tensor[[2,1]]],aux,Size,j,k,in,insim},insim=InitCounter[order];
    Size=Binomial[n+order-1,order];
    aux=Table[0,{Size}];
    Do[in=Permutations[insim];
      m=Length[in];
      Do[j=TensorComponentRowIndex[ranges,in[[k]]];
        aux[[i]]+=tensor[[1,j]];,{k,1,m}];
      aux[[i]]=aux[[i]]/m;
      insim=AddSymmetricCounter[insim,ranges];,{i,1,Size}];
    Return[{aux,tensor[[2]]}]]
  
 SymmetrizeTensor[tensor_]:=
  Module[{i,m,order=Length[tensor[[2]]],ranges=Abs[tensor[[2]]],
      n=Abs[tensor[[2,1]]],aux,aux1,Size,Size1,j,k,in,insim},
    insim=InitCounter[order];
    Size=Binomial[n+order-1,order];
    Size1=Length[tensor[[1]]];
    aux=Table[0,{Size}];
    aux1=Table[0,{Size1}];
    Do[in=Permutations[insim];
      m=Length[in];
      Do[j=TensorComponentRowIndex[ranges,in[[k]]];
        aux[[i]]+=tensor[[1,j]];,{k,1,m}];
      aux[[i]]=aux[[i]]/m;
      Do[j=TensorComponentRowIndex[ranges,in[[k]]];
        aux1[[j]]=aux[[i]];,{k,1,m}];
      insim=AddSymmetricCounter[insim,ranges];,{i,1,Size}];
    Return[{aux1,tensor[[2]]}]]
  
SymmetricTensorToTensor[symmetrictensor_]:=
  Module[{order=Length[symmetrictensor[[2]]],n=Abs[symmetrictensor[[2,1]]],
      Sizesym,Size,insim,ranges=symmetrictensor[[2]],ranges1,comp,
      comp1=symmetrictensor[[1]],isym=1,in,m,k,j,aux},
    ranges1=Abs[ranges];
    Size=n^order;
    comp=Table[0,{i,1,Size}];
    Sizesym=Binomial[n+order-1,order];
    insim=InitCounter[order];
    Do[
      in=Permutations[insim];
      m=Length[in];
      Do[j=TensorComponentRowIndex[ranges1,in[[k]]];
        comp[[j]]=comp1[[isym]];,{k,1,m}];
      insim=AddSymmetricCounter[insim,ranges1];
      isym++;
      ,{i,1,Sizesym}];
    Return[{comp,ranges}]
    ]

IsSymmetricTensor::invalidsymmetrictensor="`1` is not a symmetric tensor.";

IsSymmetricTensor[tensor_]:=
  Module[{aux=False,aux1,m,i,order=Length[tensor[[2]]]},
    If[ListQ[tensor],
      If[Length[tensor]\[Equal]2,
        If[ListQ[tensor[[1]]],
          If[ListQ[tensor[[2]]],
            
            If[Length[tensor[[1]]]==
                Binomial[Abs[tensor[[2,1]]]+order-1,order],
              aux=True;
              If[aux,
                i=1;aux1=True;
                m=Length[tensor[[2]]];
                While[i<=m && aux1,
                  
                  If[Length[tensor[[2,i]]]>0 && IntegerQ[tensor[[2,i]]],
                    aux1=False];
                  i++;
                  ];
                aux=aux&&aux1;
                ]
              ]
            ]
          ]
        ]
      ];
    If[!aux,Message[IsSymmetricTensor::invalidsymmetrictensor,tensor]];
    Return[aux]
    ]

SymmetricTensorStrictComponents[tensor_]:=
  Module[{in,i,range,ranges,ranges1,n=tensor[[2,1]],Size,comp={},count=0,
        order=Length[tensor[[2]]]},
      If[tensor[[2]]\[Equal]{},Return[{1,{0},tensor[[1]]}],
        range=Length[tensor[[2]]],aa];
      in=InitCounter[range];
      ranges=Abs[tensor[[2]]];
      ranges1=tensor[[2]];
      Size=Binomial[n+order-1,order];
      Do[count=count+1;
        aa={in,in};
        If[ranges1[[1]]>0,aa[[1,1]]=StringJoin["(",ToString[aa[[1,1]]]],
          aa[[2,1]]=StringJoin["(",ToString[aa[[1,1]]]]];
        If[ranges1[[range]]>0,
          aa[[1,range]]=StringJoin[ToString[aa[[1,range]]],")"],
          aa[[2,range]]=StringJoin[ToString[aa[[1,range]]],")"]];
        Do[If[ranges1[[j]]>0,aa[[2,j]]="o",aa[[1,j]]="o"],{j,1,range}];
        AppendTo[comp,{count,MatrixForm[aa],tensor[[1,i]]}];
        in=AddSymmetricCounter[in,ranges],{i,1,Size}];
      Return[MatrixForm[comp]]]/;IsSymmetricTensor[tensor]
    
(*CHAPTER  7*)

AddAntiSymmetricCounter[insim_List,r_List]:=
  Module[{Done=False,rr=Length[r],in1sim=insim,j},i=rr;
    While[!Done&&i\[GreaterEqual]1,
      If[in1sim[[i]]<r[[i]]+i-rr,Done=True;in1sim[[i]]+=1;
          Do[in1sim[[j]]=in1sim[[i]]+1,{j,i+1,rr}],i-=1];
      ];
    Return[in1sim]]

InitAntiSymmetricCounter[order_,n_]:=
  If[order\[LessEqual]n,Table[i,{i,1,order}],0]

 TemplateAntiSymmetricTensor[order_,n_]:=
  Module[{ranges,i},ranges=Table[n,{i,1,order}];
    If[order\[LessEqual]n,
      insim=InitAntiSymmetricCounter[order,n];
      Size=Binomial[n,order];
      Do[Print[insim];
        insim=AddAntiSymmetricCounter[insim,ranges];,{i,1,Size}],
      Print["Is order less than n"]]
    ] 
  
AntiSymmetrizeTensorCompact[tensor_]:=
  Module[{i,m,order=Length[tensor[[2]]],ranges=Abs[tensor[[2]]],
      n=Abs[tensor[[2,1]]],aux,Size,j,k,in,insim},
    insim=InitAntiSymmetricCounter[order,n];
    Size=Binomial[n,order];
    aux=Table[0,{Size}];
    Do[in=Permutations[insim];
      m=Length[in];
      Do[j=TensorComponentRowIndex[ranges,in[[k]]];
        aux[[i]]+=Signature[in[[k]]]*tensor[[1,j]];,{k,1,m}];
      aux[[i]]=aux[[i]]/m;
      insim=AddAntiSymmetricCounter[insim,ranges];,{i,1,Size}];
    Return[{aux,tensor[[2]]}]] 

AntiSymmetrizeTensor[tensor_]:=
  Module[{i,m,order=Length[tensor[[2]]],ranges=Abs[tensor[[2]]],
      n=Abs[tensor[[2,1]]],aux,aux1,Size,Size1,j,k,in,insim},
    insim=InitCounter[order];
    Size=Binomial[n+order-1,order];
    Size1=Length[tensor[[1]]];
    aux=Table[0,{Size}];
    aux1=Table[0,{Size1}];
    Do[in=Permutations[insim];
      m=Length[in];
      Do[j=TensorComponentRowIndex[ranges,in[[k]]];
        aux[[i]]+=Signature[in[[k]]]*tensor[[1,j]];,{k,1,m}];
      aux[[i]]=aux[[i]]/m;
      Do[j=TensorComponentRowIndex[ranges,in[[k]]];
        aux1[[j]]=aux[[i]];,{k,1,m}];
      insim=AddSymmetricCounter[insim,ranges];,{i,1,Size}];
    Return[{aux1,tensor[[2]]}]]

AntySymmetricTensorToTensor[antisymmetrictensor_]:=
  Module[{order=Length[antisymmetrictensor[[2]]],
      n=Abs[antisymmetrictensor[[2,1]]],Sizeanti,Size,insimanti,
      ranges=antisymmetrictensor[[2]],ranges1,comp,
      comp1=antisymmetrictensor[[1]],ianti=1,in,m,k,j,aux},
    ranges1=Abs[ranges];
    Size=n^order;
    comp=Table[0,{i,1,Size}];
    Sizeanti=Binomial[n,order];
    insimanti=InitAntiSymmetricCounter[order,n];
    Do[
      in=Permutations[insimanti];
      m=Length[in];
      Do[j=TensorComponentRowIndex[ranges1,in[[k]]];
        comp[[j]]=comp1[[ianti]];,{k,1,m}];
      insimanti=AddAntiSymmetricCounter[insim,ranges1];
      ianti++;
      ,{i,1,Sizeanti}];
    Return[{comp,ranges}]
    ]

IsAntiSymmetricTensor::invalidantisymmetrictensor=
    "`1` is not an antisymmetric tensor.";

IsAntiSymmetricTensor[tensor_]:=
  Module[{aux=False,aux1,m,i,order=Length[tensor[[2]]]},
    If[ListQ[tensor],
      If[Length[tensor]\[Equal]2,
        If[ListQ[tensor[[1]]],
          If[ListQ[tensor[[2]]],
            If[Length[tensor[[1]]]==Binomial[Abs[tensor[[2,1]]],order],
              aux=True;
              If[aux,
                i=1;aux1=True;
                m=Length[tensor[[2]]];
                While[i<=m && aux1,
                  
                  If[Length[tensor[[2,i]]]>0 && IntegerQ[tensor[[2,i]]],
                    aux1=False];
                  i++;
                  ];
                aux=aux&&aux1;
                ]
              ]
            ]
          ]
        ]
      ];
    If[!aux,
      Message[IsAntiSymmetricTensor::invalidantisymmetrictensor,tensor]];
    Return[aux]
    ]


AntiSymmetricTensorStrictComponents[tensor_]:=
  Module[{in,i,range,ranges,ranges1,n=tensor[[2,1]],Size,comp={},count=0,
        order=Length[tensor[[2]]]},
      If[tensor[[2]]\[Equal]{},Return[{1,{0},tensor[[1]]}],
        range=Length[tensor[[2]]],aa];
      in=InitAntiSymmetricCounter[order,n];
      ranges=Abs[tensor[[2]]];
      ranges1=tensor[[2]];
      Size=Binomial[n,order];
      Do[count=count+1;
        aa={in,in};
        If[ranges1[[1]]>0,aa[[1,1]]=StringJoin["(",ToString[aa[[1,1]]]],
          aa[[2,1]]=StringJoin["(",ToString[aa[[1,1]]]]];
        If[ranges1[[range]]>0,
          aa[[1,range]]=StringJoin[ToString[aa[[1,range]]],")"],
          aa[[2,range]]=StringJoin[ToString[aa[[1,range]]],")"]];
        Do[If[ranges1[[j]]>0,aa[[2,j]]="o",aa[[1,j]]="o"],{j,1,range}];
        AppendTo[comp,{count,MatrixForm[aa],tensor[[1,i]]}];
        in=AddAntiSymmetricCounter[in,ranges],{i,1,Size}];
      Return[MatrixForm[comp]]]/;IsAntiSymmetricTensor[tensor]
(*CHAPTER  8*)
DeltaK[i_,j_]:=If[i\[Equal]j,1,0] 

GeneralizedKroneckerDeltaR[upperindices_,lowerindices_]:=
  Det[Table[
      DeltaK[upperindices[[i]],lowerindices[[j]]],{i,1,
        Length[upperindices]},{j,1,Length[lowerindices]}]] 

GeneralizedKroneckerDelta[upperindices_,lowerindices_]:=
    LeviCivita[upperindices]*LeviCivita[lowerindices] ;

LeviCivita[indices_]:=Signature[indices]

(* LeviCivita1[indices_]:=
    Module[{i,j},
      If[Length[Union[indices]]\[Equal]Length[indices],
        Return[Product[
            Product[Sign[indices[[i+j]]-indices[[j]]],{i,1,
                Length[indices]-j}],{j,1,Length[indices]-1}]],Return[0]]] *)

 CreateLeviCivita[order_,n_,var_]:=
  Module[{in,i,Size,comp={}},in=InitCounter[order];
    ranges=n*in;
    Size=n^order;
    Do[AppendTo[comp,LeviCivita[in]];
      in=AddCounter[in,ranges];,{i,1,Size}];
    Return[{comp,var*ranges}]] 
  
TensorProduct[tensor1_,tensor2_]:=
  Module[{order1,order2,ranges,Size,ranges0,comp={},ranges1,ranges2,in,
      in1,\
        in2,i},order1=Length[tensor1[[2]]];
    order2=Length[tensor2[[2]]];
    ranges=Join[tensor1[[2]],tensor2[[2]]];
    ranges0=Abs[ranges];
    ranges1=Take[ranges0,{1,order1}];
    ranges2=Take[ranges0,{order1+1,order1+order2}];
    Size=Product[ranges0[[i]],{i,1,Length[ranges0]}];
    in=InitCounter[order1+order2];
    Do[in1=Take[in,{1,order1}];
      in2=Take[in,{order1+1,order1+order2}];
      AppendTo[comp,
        tensor1[[1,TensorComponentRowIndex[ranges1,in1]]]*
          tensor2[[1,TensorComponentRowIndex[ranges2,in2]]]];
      in=AddCounter[in,ranges0];,{i,1,Size}];
    Return[{comp,ranges}]] 

ListTensorProduct[TList_]:=
  Module[{n=Length[TList],aux=TList[[1]],i},
    If[n\[NotEqual]1,Do[aux=TensorProduct[aux,TList[[i]]],{i,2,n}];];
    Return[aux]] 

AdjTensor[tensor_]:=
  Module[{q=Length[tensor[[2]]],aux,var=Sign[tensor[[2,1]]],i,aux1},
    n=Abs[tensor[[2,1]]];
    aux1=CreateLeviCivita[n,n,-var];
    aux=TensorProduct[aux1,tensor];
    aux[[1]]=aux[[1]]/q!;
    Return[TensorContractionList[aux,Table[{i,n+i},{i,1,q}]]]]
(*CHAPTER  9*)
(*CHAPTER 9*)
(*CHAPTER  11*)
 IndexExchange[tensor_,changes_,Gram_]:=
  Module[{ranges,order,n,i,aux={},aux1},ranges=tensor[[2]];
    order=Length[ranges];
    n=ranges[[1]];
    Do[If[changes[[i]]\[Equal]1,aux1=IdentityMatrix[n],
        If[changes[[i]]*ranges[[i]]<0,aux1=Inverse[Gram],aux1=Gram]];
      AppendTo[aux,aux1],{i,1,order}];
    aux1=ListMatrixTensorProduct[aux];
    Return[aux1]]

 InitStrictCounter[order_]:=Table[i,{i,1,order}] 
  
AddStrictCounter[insim_List,r_List]:=
  Module[{Done=False,rr=Length[r],in1sim=insim,j},i=rr;
    While[!Done&&i\[GreaterEqual]1,
      If[in1sim[[i]]<r[[i]]-rr+i,Done=True;in1sim[[i]]+=1;
          Do[in1sim[[j]]=in1sim[[j-1]]+1,{j,i+1,rr}],i-=1];];
    Return[in1sim]]

 TemplateStrictTensor[order_,n_]:=
  Module[{ranges,i,in},ranges=Table[n,{i,1,order}];
      in=InitStrictCounter[order];
      Size=Binomial[n,order];
      Do[Print[in];
        in=AddStrictCounter[in,ranges],{i,1,Size}];]/;n\[GreaterEqual]order 
      
ExteriorProduct[vectors_List]:=
  Module[{ranges,Size,in,i,j,r,n,comp={},aux},n=Length[vectors];
    r=Length[vectors[[1]]];
    ranges=Table[n,{i,1,r}];
    Size=Binomial[n,r];
    in=InitStrictCounter[r];
    Do[aux=Table[vectors[[in[[k]],j]],{k,1,r},{j,1,r}];
      aux=Det[aux];
      AppendTo[comp,aux];
      in=AddStrictCounter[in,ranges],{i,1,Size}];
    Return[{comp,{n,r}}]] 

EuclideanTensorScalarProduct[tensor1_,tensor2_,Gram_]:=
  Module[{ranges1,ranges2,order,n,i,aux={}},n=tensor1[[2,1]];
    ranges1=tensor1[[2]];
    ranges2=tensor2[[2]];
    order=Length[ranges1];
    Do[If[ranges1[[i]]>0&&ranges2[[i]]>0,aux1=Gram,
        If[ranges1[[i]]<0&&ranges2[[i]]<0,aux1=Inverse[Gram],
          aux1=IdentityMatrix[n]]];
      AppendTo[aux,aux1],{i,1,order}];
    aux1=tensor1[[1]].ListMatrixTensorProduct[aux].Transpose[{tensor2[[1]]}];
    Return[aux1[[1]]]] 

EuclideanTensorScalarProductMatrix[tensor1_,tensor2_,Gram_]:=
  Module[{ranges1,ranges2,order,n,i,aux={}},n=tensor1[[2,1]];
    ranges1=tensor1[[2]];
    ranges2=tensor2[[2]];
    order=Length[ranges1];
    Do[If[ranges1[[i]]>0&&ranges2[[i]]>0,aux1=Gram,
        If[ranges1[[i]]<0&&ranges2[[i]]<0,aux1=Inverse[Gram],
          aux1=IdentityMatrix[n]]];
      AppendTo[aux,aux1],{i,1,order}];
    aux1=ListMatrixTensorProduct[aux];
    Return[aux1[[1]]]] 

TensorModulus[tensor1_,Gram_]:=
  Sqrt[EuclideanTensorScalarProduct[tensor1,tensor1,Gram]] 

TensorAngle[tensor1_,tensor2_,Gram_]:=
  EuclideanTensorScalarProduct[tensor1,tensor2,
      Gram]/(TensorModulus[tensor1,Gram]*TensorModulus[tensor2,Gram])
(*CHAPTER  12*)
EuclideanAdjTensor[tensor_,Gram_]:=
  Module[{q=Length[tensor[[2]]],aux,var=Sign[tensor[[2,1]]],i,aux1,aux2},
    n=Abs[tensor[[2,1]]];
    aux1=CreateLeviCivita[n,n,-var];
    aux=TensorProduct[aux1,tensor];
    aux2=Sqrt[Det[Gram]];
    If[tensor[[2,1]]<0,aux2=1/aux2];
    aux[[1]]=aux2*aux[[1]]/q!;
    aux1=TensorContractionList[aux,Table[{i,n+i},{i,1,q}]];
    Return[aux1]]
(*CHAPTER  14*)
ProjectionTensor[Direction_]:=Direction.Transpose[Direction] 
ProjectionComplementTensor[
    Direction_]:=(IdentityMatrix[Length[Direction]]-
      Direction.Transpose[Direction]) 

MomentTensor[
    r_]:={{0,-r[[3,1]],r[[2,1]]},{r[[3,1]],0,-r[[1,1]]},{-r[[2,1]],r[[1,1]],
      0}} 

AxisMomentTensor[r_,axis_]:=ProjectionTensor[axis].MomentTensor[r] 
  
RotationTensor[angle_,axis_]:=
  Cos[angle]*IdentityMatrix[Length[axis]]+(1-Cos[angle])*
      ProjectionTensor[axis]+Sin[angle]*MomentTensor[axis] 
      
ReflexionTensor[axis_]:=IdentityMatrix[Length[axis]]-2*ProjectionTensor[axis]
      
AffinityTensor[A_,L_]:=
  Append[Transpose[Append[Transpose[A],L]],
    Append[Table[0,{i,1,Length[L]}],1]]

HomoteticTensor[center_,k_]:=
  Append[Transpose[
      Append[Transpose[k*IdentityMatrix[Length[center]]],
        Table[(1-k)*center[[i]],{i,1,Length[center]}]]],
    Append[Table[0,{i,1,Length[center]}],1]] 

TranslationTensor[t_]:=
  Append[Transpose[
      Append[Transpose[IdentityMatrix[Length[t]]],
        Table[t[[i]],{i,1,Length[t]}]]],Append[Table[0,{i,1,Length[t]}],1]] 

RotationTensor2D[angle_]:={{Cos[angle],-Sin[angle]},{Sin[angle],Cos[angle]}}

 CentralSymmetryTensor[center_]:=HomoteticTensor[center,-1] 
  
AxialSymmetryTensor2D[A_,B_,
    C1_]:={{-A^2+B^2,-2*A*B,-2*A*C1},{-2*A*B,A^2-B^2,-2*B*C1},{0,0,
        A^2+B^2}}/(A*A+B*B)

TensorContractionHomogeneous[homogeneoustensor_,p_,q_]:=
  Module[{homogeneoustensorc,i,j,r,s,order,in,Size},
    (* Homogeneoustensor has two components: 
          the list of components and a list with the dimension of the linear \
space and the tensor order *)
    If[Length[homogeneoustensor]==2,
      n=homogeneoustensor[[2,1]];
      order=homogeneoustensor[[2,2]];
      rr=order-2;
      Size=n^rr;
      homogeneoustensorc=Table[0,{i,1,Size}];
      If[p<=order && q<=order && p>=1 &&q>=1 && p\[NotEqual]q,
        increments=Table[n^(order-i),{i,1,order}];
        p1=p;
        q1=q;
        in=InitContractCounter;
        counter=1;
        Do[
          homogeneoustensorc[[i]]=
            Sum[homogeneoustensor[[1,r]],{r,counter,counter+(n-1)*jump,
                jump}];
          AddContractCounter;
          ,{i,1,Size}],
        Print["Erroneous contraction indices"];
        ],
      Print["homogeneous tensors must have two components"]
      ];
    Return[{homogeneoustensorc,{n,rr}}];
]
        
TensorContractionListHomogeneous[homogeneoustensor_,ContList_]:=
  Module[{aux,i,j,k,k1,cL,aa},
    aa=Table[0,{k,1,2}];
    cL=ContList;
    aux=homogeneoustensor;
    Do[aux=TensorContraction[aux,cL[[i,1]],cL[[i,2]]];
      Do[Do[If[cL[[i,k1]]<cL[[j,k]],aa[[k]]=aa[[k]]+1],{k,1,2},{k1,1,2}];
        cL[[j,1]]=cL[[j,1]]-aa[[1]];
        cL[[j,2]]=cL[[j,2]]-aa[[2]],{j,i+1,Length[cL]}];
      ,{i,1,Length[ContList]}];
    Return[aux]
]

VectorTensorProduct[VList_]:=Module[{in,i,ranges,Size,comp={},count=0,aux},
    in=InitCounter[Length[VList]];
    ranges=CalculateRangesTensProd[VList];
    Size=CalculateSizeTensProd[VList];
    Do[count+=1;
      aux=Product[VList[[j,in[[j]]]],{j,1,Length[VList]}];
      AppendTo[comp,{aux}];
      in=AddCounter[in,ranges];
      ,{i,1,Size}];
    Return[{comp,ranges}]
    ]
    
TensorProductComponents[VList_]:=
  Module[{in,i,ranges,Size,comp={},count=0,aux},
    in=InitCounter[Length[VList]];
    ranges=CalculateRangesTensProd[VList];
    Size=CalculateSizeTensProd[VList];
    Do[count+=1;
      aux=Product[VList[[j,in[[j]]]],{j,1,Length[VList]}];
      AppendTo[comp,{count,in,aux}];
      in=AddCounter[in,ranges];
      ,{i,1,Size}];
    Return[comp]
    ]

HomogeneousTensorComponents[homogeneoustensorc_]:=
  Module[{in,i,range,ranges, n,Size,comp={},count=0},
    range=homogeneoustensorc[[2,2]];
    in=InitCounter[range];
    n=homogeneoustensorc[[2,1]];
    ranges=Table[n,{i,1,range}];
    Size=n^range;
    Do[count=count+1;
      AppendTo[comp,{count,in,homogeneoustensorc[[1,i]]}];
      in=AddCounter[in,ranges];
      ,{i,1,Size}];
    Return[MatrixForm[comp]]
    ]

TensorMatrixForm[tensor_,nrowindices_]:=
  Module[{rows=Product[Abs[tensor[[2,j]]],{j,1,nrowindices}],
      columns=Product[Abs[tensor[[2,j]]],{j,nrowindices+1,Length[tensor[[2]]]}]},
    A=Table[Take[tensor[[1]],{(i-1)*columns+1,i*columns}],{i,1,rows}];
    Return[A]
    ]
 
 Permute[in_, perm_]:=Table[in[[perm[[i]]]],{i,1,Length[perm]}];

Permutation[tensor_,perm_]:=
  Module[{tensor1=tensor[[1]],Size,i,in,ranges=tensor[[2]],m,ranges1,
      order=Length[tensor[[2]]]},
    ranges1=Permute[ranges,perm];
    Size=Length[tensor1];
    in=InitCounter[order];
    in1=in;
    Do[m=TensorComponentRowIndex[ranges1,in1];
      tensor1[[m]]=tensor[[1,i]];
      in=AddCounter[in,ranges];
      in1=Permute[in,perm];
      ,{i,1,Size}];
    Return[{tensor1,ranges1}]
    ]  
    
End[]

EndPackage[]

























