Ideas to introduce indexed objects (possibly representing tensors)

Hi there!
I'm currently writing a system of indexed n-dimensional arrays to be used for the Physics module.

Before starting, I went through other proposals, and I saw suggestions about using GA or diffgeom for the physics module.

GA or diffgeom would be great ideas, but actually their implementation is very complicated and would require a lot of work before achieving any result.

I thought that the best way to write modern physics would be to do it the same way all physics textbooks do it: indexed multidimensional arrays with contractions.

This approach has many advantages:
provide a framework for multidimensional arrays with covariant/contravariant parts [for now basis-dependent, but maybe in the future this work could be integrated with the diffgeom module].
provide the same approach to physical formulae as usual textbooks and papers.
Example: you could be reading a paper about some non-standard physics, and be willing to do the calculations by yourself. It would be an enormous amount of work to start from the diffgeom/GA modules! Better something simpler.
using subclassing to define how these multidimensional arrays transform under various physical transformations (Lorentz group, isospin, color symmetry).
I will provide a Pull Request as soon as I can [not this week]. For now I have implemented the following:
class _PartlyContractedSymbol
class IndexedSymbol
class RelativitySymbol(IndexedSymbol)
indexes are either integers or Sympy's Symbols, square brackets means covariant, round brackets contravariant.

Examples:

minkowsky_metric = [
    (1, 0, 0, 0),
    (0,-1, 0, 0),
    (0, 0,-1, 0),
    (0, 0, 0,-1),
]

>>> CC = RelativitySymbol("CC", ['lorentz', 'lorentz'], minkowsky_metric)
>>> assert CC.get_matrix() == Matrix(minkowsky_metric)


In this one, U has three indices, 1st and 2nd depend on Minkowski metric, i.e. Matrix(1,-1,-1,-1), while the second is Euclidean in 4 dimensions:
>>> assert U[1, 2, 1] == 2 
>>> assert U[1](2)[1] == 2 # 2nd index obeys Euclidean metric
>>> assert U[1, 2](1) == -2 # 3rd index obeys Minkowsky metric.

I think I'm going to provide a pull request next week. Let me know what you think.

By the way, I also started a FourVector and a "contraction" system for symbolic indices:

>>> momentum = FourMomentum('momentum', [E, 0, 0, pz])
>>> assert contract(momentum[mu]*momentum(mu)) == E**2 - pz**2
>>> assert momentum**2 == E**2 - pz**2
 
and a Lorentz boost using multiplication with Matrix.exp() of the basis of the boost generators:
>>> px_boosted_momentum = momentum.boost([px, 0, 0])