## FileName: Riemann.mpl
## Contents: maple 18 program to calculate the curvature tensor of a metric
## CreatedBy: Kodama Hideo 
## CreationDate: 2004/2/16
## 2020/9/21 -10/2 : the program is modifiled to the package type
## 2020/10/3: pds version 1
## 2022/11/24 : minor bug fix => version 1.1
## LastUpdate: 2022/11/24

######>TOP
#####>packages
with(LinearAlgebra):

#####>Joblist

################################################################
######>Data to be given
#####>varlist
#% Coordinates
#phi:='phi': theta:='theta': psi:='psi':
#t:='t': r:='r': 
#x:='x': y:='y':
#varlist:=[t,r,theta,phi];
#dim:=nops(varlist);
#####>metric
#% Parameters
#a:='a':
#% metric
#ds2:=-(1-a/4/r)^2/(1+a/4/r)^2*dt^2+(1+a/4/r)^4*(dr^2+r^2*(dtheta^2+sin(theta)^2*dphi^2)):
#%
###############################################################

######>Global variables
_u:='_u': _d:='_d':  #% index position values
FB:=table(antisymmetric):

######>conventions
#% T:: tensor table =table([(1,2,..)=..., ..., 
#%         "indextype"=[_u,_d,...], "dim"=... ])

######>Curvature calculation
#####>proc:df
#% df(eq, varlist)
df:=proc()
 local eq,varlist, dim,
    var, dvarlist,i:
    
  if _npassed<2 then
    error("Usage: df(eq,varlist");
  end if:
  eq:=_passed[1]:
  varlist:=_passed[2]:
  if not is(varlist,list) then
    error("Usage: df(eq,varlist::list)");
  end if:
  dim:=nops(varlist):
  
  var:='var':
  dvarlist:=map(var -> d||var, varlist):
  add(diff(eq,varlist[i])*dvarlist[i],i=1..dim);
end proc:


#####>proc:Riemann
#% Procedure to calculate the connection and the Riemann tensor.
#% The output is used to calculate derived quantities such as 
#% the Ricci tensor, covariant derivatives and so on.
##% Riemann(varlist,ds2[,sw]) => Gdata = geometry deta table (ds2, Ch,Riem, Ric, ET, Weyl)
#% sw=0 => no simplification
#% sw=1 => no simplification with stage messages
#% sw=2 -> simplify Riem with stage messages 
Riemann:=proc()
local varlist, ds2, sw,
      dvarlist,  xx,dxx, G, GI, 
      var,ii, jj, jj1,jj2, jj3, kk,ll,mm,
      dim,Ch, Riem13, Riem04,Riem22, Ric02, Ric11, RS, ET02, ET11,
      Weyl04, Gdata;

##% read the input data
if _npassed<2 then
  error("Usage: Riemann(varlist, ds2 [,sw]");
else
  varlist:=_passed[1]:
  ds2:=_passed[2]:
end if:
sw:=0:
if _npassed>=3 then
  sw:=_passed[3]:
end if:
dim:=nops(varlist):

##% dvarlist
var:='var':
dvarlist:=map(var -> d||var, varlist):

#% general coordinates
xx:=Array(1..dim):
dxx:=Array(1..dim):
ii:='ii':
for ii from 1 to dim do
  xx[ii]:=varlist[ii]:
  dxx[ii]:=dvarlist[ii]:
end do:
####>G[i,j] & GI[i,j]
#%
G:=Matrix(1..dim,1..dim,shape=symmetric):
ii:='ii':
for ii from 1 to dim do
  G[ii,ii]:=factor(coeff(ds2,dxx[ii]^2)):
end do:
ii:='ii': jj:='jj':
for ii from 1 to dim do
  for jj from ii+1 to dim do
     G[ii,jj]:=1/2*factor(coeff(coeff(ds2,dxx[ii]),dxx[jj])):
  end do:
end do:
GI:=G^(-1):

####>Ch, Riem13, Riem04, Riem22
Ch:=Array(1..dim,1..dim,1..dim): # Ch[i,j,k]=Gamma^i_{jk}
Riem13:=Array(1..dim,1..dim,1..dim,1..dim):
Riem04:=Array(1..dim,1..dim,1..dim,1..dim):
Riem22:=Array(1..dim,1..dim,1..dim,1..dim):

####> Christoffel symbols
if sw>0 then
  print("Calculating the  Christoffel symbols");
end if:
#%
ii:='ii': 
for ii from 1 by 1 to dim do 
 jj1:='jj1':
 for jj1 from 1 by 1 to dim do 
  jj2:='jj2':
  for jj2 from 1 by 1 to dim do 
    jj3:='jj3':
    Ch[ii,jj1,jj2]:=add(1/2*GI[ii,jj3]*(diff(G[jj1,jj3],xx[jj2])\ 
                   +diff(G[jj2,jj3],xx[jj1])\ 
                   -diff(G[jj1,jj2],xx[jj3])),jj3=1..dim):
    Ch[ii,jj1,jj2]:=simplify(Ch[ii,jj1,jj2]):
  end do:
 end do:
end do;

####>Riemann curvature
if sw>0 then
  print("Calculating Curvatures");
end if:
###>Riem13 & Riem04& Riem22
if sw>0 then
  print("Riem13[*,*,*,*], Riem04[*,*,*,*], Riem22[*,*,*,*]");
end if:
#%
 kk:='kk': for kk from 1 to dim do
  ll:='ll': for ll from kk to dim do
    ii:='ii': for ii from 1 to dim do
      jj:='jj': for jj from 1 to dim do 
        mm:='mm':       
        Riem13[ii,jj,kk,ll]:=diff(Ch[ii,jj,ll],xx[kk])-diff(Ch[ii,jj,kk],xx[ll]) 
       +add(Ch[ii,kk,mm]*Ch[mm,jj,ll]-Ch[ii,ll,mm]*Ch[mm,jj,kk],mm=1..dim):
        if sw>1 then
          Riem13[ii,jj,kk,ll]:=simplify(Riem13[ii,jj,kk,ll]);
        end if:
        Riem13[ii,jj,ll,kk]:=-Riem13[ii,jj,kk,ll]:
      end do:
    end do:
    ii:='ii': for ii from 1 to dim do
      jj:='jj': for jj from 1 to dim do 
        mm:='mm':
        Riem04[ii,jj,kk,ll]:=add(G[ii,mm]*Riem13[mm,jj,kk,ll],mm=1..dim):
       mm:="mm":
        Riem22[ii,jj,kk,ll]:=add(GI[jj,mm]*Riem13[ii,mm,kk,ll],mm=1..dim):
        if sw>1 then
           Riem04[ii,jj,kk,ll]:=simplify(Riem04[ii,jj,kk,ll]):
           Riem22[ii,jj,kk,ll]:=simplify(Riem22[ii,jj,kk,ll]):
        end if:
        Riem04[ii,jj,ll,kk]:=-Riem04[ii,jj,kk,ll]:
        Riem22[ii,jj,ll,kk]:=-Riem22[ii,jj,kk,ll]:
    end do:
  end do:
 end do:
end do:


####>Ricci curvatures
if sw>0 then
  print("Ricci curvatures: Ric02[*,*], Ric11[*,*], RS");
end if:
###>Ric02
Ric02:=Matrix(1..dim,1..dim):
ii:='ii': 
for ii from 1 by 1 to dim do
 jj:='jj':
 for jj from ii by 1 to dim do
  kk:='kk':
  Ric02[ii,jj]:=add(Riem13[kk,ii,kk,jj],kk=1..dim):
  if sw>1 then
    Ric02[ii,jj]:=simplify(Ric02[ii,jj]):
  end if:
  if ii<jj then
    Ric02[jj,ii]:=Ric02[ii,jj]:
  end if:
 end do:
end do:

###>Ric11
Ric11:=Matrix(1..dim,1..dim): # Ric11[i,j]=R^i{}_j
ii:='ii': 
for ii from 1 by 1 to dim do
 jj:='jj':
 for jj from 1 by 1 to dim do
  kk:='kk':
  Ric11[ii,jj]:=add(GI[ii,kk]*Ric02[kk,jj],kk=1..dim):
  if sw>1 then
    Ric11[ii,jj]:=simplify(Ric11[ii,jj]):
  end if:
 end do:
end do:

###>Rs
ii:='ii': 
RS:=add(Ric11[ii,ii],ii=1..dim):
if sw>1 then
  RS:=simplify(RS):
end if:
#%

####>Einstein tensor
if sw>0 then 
  print("Einstein tensors: ET02[*,*], ET11[*,*]");
end if:
#%
ET02:=Matrix(1..dim,1..dim): # Einstein tensor
ET11:=Matrix(1..dim,1..dim): # ET11[i,j]=ET^i{}_j:
ii:='ii': for ii from 1 to dim do
  jj:='jj': for jj from ii to dim do
     ET02[ii,jj]:=Ric02[ii,jj]-1/2*RS*G[ii,jj]:
     if sw>1 then
       ET02[ii,jj]:=simplify(ET02[ii,jj]):
     end if:
     if jj=ii then 
        ET11[ii,jj]:=Ric11[ii,jj]-1/2*RS:
     else
        ET02[jj,ii]:=ET02[ii,jj];
        ET11[ii,jj]:=Ric11[ii,jj]:
     end if:
     if sw>1 then
        ET11[ii,jj]:=simplify(ET11[ii,jj]):
     end if:
     if ii<jj then
          ET11[jj,ii]:=ET11[ii,jj]:
     end if:
  end do:
end do:

####>Weyl tensor
if sw>0 then
  print("Weyl tensor: Weyl04[*,*,*,*]");
end if:
#%
Weyl04:=Array(1..dim,1..dim,1..dim,1..dim):
#%
if dim>2 then
ii:='ii': 
for ii from 1 to dim do
 jj:='jj':
 for jj from ii to dim do
  kk:='kk':
  for kk from 1 to dim do
   ll:='ll':
   for ll from kk to dim do
     Weyl04[ii,jj,kk,ll]:=Riem04[ii,jj,kk,ll]
            -1/(dim-2)*(G[ii,kk]*Ric02[jj,ll]-G[ii,ll]*Ric02[jj,kk]
                       +G[jj,ll]*Ric02[ii,kk]-G[jj,kk]*Ric02[ii,ll])
       +1/(dim-1)/(dim-2)*(G[ii,kk]*G[jj,ll]-G[ii,ll]*G[jj,kk])*RS:
     if sw>1 then
       Weyl04[ii,jj,kk,ll]:=simplify(Weyl04[ii,jj,kk,ll]):
     end if:
     Weyl04[ii,jj,ll,kk]:=-Weyl04[ii,jj,kk,ll]:
     Weyl04[jj,ii,kk,ll]:=-Weyl04[ii,jj,kk,ll]:
     Weyl04[jj,ii,ll,kk]:=Weyl04[ii,jj,kk,ll]:
   end do:
  end do:
 end do:
end do:
end if:

####>array to table
Riem13:=convert(Riem13,table):
Riem13["indextype"]:=[_u,_d,_d,_d]:
Riem13["dim"]:=dim:
Riem04:=convert(Riem04,table):
Riem04["indextype"]:=[_d,_d,_d,_d]:
Riem04["dim"]:=dim:
Riem22:=convert(Riem22,table):
Riem22["indextype"]:=[_u,_u,_d,_d]:
Riem22["dim"]:=dim:
if dim>2 then
  Weyl04:=convert(Weyl04,table):
  Weyl04["indextype"]:=[_d,_d,_d,_d]:
  Weyl04["dim"]:=dim:
end if:
Ric02:=convert(Ric02,table):
Ric02["indextype"]:=[_d,_d]:
Ric02["dim"]:=dim:
Ric11:=convert(Ric11,table):
Ric11["indextype"]:=[_u,_d]:
Ric11["dim"]:=dim:
ET02:=convert(ET02,table):
ET02["indextype"]:=[_d,_d]:
ET02["dim"]:=dim:
ET11:=convert(ET11,table):
ET11["indextype"]:=[_u,_d]:
ET11["dim"]:=dim:

####>output
Gdata:=table(["varlist"=varlist,"dim"=dim, "metric"=G,"invmetric"=GI,
  "Christoffel"=Ch, "Riem13"=Riem13, "Riem04"=Riem04, "Riem22"=Riem22,
  "Ric02"=Ric02, "Ric11"=Ric11, "RS"=RS, "ET02"=ET02, "ET11"=ET11]):
if dim>2 then
  Gdata["Weyl04"]:=Weyl04:
end if:
if sw>0 then
  print("Done.");
end if:
return copy(Gdata);
end proc:



#
######>NP coefficients
#####>proc:NPC
#% NPC(Gdata,NullTetrad[,sw]]) => NPCtbl
#% NullTetrad=[FB[1]=..., FB[4]=..]
#% NPCtbl:=table(["SC"=SCNP, "Weyl"=WeylNP, "Ric"=RicNP, "NullTetrad"=NullTerad, "dim"=dim, "varlist"=varlist, "metric"=G])
#% sw=0 => no output display, sw=1 => display the simplified Psi0,..,Psi4

NPC:=proc()

  local Gdata, NullTetrad, NB, NV, sw ,G, dim,xx,dxx,var, 
        f, rho, eta,
        OmegaM, OmegaM1, SCNP, Weyl04, WeylNB, WeylNP,
        Ric02,RS,RicNB, RicNP, NPCtbl:
  
  if _npassed=0 then
    error("Usage: NPC(Gdata,NullTetrad[, sw])");
  end if:
  Gdata:=_passed[1]:
  dim:=Gdata["dim"]:
  xx:=Gdata["varlist"]:
  if dim<>4 then 
    error("NP coefficients are defined only for four dimensions");
  end if:
  NullTetrad:=_passed[2]:
  var:='var':
  NB:=map(var->op(1,var),NullTetrad):
  sw:=0:
  if _npassed>=3 then
    sw:=_passed[3]:
  end if:

  
  ##% Spin coefficients: SCNP["*"]
  OmegaM:=ConnectionForm(NullTetrad, Gdata,1):
  eta:=DiagonalMatrix([<<0,-1>|<-1,0>>, <<0,1>|<1,0>>]):
  OmegaM1:=eta.OmegaM:
    
  SCNP["alpha"]:=1/2*coeff(-OmegaM1[2,1]+OmegaM1[4,3],NB[4]):
  SCNP["beta"]:=1/2*coeff(-OmegaM1[2,1]+OmegaM1[4,3],NB[3]):
  SCNP["rho"]:=-coeff(OmegaM1[3,1],NB[4]):
  SCNP["sigma"]:=-coeff(OmegaM1[3,1],NB[3]):
  SCNP["lambda"]:=coeff(OmegaM1[4,2],NB[4]):
  SCNP["mu"]:=coeff(OmegaM1[4,2],NB[3]):
  SCNP["kappa"]:=-coeff(OmegaM1[3,1],NB[3]):
  SCNP["pi"]:=coeff(OmegaM1[4,2],NB[1]):
  SCNP["epsilon"]:=1/2*coeff(-OmegaM1[2,1]+OmegaM1[4,3],NB[1]);
  SCNP["nu"]:=coeff(OmegaM1[4,2],NB[2]):
  SCNP["tau"]:= -coeff(OmegaM1[3,1],NB[2]):
  SCNP["gamma"]:=1/2*coeff(-OmegaM1[2,1]+OmegaM1[4,2],NB[2]):
print(SCNP);
 
  ##% Weyl: Psi[j]
  Weyl04:=Gdata["Weyl04"]:
  WeylNB:=Array(1..dim,1..dim,1..dim, 1..dim,
     CBtoFBT(Weyl04,xx,NullTetrad)):
  WeylNP["Psi[0]"]:=WeylNB[1,3,1,3]:
  if sw>0 then
    WeylNP["Psi[0]"]:=simplify(WeylNP["Psi[0]"]):
    print(Psi||"0"=simplify(WeylNP["Psi[0]"]));
  end if:
  WeylNP["Psi[1]"]:=WeylNB[1,2,1,3]:
  if sw>0 then
    WeylNP["Psi[1]"]:=simplify(WeylNP["Psi[1]"]):
    print(Psi||"1"=simplify(WeylNP["Psi[1]"]));
  end if:
 WeylNP["Psi[2]"]:= WeylNB[1,3,4,2]:
  if sw>0 then
    WeylNP["Psi[2]"]:=simplify(WeylNP["Psi[2]"]):
    print(Psi||"2"=simplify(WeylNP["Psi[2]"]));
  end if:
   WeylNP["Psi[3]"]:=WeylNB[2,1,2,4]:
  if sw>0 then
    WeylNP["Psi[3]"]:=simplify(WeylNP["Psi[3]"]):
    print(Psi||"3"=simplify(WeylNP["Psi[3]"]));
  end if:
  WeylNP["Psi[4]"]:=WeylNB[2,4,2,4]:
  if sw>0 then
    WeylNP["Psi[4]"]:=simplify(WeylNP["Psi[4]"]):
    print(Psi||"4"=simplify(WeylNP["Psi[4]"]));
  end if:
  if sw>0 then
      print(WeylNP);
  end if:
  
  ##% Ricci: Phi[i,j]
  Ric02:=Gdata["Ric02"]:
  RS:=Gdata["RS"]:
  RicNB:=Array(1..dim, 1..dim, CBtoFBT(Ric02,xx,NullTetrad)):

  RicNP["Phi[0,0]"]:=1/2*RicNB[1,1]:
  if sw>0 then
    RicNP["Phi[0,0]"]:=simplify(RicNP["Phi[0,0]"]):
    print(Phi||"00"=simplify(RicNP["Phi[0,0]"]));
  end if:
  RicNP["Phi[1,1]"]:=1/2*(RicNB[1,2]+RicNB[3,4]):
  if sw>0 then
    RicNP["Phi[1,1]"]:=simplify(RicNP["Phi[1,1]"]):
    print(Phi||"11"=simplify(RicNP["Phi[1,1]"]));
  end if:
  RicNP["Phi[2,2]"]:=1/2*RicNB[2,2]:
  if sw>0 then
    RicNP["Phi[2,2]"]:=simplify(RicNP["Phi[2,2]"]):
    print(Phi||"22"=simplify(RicNP["Phi[2,2]"]));
  end if:
  RicNP["Phi[0,1]"]:=1/2*RicNB[1,3]:
  if sw>0 then
    RicNP["Phi[0,1]"]:=simplify(RicNP["Phi[0,1]"]):
    print(Phi||"01"=simplify(RicNP["Phi[0,1]"]));
  end if:
  RicNP["Phi[1,0]"]:=1/2*RicNB[1,4]:
  if sw>0 then
    RicNP["Phi[1,0]"]:=simplify(RicNP["Phi[1,0]"]):
    print(Phi||"10"=simplify(RicNP["Phi[1,0]"]));
  end if:
  RicNP["Phi[0,2]"]:=1/2*RicNB[3,3]:
  if sw>0 then
    RicNP["Phi[0,2]"]:=simplify(RicNP["Phi[0,2]"]):
    print(Phi||"02"=simplify(RicNP["Phi[0,2]"]));
  end if:
  RicNP["Phi[2,0]"]:=1/2*RicNB[4,4]:
  if sw>0 then
    RicNP["Phi[2,0]"]:=simplify(RicNP["Phi[2,0]"]):
    print(Phi||"20"=simplify(RicNP["Phi[2,0]"]));
  end if:
  RicNP["Phi[1,2]"]:=1/2*RicNB[2,3]:
  if sw>0 then
    RicNP["Phi[1,2]"]:=simplify(RicNP["Phi[1,2]"]):
    print(Phi||"12"=simplify(RicNP["Phi[1,2]"]));
  end if:
  RicNP["Phi[2,1]"]:=1/2*RicNB[2,4]:
  if sw>0 then
    RicNP["Phi[2,1]"]:=simplify(RicNP["Phi[2,1]"]):
    print(Phi||"21"=simplify(RicNP["Phi[2,1]"]));
  end if:

  NPCtbl:=table(["SCNP"=SCNP, "WeylNP"=WeylNP, "RicNP"=RicNP]):
  
  return copy(NPCtbl);
end proc:


######>print procedures for Riemann
#####>proc:printCh
#% printCh(Gdata) => printing sum_{j1,j2} Ch[i,j1,j2] dx[j1] dx[j2]
printCh:=proc()
  local Ch, varlist,dim, dx,var, ii,jj1,jj2,tmp:
  
  if _npassed=0 then
    error("Usage: printCh(Gdata)");
  end if:
  Ch:=_passed[1]["Christoffel"]:
  dim:=_passed[1]["dim"]:
  varlist:=_passed[1]["varlist"]:
  var:='var':
  dx:=map(var -> d||var, varlist):
  
  ii:='ii':
  for ii from 1 to dim do
    jj1:='jj1': jj2:='jj2':
    tmp:=add(add(Ch[ii,jj1,jj2]*dx[jj1]*dx[jj2],jj1=1..dim),jj2=1..dim):
    jj1:='jj1':
    tmp:=collect(tmp,[seq(dx[jj1],jj1=1..dim)],distributed,simplify):
    print(Gamma||"^"||ii||"_** "=tmp);
#    printf("  Gamma^%a_(jk)dx^jdx^k = %a\n",ii, tmp);
  end do:
  
end proc:
    
#####>proc:printRiem
#% printRiem(Gdata[,Ttype]) => printing Riem04/13/22[*,*,*,*] <>0
#% Ttype = "04", "13", "22"
printRiem:=proc()
  local dim,Ttype, RiemT, ii1,ii2,jj1,jj2,tmp, outcount:
  
  if _npassed=0 then
    error("Usage: printRiem04(Gdata)");
  end if:
  dim:=_passed[1]["dim"]:
  Ttype:="13":
  if _npassed>1 then
    Ttype:=_passed[2]:
    if Ttype<>"04" and Ttype<>"13" and Ttype<>"22" then
      error("Tensor type value should be \"04\" or \"13\" or \"22\". ");
    end if: 
  end if:
  RiemT:=_passed[1][cat("Riem",Ttype)]:
 
  outcount:=0:
  ii1:='ii1': 
  for ii1 from 1 to dim do 
    ii2:='ii2': 
    for ii2 from ii1+1 to dim do 
      jj1:='jj1': 
      for jj1 from 1 to dim do 
        jj2:='jj2': 
        for jj2 from jj1+1 to dim do 
          if assigned(RiemT[ii1,ii2,jj1,jj2]) then
            tmp:=simplify(RiemT[ii1,ii2,jj1,jj2]):
            if tmp<>0 then 
              print(Riem||Ttype||"["||ii1||","||ii2||","||jj1||","||jj2||"]"
              =tmp);  outcount:=outcount+1:
            end if:
          end if: 
        end do: 
      end do: 
    end do: 
  end do:
  if outcount=0 then
    print("All components vanish.");
  end if:
end proc:

#####>proc:printWeyl
#% printWeyl(Gdata) => printing Weyl04[*,*,*,*] <>0
printWeyl:=proc()
  local dim,Weyl04, ii1,ii2,jj1,jj2,tmp, outcount:
  
  if _npassed=0 then
    error("Usage: printWeyl04(Gdata)");
  end if:
  dim:=_passed[1]["dim"]:
  Weyl04:=_passed[1]["Weyl04"]:
 
  outcount:=0:
  ii1:='ii1': 
  for ii1 from 1 to dim do 
    ii2:='ii2': 
    for ii2 from ii1+1 to dim do 
      jj1:='jj1': 
      for jj1 from 1 to dim do 
        jj2:='jj2': 
        for jj2 from jj1+1 to dim do 
          if assigned(Weyl04[ii1,ii2,jj1,jj2]) then
            tmp:=factor( Weyl04[ii1,ii2,jj1,jj2]):
            if tmp<>0 then 
              print(Weyl04||"["||ii1||","||ii2||","||jj1||","||jj2||"]" = tmp);  outcount:=outcount+1:
            end if:
          end if: 
        end do: 
      end do: 
    end do: 
  end do:
  
  if outcount=0 then
    print("All components vanish.");
  end if:

end proc:


#####>proc:printRic
#% printRic(Gdata [, Ttype]) => print Ric02/11[*,*]<>0
printRic:=proc()
  local Ttype,dim, RicT,ii1,ii2,tmp,outcount:
  
  if _npassed=0 then
    error("Usage: printRic02(Gdata)");
  end if:
  dim:=_passed[1]["dim"]:
  Ttype:="02":
  if _npassed>1 then
    Ttype:=_passed[2]:
    if Ttype<>"02" and Ttype<>"11" then
      error("Tensor type value should be \"02\" or \"11\". ");
    end if: 
  end if:
  RicT:=_passed[1][cat("Ric",Ttype)]:

  outcount:=0:
  ii1:='ii1': 
  for ii1 from 1 to dim do 
    ii2:='ii2': 
    for ii2 from ii1 to dim do 
       if assigned(RicT[ii1,ii2]) then
         tmp:=simplify( RicT[ii1,ii2]):
         if tmp<>0 then 
            print(Ric||Ttype||"["||ii1||","||ii2||"]" = tmp); 
            outcount:=outcount+1;
         end if:
       end if: 
    end do: 
  end do:
  
  if outcount=0 then
    print("All components vanish.");
  end if:

end proc:



#####>proc:printET
#% printET(Gdata[, Ttype]) => print ET02/11[*,*] <>0  
printET:=proc()
  local Ttype,dim,ETT,ii1,ii2,tmp,outcount:
  
  if _npassed=0 then
    error("Usage: printET02(Gdata)");
  end if:
  dim:=_passed[1]["dim"]:
  Ttype:="02":
  if _npassed>1 then
    Ttype:=_passed[2]:
    if Ttype<>"02" and Ttype<>"11" then
      error("Tensor type value should be \"02\" or \"11\". ");
    end if: 
  end if:
  ETT:=_passed[1][cat("ET",Ttype)]:

  outcount:=0:
  ii1:='ii1': 
  for ii1 from 1 to dim do 
    ii2:='ii2': 
    for ii2 from ii1 to dim do 
       if assigned( ETT[ii1,ii2]) then
         tmp:=simplify( ETT[ii1,ii2]):
         if tmp<>0 then 
            print(ET||Ttype||"["||ii1||","||ii2||"]" = tmp); 
            outcount:=outcount+1;
         end if:
       end if: 
    end do: 
  end do:
 
  if outcount=0 then
    print("All components vanish.");
  end if:

end proc:

#####>proc:printTensor
#% printTensor(T) => T[i,j,...]=...
printTensor:=proc()
  local T, dim,Tdim,indextype,
        indset,ind, nonzerocount,i;
  
  if _npassed=0 then
    error("Usage: printTensor(T::tensor table)");
  end if:
  T:=_passed[1]:
  if type(T,table) then
    dim:=T["dim"]:
    indextype:=T["indextype"]:
    Tdim:=nops(indextype):
  else 
    Tdim:=0:
  end if:
  
  if Tdim=0 then
    print("Tensor type is scalar.");
    print("value"=T);
  else
    print("Tensor type is tensor", 'indextype'=indextype);
    nonzerocount:=0:
    indset:=mklattice(Tdim,[1,dim]):
    ind:='ind': for ind in indset do
      if assigned(T[op(ind)]) and eval(T[op(ind)])<>0 then
         if nonzerocount=0 then
           print("Non-zero components are as follows:");
         end if:
         print(cat(T,convert(ind,string)) =eval( T[op(ind)]));
         nonzerocount:=nonzerocount+1:
       end if:
    end do:
    if nonzerocount=0 then
      print("All components vanish!");
    end if:
  end if: 
end proc:

######>Kretschmann invariants 
#print("Kretschmann invariants");
#%
#####>proc:Kretschmann
#% Kretschmann(Gdata) => the Kretschmann invariant
Kretschmann:=proc()
local Gdata,dim, Riem22,
      i,j,k,l, val:

if _npassed=0 then
  error("Usage: Kretschmann(Gdata)");
end if:
Gdata:=_passed[1]:
dim:=Gdata["dim"]:
Riem22:=Array(1..dim,1..dim,1..dim,1..dim,Gdata["Riem22"]):

 val:=0:
 i:='i': 
 for i from 1 to dim-1 do
  j:='j':
  for j from i+1 to dim do
    k:='k';
    for k from 1 to dim-1 do
       l:='l':
       for l from k+1 to dim do
          val:=val+4*Riem22[i,j,k,l]*Riem22[k,l,i,j]:
       end do:
    end do:
  end do:
 end do:
 return val;
end proc:

#####>proc:KretschmannS
#% KretschmannS(Gdata) => Kretschmann1, Kretschmann2 
#% for static metric with xx[1] in the static direction
KretschmannS:=proc()
local Gdata,dim, Riem22,
      ii, jj,
      Kretschmann1, Kretschmann2:

if _npassed=0 then
  error("Usage: KretschmannS(Gdata)");
end if:
Gdata:=_passed[1]:
dim:=Gdata["dim"]:
Riem22:=Array(1..dim,1..dim,1..dim,1..dim,Gdata["Riem22"]):


 Kretschmann1:=0:
 ii:='ii': 
 for ii from 2 to dim do
  Kretschmann1:=Kretschmann1+4*Riem22[1,ii,1,ii]^2:
  jj:='jj':
  for jj from ii+1 to dim do
    Kretschmann1:=Kretschmann1+8*Riem22[1,ii,1,jj]*Riem22[1,jj,1,ii]:
  end do:
 end do:
 Kretschmann2:=Kretschmann(Gdata)-Kretschmann1:

return [Kretschmann1, Kretschmann2];

end proc:

######>Tensor algebra
#####>proc:islinear
#% islinear (eq, basis) => true/false (eq=dform?)
islinear:=proc()
  local eq,basis:
  
  if _npassed<2 then
    error("Usage: islinear(eq, basis) ");
  end if:
  eq:=expand(_passed[1]):
  basis:=_passed[2]:
  if degree(eq,basis)=1 and ldegree(eq,basis)=1 then
     return true;
  else
     return false:
  end if:
end proc:

#####>proc:ContractT
#% ContractT(T::tensor table,[a,b],Gdata[,sw]) => CT::tensor table
#%
ContractT:=proc()

  local T, indpair, Gdata, sw, dim, indextype, Tdim, G, GI,
        TA, T0, s1,s2, Q, IDM,
        i,j, indset, ind, newind, CTtbl;
  
  if _npassed<3 then
    error("Usage: ContractT(T::tensor table, [a,b], Gdata[, sw]) ");
  end if:
  T:=_passed[1]:
  dim:=T["dim"]:
  indextype:=T["indextype"]:
  Tdim:=nops(indextype):
  if Tdim<2 then
     error("The rank of the tensor must be greater than 1");
  end if:
  indpair:=sort(_passed[2]):
  i:='i':
  if not type(indpair, list) or nops(indpair)<>2 or 
     not ({op(indpair)} subset {seq(i,i=1..dim)}) then
     error("Invalid index positions to contract");
  end if:
  Gdata:=_passed[3]:
  G:=Gdata["metric"]:
  GI:=Gdata["invmetric"]:
  sw:=0:
  if _npassed>=4 then
    sw:=_passed[4]:
  end if:
  
  indset:=mklattice(Tdim-2,[1,dim]):
  CTtbl:=table():
  IDM:=IdentityMatrix(dim):
  i:='i': TA:=Array(seq(1..dim,i=1..Tdim),T):
  s1:=indpair[1]: s2:=indpair[2]:
  
  newind:='newind': for newind in indset do
     i:='i': j:='j': 
     ind:=[op(newind[1..(s1-1)]),i,op(newind[s1..(s2-2)]),j,op(newind[(s2-1)..(Tdim-2)])]:
     Q:=add(add(eval(M0[i,j]*T0[op(ind)]),i=1..dim),j=1..dim):
     if indextype[s1]='_u' and indextype[s2]='_u' then
        CTtbl[op(newind)]:= eval(Q,{M0=G,T0=TA}):
      elif indextype[s1]='_d' and indextype[s2]='_d' then
        CTtbl[op(newind)]:= eval(Q,{M0=GI,T0=TA}):
      else 
        CTtbl[op(newind)]:= eval(Q,{M0=IDM,T0=TA}):
      end if:
   end do:
   
   if sw>0 then
      CTtbl:=simplify(CTtbl):
   end if:
   CTtbl["dim"]:=dim:
   CTtbl["indextype"]:=[op(indextype[1..(s1-1)]),op(indextype[(s1+1)..(s2-1)]),op(indextype[(s2+1)..Tdim])]:
   
   return copy(CTtbl);
   
end proc:
  
  
#####>proc:IP
#% Inner product of two tensors/vectors/forms
#% IP(V,W, Gdata) => val::scalar
#% V,W:: tensor tables/vectors/dforms
IP:=proc()
   local Gdata,dim,G,GI,V,W, xx, dxx,IDM, IDM0,FB, 
       Q1,Q2,Q, indset,ind1,ind2, indextypeV, indextypeW,
       G0,GI0, Tdim,T1,T2, FBlist,
       i,j,var,mu,nu:
   
   if _npassed<3 then 
     error("Usage: IP(V,W, Gdata)");
   end if:
   V:=simplify(_passed[1]):
   W:=simplify(_passed[2]):
   Gdata:=_passed[3]:
   G:=Gdata["metric"];
   GI:=Gdata["invmetric"]:
   dim:=Gdata["dim"]:
   IDM:=IdentityMatrix(dim):
   xx:=Gdata["varlist"];
   var:='var': dxx:=map(var->d||var,xx):
   
   if type(V,Vector) and type(W,Vector) then  # vectors
      V:=convert(V,list):
      W:=convert(W,list):
     if nops(V)<> dim or nops(W) <> dim then
        error( "Invalid vectors");
     end if:
     return add(add(G[i,j]*V[i]*W[j],i=1..dim),j=1..dim);
   elif type(V,Matrix) and type(W,Matrix) then # covectors
     if op(1,V)<>(1,dim) or op(1,W) <> (1,dim) then
        error( "Invalid covectors");
     end if:
     return add(add(GI[i,j]*V[1,i]*W[1,j],i=1..dim),j=1..dim);
   elif islinear(V,dxx) and islinear(W,dxx) then  # 1-forms
     V:=Matrix(map2(coeff,V,dxx)):
     W:=Matrix(map2(coeff,W,dxx)):
     return IP(V,W,Gdata);
   elif type(V,table) and type(W,table) then  # tensors
     indextypeV:=V["indextype"]:
     indextypeW:=W["indextype"]:
     Tdim:=nops(indextypeV):
     if nops(indextypeW)<>Tdim then
        error("Two tensors should have the same rank");
     end if:
     indset:=mklattice(Tdim,[1,dim]):
     mu:='mu': i:='i': Q1:=T1[seq(mu[i],i=1..Tdim)]:
     nu:='nu': i:='i': Q2:=T2[seq(nu[i],i=1..Tdim)]:
     Q:=Q1*Q2:
     i:='i': for i from 1 to Tdim do
       if indextypeV[i] ='_u' and indextypeW[i] ='_u' then
          Q:=G0[mu[i],nu[i]]*Q:
       elif indextypeV[i] ='_d' and indextypeW[i] ='_d' then
          Q:=GI0[mu[i],nu[i]]*Q:
       else
          Q:=IDM0[mu[i],nu[i]]*Q:
       end if:
     end do:
     i:='i': V:=Array(seq(1..dim,i=1..Tdim),V):
     i:='i': W:=Array(seq(1..dim,i=1..Tdim),W):
    ind1:='ind1': ind2:='ind2': 
     Q:=add(add(eval(Q,{mu=ind1, nu=ind2}), ind1 in indset),
            ind2 in indset):
     Q:= eval(Q,{T1=V,T2=W,G0=G,GI0=GI, IDM0=IDM}):
     return Q;
   else
     error("Invalid vectors/covectors/dforms/tensors");
   end if:
 
end proc:

######>General frame & Covariant derivative
#####>proc:Covder0
#% Covder0(T,i,[j,k,..],Gdata) => val::scalar
#% T is a list representing a tensor of type T["indextype"]
#Covder0:=proc(T::list,derind,indlist::list,Gdata::table)
Covder0:=proc()

  local T, derind, Tindex, indextype, Gdata, dim, Ch, xx,
        Tdim, ii, jj, tmpindex, val:
  
  if _npassed<4 then 
    error("Usage: Covder0(Tensor table, i, [j,k,..], Gdata)");
  end if:
  T:=_passed[1]: 
  Tdim:=0:
  if type(T,table) then 
     indextype:=T["indextype"]:#% [_u, _d, ..]
    if nops(indexype)=0 then
       T:=T[]:
    else 
       Tindex:=_passed[3]:
       Tdim:=nops(indextype):
    end if:
  else 
     Tindex:=():
  end if:
  derind:=_passed[2]:
  Gdata:=_passed[4]:
  dim:=Gdata["dim"]:
  if T["dim"]<>dim then
    error("The dimensions of the tensor and space do not match");
  end if:
  Ch:=Gdata["Christoffel"]:
  xx:=Gdata["varlist"]:
  
  #% T=scalar case
  if Tdim=0 then
     if indextype=[] then
       T:=T[]:
     end if:
     return diff(T,xx[derind]);
  end if:
  
  #% T: non-scalar tensor case
  ii:='ii': 
  T:=Array(seq(1..dim,ii=1..Tdim),T):
  
  if nops(Tindex)<>Tdim then 
    print("Tdim"=Tdim,"Tindex"=Tindex);
    error("Type of Tensor does not match the index list.");
  end if:
 
  val:=diff(T[op(Tindex)],xx[derind]):
  ii:='ii':
  for ii from 1 to Tdim do
      tmpindex:=Tindex:
      jj:='jj': for jj from 1 to dim do
        tmpindex[ii]:=jj:
        if indextype[ii]='_d' then
          val:=val-Ch[jj,derind,Tindex[ii]]*T[op(tmpindex)]:
        elif indextype[ii]='_u' then
          val:=val+Ch[Tindex[ii],derind, jj]*T[op(tmpindex)]:
        else
          error("Invalid indextype of T.");
        end if:
      end do:
    end do:
    return val;

end proc:

#####>proc:Covder
#% Covder(T,[i1,...,in],[j, k, ...],Gdata)=> D_i1 ... D_in T[j,k,,...]
#% T_{j,...}^{k...} is a tensor
#Covder:=proc(T,derind::list,Tindex::list,Gdata::table)

Covder:=proc()

  local T, Tarray,Tdim, derind, Tind, Gdata, indextype, dim, xx, Ch,
  dinddim,newderind,ii, jj, tmpindex, val:

  
  if _npassed<4 or (not type(_passed[2],list) )
      or (not type(_passed[3],list)) then 
    print("arg[3]"=_passed[3]):
    error("Usage: Covder(Ttbl, [i,..], [j,k,..], Gdata)");
  end if:
  T:=copy(_passed[1]): 
  
  Tdim:=0:
  if type(T,table) then 
     indextype:=T["indextype"]:
    if nops(indexype)=0 then
       Tarray:=T[]:
    else 
       Tind:=_passed[3]:
       Tdim:=nops(indextype):
    end if:
  else 
     Tind:=[]:
     Tarray:=T:
  end if:
  if Tdim>0 and nops(indextype)<>Tdim then 
    error("Type of Tensor does not match the index list.");
  end if:
  derind:=_passed[2]:
  Gdata:=_passed[4]:
  dim:=Gdata["dim"]:
  if T["dim"]<>dim then
    error("The dimensions of the tensor and space do not match");
  end if:
  Ch:=Gdata["Christoffel"]:
  xx:=Gdata["varlist"]:
  if Tdim>0 then
    ii:='ii':  Tarray:=Array(seq(1..dim,ii=1..Tdim),T):
  end if:

  dinddim:=nops(derind):
  if dinddim=0 then
     if Tdim=0 then
        return Tarray;
     else
        return Tarray[op(Tind)];
     end if:
  end if:
  
  newderind:=derind[2..dinddim]:
  val:=diff(Covder(T,newderind,Tind,Gdata),xx[derind[1]]):
  ii:='ii': for ii from 1 to dinddim-1 do
    jj:='jj': for jj from 1 to dim do
        tmpindex:=newderind:
        tmpindex[ii]:=jj:
        val:=val-Ch[jj,derind[1],newderind[ii]]*Covder(T,tmpindex,Tind,Gdata):
    end do:
  end do:
  ii:='ii': for ii from 1 to Tdim do
    jj:='jj': for jj from 1 to dim do
      tmpindex:=Tind:
      tmpindex[ii]:=jj:
      if indextype[ii]='_d' then 
        val:=val-Ch[jj,derind[1],Tind[ii]]*Covder(T,newderind,tmpindex,Gdata):
      else
        val:=val+Ch[Tind[ii],derind[1],jj]*Covder(T,newderind,tmpindex,Gdata):
      end if:
    end do
  end do:
   
  return val;
  
end proc:


#####>proc:Laplacian
#% Laplacian of tensors
#% T ; tensor table => \triangle T: tensor table
#% Laplacian(T,indlist,Gdata)
Laplace:=proc()

  local T, indlist, Gdata,GI, dim,
        i1,i2:

  if _npassed<3 then
    error("Usage: Laplace(T::tensor table, indexlist, Gdata)");
  end if:
  T:=_passed[1]:
  indlist:=_passed[2]:
  Gdata:=_passed[3]:
  dim:=Gdata["dim"]:
  if T["dim"]<>dim then
    error("The dimensions of the tensor and space do not match");
  end if:

  GI:=Gdata["invmetric"]:

  i1:='i1': i2:='i2':
  add(add(GI[i1,i2]*Covder(T,[i1,i2],indlist,Gdata),i1=1..dim),i2=1..dim);

end proc:

 

#####>CovDnT
####>proc:mklattice
#% mklattice(rn::integer,[n1,n2]) => lattice point list
mklattice:=proc()
local rn,rng,
      sublattice,newlattice,x,i:

rn:=_passed[1]:
if _npassed<2 or not type(rn,integer) or rn<1 then
   error("usage: mklattice(rn::positive integer, [n1::integer,n2::integer])");
end if:
rng:=_passed[2]:
if not type(rng,list) or nops(rng)<2 or not type(rng[1],integer) 
   or not type(rng[2],integer) then
   error("usage: mklattice(rn::positive integer, [n1::integer,n2::integer])");
end if:

if rn=1 then
  i:='i': newlattice:=[seq([i],i=rng[1]..rng[2])]:
else
  sublattice:=mklattice(rn-1,rng):
  newlattice:=[]:
  for x in sublattice do
    i:='i':newlattice:=[op(newlattice),seq([op(x),i],i=rng[1]..rng[2])]:
  end do:
end if:
return newlattice;
end proc:


####>proc:CovDT
#% CovDT(T::table,Gdata::table[,sw]) => DT::table
CovDT:=proc()
  local T, Gdata,sw, dim,xx, Tdim, indextype, Tval,DT,
        i, p, ind, indexset;
  
  if _npassed<2 then
     error("Usage: CovDT(T, Gdata)");
  end if:
  T:=_passed[1]:
  Gdata:=_passed[2]:
  dim:=Gdata["dim"]:
  xx:=Gdata["varlist"]:
  if not type(T,table) then
    Tdim:=0: Tval:=T: #% regarded as a scalar
  else
    indextype:=T["indextype"]:
    Tdim:=nops(indextype):
    if Tdim=0 then
       Tval:=T[]:
    end if:
  end if:
  sw:=0:
  if _npassed>=3 then
    sw:=_passed[3]:
  end if:
  if Tdim>0 and T["dim"]<>dim then
    error("The dimensions of the tensor and space do not match");
  end if:


  if Tdim=0 then  #% scalar case
    DT:=Array(1..dim):
    i:='i': for i from 1 to dim do 
      DT[i]:=diff(Tval,xx[i]): 
      if sw>0 then
        DT[i]:=simplify(DT[i]):
      end if:
    end do:
    DT:=convert(DT,table):
    DT["indextype"]:=['_d']:
  else #% tensor case
    i:='i': DT:=Array(seq(1..dim,i=1..(Tdim+1))):
    indexset:=mklattice(Tdim,[1,dim]):
    for ind in indexset do 
        i:='i': for i from 1 to dim do
          DT[i,op(ind)]:=Covder0(T,i,ind,Gdata):
          if sw>0 then
            DT[i,op(ind)]:=simplify(DT[i,op(ind)]):
           end if:
         end do:
    end do:
    DT:=convert(DT,table):
    DT["indextype"]:=[_d,op(indextype)]:
    DT["dim"]:=dim:
  end if:
  return copy(DT):

end proc:

####>proc:CovDnT
#% CovDnT(T::tensor table, n::power, Gdata[,sw]) => D^nT
CovDnT:=proc()

  local T, n, Gdata, sw:
  
  if _npassed<3 then 
    error("Usage: CovDnT(T, n, Gdata[,sw])");
  end if:
  T:=_passed[1]:
  n:=_passed[2]:
  Gdata:=_passed[3]:
  sw:=0:
  if _npassed>=4 then
    sw:=_passed[4]:
  end if:
  
  if n=0 then
    return T;
  elif n=1 then
    return CovDT(T,Gdata,sw);
  else
    return CovDT(CovDnT(T,n-1,Gdata,sw),Gdata);
  end if:
end proc:

 
#####>proc:DV
#% DV(X,Y,Gdata): D_X Y ; X,Y vectors by lists or arrays => list
DV:=proc()

   local X,Y,Gdata,dim, Z,ii,jj:
   
   if _npassed<3 then
     error("Usage: DV(X, Y, Gdata)");
   end if:
   Gdata:=_passed[3]:
   dim:=Gdata["dim"]:
   X:=convert(_passed[1],list):
   Y:=convert(_passed[2],list):
   
   if nops(X)<>dim or nops(Y)<>dim then
     error("Invalid dimension of the vectors");
   end if:
   
   Y:=convert(Y,table): Y["indextype"]:=[_u]: Y["dim"]:=dim:
   Z:=Array(1..dim):
   ii:='ii':
   for ii from 1 to dim do
      jj:='jj':
      Z[ii]:=add(X[jj]*Covder0(Y,jj,[ii],Gdata),jj=1..dim):
   end do:
   ii:='ii':
   return [seq(Z[ii],ii=1..dim)];
end proc:


#####>proc:Div
#% Div(T,indexlist, Gdata) =>  D_i T^i{}_j...
Div:=proc()
  local T, indextype,indexlist, Gdata, GI, dim,
  i1, i2, val:

  if _npassed<3 then
    error("Usage: Div(T, indexlist, Gdata[, sw])");
  end if:
  T:=_passed[1]:
  indextype:=T["indextype"]:
  indexlist:=_passed[2]:
  Gdata:=_passed[3]:
  dim:=Gdata["dim"]:
  if T["dim"]<>dim then
    error("The dimensions of the tensor T and space do not matsch.");
  end if:
  GI:=Gdtaa["invmetric"]:
  
  if indextype[1]='_u' then
    i1:='i1': val:=add(Covder(T,[i1],[i1,op(indexlist)]),i1=1..dim):
  else
    i1:='i1': i2:='i2':
    val:=add(add(GI[i1,i2]*Covder(T,[i1],[i2,op(indexlist)]),i1=1..dim),i2=1..dim):
  end if:
  
  return val;

end proc:

#####>proc:DivT
#% DivT(T, Gdata[,sw]) => DivTtbl::tensor table: DivTtbl[a,b,...]=add(D_p T[p,a,b,..],a=1..dim) 
DivT:=proc()
  local T,  Gdata, indextype,sw, GI, dim,Tdim,
  i1, i2, DT, CDT, val:

  if _npassed<2 then
    error("Usage: DivT(T, Gdata[, sw])");
  end if:
  T:=_passed[1]:
  indextype:=T["indextype"]:
  Tdim:=nops(indextype):
  if Tdim<1 then
    error("The tensor should have a positive rank");
  end if:
  Gdata:=_passed[2]:
  dim:=Gdata["dim"]:
  if T["dim"]<>dim then
    error("The dimensions of the tensor T and space do not matsch.");
  end if:
  GI:=Gdtaa["invmetric"]:
  sw:=0:
  if _npassed>=3 then
    sw:=_passed[3]:
  end if:
  
  CDT:=ContractT(CovDT(T, Gdata),[1,2],Gdata):
  CDT["dim"]:=dim:
  CDT["indextype"]:=indextype[2..Tdim]:
  if sw>0 then
    CDT:=simplify(CDT):
  end if:
  return copy(CDT);

end proc:

#####>proc:CordBV
#% CordBV(i,dim)=> [0,..0,1,0,...0] (i-th component)
#% coordinate basis vector 
CordBV:=proc()

   local i,dim,j,k:
   
   if _npassed<2 then
     error("Usage: CordBV(i::positive integer, dim::positive integer)");
   end if:
   
   i:=_passed[1]: dim:=_passed[2]:
   if not (type(i,integer) and i>0 and type(dim,integer) and dim>0 
       and i<=dim ) then
      error("Usage: CordBV(i::positive integer, dim::positive integer), i<=dim) ");
   end if:
   
   return [seq(0,j=1..(i-1)),1,seq(0,k=(i+1)..dim)]:
   
end proc:

######>Connection and curvature forms
#####>proc: ConnectionForm
#% ConnectionForm(FormBasis::list, Gdata[, sw]) => OmegaM=(omega[a,b])
#% FormBasis:[FB[1]=theta[1],...,FB[n]=theta[n]] => Spin coefficients: omega[a,b]
#% theta[i]= sum_m theta[i,m] dxx[m]
#% sw=0 => omega[a,b]=omega[a,b][mu] dx^\mu
#% sw=1 => omega[a,b]=omega[a,b][c] FB[c]

ConnectionForm:=proc()

   local   FormBasis,Gdata, sw, dim,xx, dxx,
           var, EM, ThetaM, FB1, 
           i,j, mu,nu, dxbyFB, DEM,OmegaM:
   
   if _npassed<2 then 
     error("Usage: SpinC(FormBasis::list, Gdata[, sw])");
   end if:
   FormBasis:=convert(_passed[1],list):   
   i:='i': FB1:=[seq(op(1,FormBasis[i]),i=1..nops(FormBasis))]:
   Gdata:=_passed[2]:
   dim:=Gdata["dim"]:
   xx:=Gdata["varlist"]:
   var:='var':
   dxx:=map(var -> d||var, xx):
   if nops(FormBasis)<>dim then
     print("nops(FormBasis)"=nops(FormBasis),"dim"=dim);
     error("Dimension of the form basis does not match Gdata") ;
   end if:
   sw:=0:
   if _npassed>=3 then
     sw:=_passed[3]:
   end if:
   
   ThetaM:=Matrix(dim,dim): #=(\theta^a{}_\mu)
   i:='i': for i from 1 to dim do
     mu:='mu': for mu from 1 to dim do
        ThetaM[i,mu]:=coeff(op(2,FormBasis[i]),dxx[mu]):
     end do:
   end do:
   EM:=Transpose(ThetaM^(-1)):  #=(e_a{}^\mu)
   mu:='mu':i:='i': 
   DEM:= < seq(
     add(dxx[mu]*Matrix(DV(CordBV(mu,dim),EM[i],Gdata)),mu=1..dim),
     i=1..dim)>: # DEM[i,nu]=e_j^nu omega^j{}_i
   OmegaM:=ThetaM.Transpose(DEM):
   OmegaM:=simplify(OmegaM):
   
   if sw>0 then
     i:='i': mu:='mu': 
     dxbyFB:=[seq(dxx[mu]=add(FB1[i]*EM[i,mu],i=1..dim), mu=1..dim)]:
     OmegaM:=map(collect,eval(OmegaM,dxbyFB),FB1,simplify):      
   end if:
   
   return OmegaM:
       
end proc:
 
#####>proc:RiemForm
#% RiemForm(FormBasis,Gdata[, sw]) => Matrix( RF[a,b])
#% Calculating the Curvature form
#% FormBasis:[FB[1]=theta[1],...,FB[n]=theta[n]] 
#% theta[i]= sum_m theta[i,m] dxx[m]
#% sw=0 => RF[a,b]=(1/2)R^a_{b mu nu} dxdx[mu,nu]
#% sw=1 => RF[a,b]=(1/2)R^a_{b c d} FB2[c,d]

RiemForm:=proc()

   local   FormBasis, Gdata,sw,FB1, dim, xx, var, ThetaM, EM,  
           dxx, dx2B, Riem13,RF, FB2list,dx2BbyFB,
     i,j, k,l,a,b,c,d,mu:
   
   #%% set data
   if _npassed<2 then
     error("Usage: RiemForm(FormBasis::list, Gdata[, sw])");
   end if:
   
   FormBasis:=convert(_passed[1],list):
   i:='i': FB1:=[seq(op(1,FormBasis[i]),i=1..nops(FormBasis))]:
   Gdata:=_passed[2]:
   dim:=Gdata["dim"]:
   xx:=Gdata["varlist"]:
   Riem13:=Array(1..dim,1..dim,1..dim,1..dim,Gdata["Riem13"]):
   dxx:=table(antisymmetric):
   i:='i': for i from 1 to dim do
     dxx[i]:=d||(xx[i]):
   end do:
   i:='i': for i from 1 to dim do
     j:='j': for j from i+1 to dim do
       dxx[i,j]:=FB[dxx[i],dxx[j]]:
     end do:
   end do:
   
   if nops(FormBasis)<>dim then
     error("Dimension of the form basis does not match Gdata") ;
   end if:
   sw:=0:
   if _npassed>=3 then
     sw:=_passed[3]:
   end if:
   
   #%% define basis matrices
   ThetaM:=Matrix(dim,dim): #=(\theta^a{}_\mu)
   i:='i': for i from 1 to dim do
     mu:='mu': for mu from 1 to dim do
        ThetaM[i,mu]:=coeff(op(2,FormBasis[i]),dxx[mu]):
     end do:
   end do:
   EM:=Transpose(ThetaM^(-1)):  #=(e_a{}^\mu)

   RF:=Matrix(dim,dim):
   a:='a': for a from 1 to dim do
     b:='b': for b from 1 to dim do
        i:='i': j:='j': k:='k': l:='l':
        RF[a,b]:=1/2*add(add(add(add(Riem13[i,j,k,l]*ThetaM[a,i]*EM[b,j]
                *dxx[k,l],i=1..dim),j=1..dim),k=1..dim),l=1..dim):
     end do:
   end do:
 
   if sw=0 then
     FB2list:=[]:
     i:='i': for i from 1 to dim-1 do
       j:='j': for j from i+1 to dim do
         if op(1,dxx[i,j])=-1 then
            FB2list:=[op(FB2list),-dxx[i,j]]:
         else
            FB2list:=[op(FB2list),dxx[i,j]]:
         end if:
       end do:
     end do:
   else 
     dx2BbyFB:=[]:
     i:='i': for i from 1 to dim-1 do
       j:='j': for j from i+1 to dim do
         a:='a': b:='b':
         if op(1,dxx[i,j])=-1 then
             dx2BbyFB:=[op(dx2BbyFB),
                      -dxx[i,j]=-add(add(EM[a,i]*EM[b,j]*FB[a,b],
                        a=1..dim),b=1..dim)]:
         else
             dx2BbyFB:=[op(dx2BbyFB),
                      dxx[i,j]=add(add(EM[a,i]*EM[b,j]*FB[a,b],
                        a=1..dim),b=1..dim)]:
         end if:
       end do:
     end do:
     RF:=eval(RF,dx2BbyFB):
     i:='i': j:='j':
     FB2list:=[seq(seq(FB[i,j],j=(i+1)..dim),i=1..(dim-1))]:
   end if:
   RF:=map(collect,RF,FB2list,simplify):
   return RF:
   
end proc:
    
######>Basis change
#####>proc:CBtoFBT
#% CBtoFBT(T,varlist,FormBasis[,sw])
#% transforming the component rep of a tensor T from the coord basis to a general basis
#% T::tensor table with indextype info
#% varlist= the list of coordinate names
#% FormBasis=[FB[1]=theta[1], ...]
#% sw=0 => no simplification
#% sw=1 => simplification

CBtoFBT:=proc()

  local T, FormBasis, sw,indextype,dim,xx,dxx,Tdim, 
         ThetaM, EM, ThetaMtbl, EMtbl,Ttbl,
        FB1,T1,indset,Q, Tnew, ind, sind, var,i,j,mu,nu ;
  
  #%% set dat  
  if _npassed<3 then
    error("Usage: CBtoFBT(T::tensor table, xx::varlist, FormBasis[,sw])");
  end if:
  T:=copy(_passed[1]):
  indextype:=T["indextype"]:
  if not type(indextype,list) then
    error("The tensor table T should increase the indextype info");
  end if:
  Tdim:=nops(indextype):
  if Tdim=0 then
    return T;
  end if:
  Ttbl:=table():
  
  xx:=_passed[2]:
  dim:=nops(xx):
  if T["dim"]<>dim then
    error("The dimensions of the tensor and space do not match");
  end if:
  var:='var':
  dxx:=map(var->d||var,xx):
  FormBasis:=convert(_passed[3],list):
  i:='i': FB1:=[seq(op(1,FormBasis[i]),i=1..nops(FormBasis))]:
  sw:=0:
  if _npassed>=4 then
    sw:=_passed[4]:
  end if:
  i:='i': T:=Array(seq(1..dim,i=1..Tdim),T):
 
  
  #%% define basis matrices
   ThetaM:=Matrix(dim,dim): #=(\theta^a{}_\mu)
   i:='i': for i from 1 to dim do
     mu:='mu': for mu from 1 to dim do
        ThetaM[i,mu]:=coeff(op(2,FormBasis[i]),dxx[mu]):
     end do:
   end do:
   EM:=Transpose(ThetaM^(-1)):  #=(e_a{}^\mu)
   ThetaMtbl:=convert(ThetaM,table):
   EMtbl:=convert(EM,table):

  
  #%% trf to the general FB
  indset:=mklattice(Tdim,[1,dim]):
  mu:='mu': i:='i': Q:=Ttbl[seq(mu[i],i=1..Tdim)]:
  i:='i': for i from 1 to Tdim do
    if indextype[i] ='_u' then
       Q:=ThetaMtbl[a[i],mu[i]]*Q:
    else
       Q:=EMtbl[a[i],mu[i]]*Q:
    end if:
  end do:
  i:='i': Tnew:=Array(seq(1..dim,i=1..Tdim)):
  ind:='ind': for ind in indset do
    sind:='sind':
    Tnew[op(ind)] :=add(eval(Q,{a=ind, mu=sind}),sind in indset):
    Tnew[op(ind)] := eval(Tnew[op(ind)],{Ttbl=T,ThetaMtbl=ThetaM,EMtbl=EM}):
    if sw>0 then
      Tnew[op(ind)]:=simplify(Tnew[op(ind)]):
    end if:
  end do:
  Tnew:=convert(Tnew,table);
  Tnew["indextype"]:=indextype:
  Tnew["dim"]:=dim:
  return copy(Tnew):
end proc:

######>help
#####>proc:procinfo_Riemann
procinfo_Riemann:=proc()

  
  ##% General info
  printf("* General information:\n");
  printf("  When you use this package, you have to prepare a coordinate name list, say varlist=%a, \n", [t,r,theta,phi]);
  printf("  and a metric ds2 expressed as a quadratic equation in the differential of the coordinate names,\n");
  printf("  say dvarlist=%a and \n",[dt,dr,dtheta,dphi]);
  printf("  ds2 = %a\n", -(1-2*M/r)*dt^2+dr^2/(1-2*M/r)+r^2*(dtheta^2+dphi^2*sin(theta)^2));
  printf("\n");
  
  printf("* Global variables\n");
  printf("  _u, _d, FB\n");
  printf("\n");
  
  printf("* Variable name convention in this help message\n");
  printf("  + <dim> means the dimension of the space(time) under consideration.\n");
  printf("  + <T::tensor table> means that the tensor T is expressed by a table with \n");
  printf("        the entries specifying the value of components as T[1,2]=a => (1,2)=a, \n"); 
  printf("        the entry \"indextype\"=[_u,_d,...] specifying the index positions, and \n");
  printf("        the entry specifying the space dimension, say \"dim\"=4.  \n");
  printf("     When you read out the component values , you must first convert T to Array TA by \n");
  printf("        TA=Array(1..dim, ..., 1..dim, T) \n");
  printf("     Then, you can get the component value as TA[1,2]. This is because only the information on the components \n");
  printf("      with non-vanishing values are left in the table obtained by type conversion from an Array. \n");  
  printf("   + <varlist> means the coordinate name list.\n");
  printf("   + <dvarlist> means the list of the differential of the coordinate names\n");
  printf("   + <ds2> means the metric expressed as an equation quadratic in dvarlist\n");
  
  printf("   + <Gdata> means a table of the geometric data\n");
  printf("         (varlist, dim, metric, invmetric, Christoffel, Riem13, Riem04, Riem22, Weyl04, Ric02, Ric11, RS, ET02, ET11)\n");
  printf("      Here, \"metric\" and \"invmetric\" are the matrix representing the space metric and its inverse matrix, respectively.\n");
  printf("      \"Christoffel\" refers to Ch:=Array(1..dim,1..dim, 1..dim) repsenting the Christoffel symbol (Levi-Civita connection coefficients).\n");
  printf("      \"Riem13\", \"Riem04\" and \"Riem22\" represent the tensor tables corresponding to the (1,3)-type, (0,4)-type, \n");
  printf("      and (2,2)-type expressions of the Riemann curvature tensors, respectively, and \"Weyl04\" means \n");
  printf("      the tensor table for the Weyl tensor of the type (0,4). \n");
  printf("      \"Ric02\" and \"Ric11\" represent the tensor tables of the Ricci tensors with index type (0,2) and (1,1), respectively.\n");
  printf("      \"ET02\" and \"ET11\" are the tensor tables of the Einstein tensors, and \"RS\" is the scalar curvature.\n");
  printf("      Each entry object can be referred to as, say, Gdata[\"dim\"]. \n");
  printf("  \n");
  
  ##* The central procedure
  printf("* Central procedures\n");
  #% Riemann
  printf("  + Riemann(varlist,ds2[,sw]) => Gdata \n");
  printf("    procedure to calculate the connection coefficients, Riemann curvature tensors, Weyl tensor, \n");
  printf("    Ricci tensor, Ricci scalar and Einstein tensor for a given metric. \n");
  printf("    sw=0 => no simplification of results and no message (default)\n");
  printf("    sw=1 => no simplification with diplay of progress messages\n")+
  printf("    sw=2 => simplify results with progress messages \n");
  printf("    For example,\n");
  printf("       Gdata:=Riemann(varlist,ds2);  Riem13tbl;=Gdata[\"Riem13\"]\n");
  printf("    gives the tensor table Riem13tlb for the (1,3)-type Riemann curvature tensor of the metrid ds2, and by the conversion\n");
  printf("       Riem13:=convert(1..dim,1..dim,1..dim,1..dim, Riem13tbl)\n");
  printf("    you can get the (1,2,1,2) component of the tensor as Riem13[1,2,1,2]\n");
  
  ##% Printing procedures
  printf("* Printing procedures\n");
  printf("  The following are procedures to show the data in Gdata on the monitor. \n");
  #% printCh
  printf("  + printCh(Gdata) =>add(add(Ch[a,b,c]*dx[b]*dx[c],b=1..dim),c=1..dim)\n");  
  #% printRiem
  printf("  + printRiem(Gdata, Ttype) => Riem04/Riem13/Riem22[*,*,*,*]<>0 \n");
  printf("    Ttype = \"13\", \"04\", \"22\" \n");
  #% printWeyl
  printf("  + printWeyl(Gdata) => Weyl04[*,*,*,*] <>0 \n");
  #% printRic
  printf("  + printRic(Gdata,Ttype) => Ric02/Ric11[*,*] <>0 \n");
  printf("    Ttype = \"02\", \"11\" \n");
  #% printET
  printf("  + printET(Gdata,Ttype) => ET02/ET11[*,*] <>0 \n");
  printf("    Ttype = \"02\", \"11\" \n");
  #% printTensor
  printf("  + printTensor(T::tensor table) => T[*,*,..]<>0 \n");
  printf("    This is a procedure to display all non-vanishing components of a given tensor table T of any type. \n");
  
  ##% Special procedures
  printf("* Special procedures\n");
  #% Kretschmann/KreschmannS
  printf("  + Kretschmann( Gdata ) => the Kretschmann invariant K \n");
  printf("      K=add(add(add(add(Riem22[a,b,c,d]*Riem22[c,d,a,b], a=1..dim),b=1..dim),c=1..dim),d=1..dim) \n");
  printf("  + KretschmannS( Gdata ) => (K1,K2) [only for a static metric]\n");
  printf("      K1=4*add(add(Riem22[1,i,1,j]*Riem22[1,j,1,i],i=2..dim),j=2..dim)\n ");
  printf("      K2=add(add(add(add(Riem22[i,j,k,l]*Riem22[k,l,i,j], i=2..dim),j=2..dim),k=2..dim),l=2..dim) \n "):
  #% Newman-Penrose
  printf("  + NPC(Gdata, NullTetrad[, sw]) => NPCtbl [dim=4 only!!]\n");
  printf("      procedure to calculate the spin coefficients and the Newman-Penrose coefficients \n");
  printf("      for the Weyl curvature and the Ricci curvature w.r.t. a given null frame in four space dimensions. \n");
  printf("          NullTetrad::list = [NB[1]=k_*, NB[2]=l_*, NB[3]=m_* NB[4]=conjugate(m_*)] \n");
  printf("      where k, l, m, conjugate(m) are a linearly independent null 1-forms reprensenting a null tetrad \n");
  printf("      normalized as g(k,k)=g(l,l)=g(m,m)=0, g(k.l)=-1, g(m,conjugate(m))=1. \n");
  printf("      NPCtbl= table([\"SCNP\"=SCNP, \"WeylNP\"=WeylNP, \"RicNP\"=RicNP, \"NullTetrad\"=NullTerad, \"dim\"=dim, \"varlist\"=varlist, \"metric\"=G])\n");
  printf("      where SCNP is a table of the NP spin coefficients which can be referred to as e.g. NPCtbl[\"SCNP\"][\"alpha\"]. \n");
  printf("      WeylNP is the NP coefficients Psi[0], .., Psi[4] for the Weyl tensor, referred to as NPCtbl[\"WeylNP\"][\"Psi[0]\"].\n");
  printf("      RicNP is the NP coefficients Phi[0,0],...,Phi[2,2] for the Ricci tensor, referred to as NPCtbl[\"RicNP\"][\"Phi[0,0]\"] \n");
  printf("         sw=0 => no output display \n");
  printf("         sw=1 => display the simplified Psi[*] and Phi[*,*] \n");
  
  ##% General frame
  printf("* General frame \n");
  #% Frame change
  printf("  + CBtoFBT(T::tensor table,varlist,FormBasis[,sw]) => T1::tensor table w.r.t. the frame FormBasis \n");
  printf("    transforming the tensor T from the coord basis to a general basis.\n");
  printf("    T = a tensor table with components in the coordinate basis\n");
  printf("    FormBasis::list = [FB[1]=theta[1], ..., FB[dim]=theta[m]] \n");
  printf("    where theta[i]'s are  linearly independent 1-forms defining a covector basis.\n"); 
  printf("        sw=0 => no simplification \n");
  printf("        sw=1 => simplification of the results \n");
  
  ##% Connection form and curvature form
  printf("* Connection and curvature forms\n");
  printf("  + ConnectionForm(FormBasis::list, Gdata[, sw]) => CFM=(omega[a,b]) \n");
  printf("    CFM is the matrix whose component omega[a,b] is the connection form w.r.t. the FormBasis\n");
  printf("      FormBasis = [FB[1]=theta[1],..,FB[dim]=theta[dim]] \n");
  printf("      theta[a]= add(theta[a,mu]*dx[mu], mu=1..dim) \n");
  printf("           sw=0 => omega[a,b]=add(omega[a,b][mu]*dx[mu],mu=1..dim) \n");
  printf("           sw=1 => omega[a,b]=add(omega[a,b][c]*FB[c],c=1..dim) \n");
  printf("  + RiemForm(FormBasis::list,Gdata[, sw]) => RFM::Matrix=( RF[a,b]) \n");
  printf("    Calculating the curvature form RFa,b] w.r.t. the FormBasis\n");
  printf("      sw=0 => RF[a,b]=(1/2)add(add(Riem13[a,b mu nu]*FB[dx[mu],dx[nu]],mu=1..dim),nu=1..dim) \n");
  printf("      sw=1 => RF[a,b]=(1/2)add(add(Riem13[a,b c d]*FB[c,d],c=1..dim),d=1..dim) \n");
  printf("    where mu and nu refer to coordinate labels and a,b,c,d refer to frame labels of the FormBasis. \n");  
  
  ##% Covariant derivative of tensors
  printf("* Covariant derivative of tensors \n");
  #% Covder
  printf("  + Covder(T,[a[1],...,a[n]],[b, c, ...],Gdata)=> (D_a[1] ... D_a[n] T)[b,c,,...] \n");
  printf("    calulating each components of the n-th covariant derivative of a tensor T.\n");
  printf("    T = a tensor table or a scalar \n");
  printf("    No simplification is applied. \n");
  #% CovDT
  printf("  + CovDT(T:: tensor table, Gdata[, sw]) => DT \n");
  printf("    constructing a tensor table DT correspoinding to the covariant derivative of the given tensor table T \n");
  printf("        sw=0 => no simplification \n");
  printf("        sw=1 => simplification of the results \n");
  #% CovnDT
  printf("  + CovnDT(T::tensor table, n::positive integer, Gdata[, sw]) => DnT \n");
  printf("    constructing a tensor table DnT correspoinding to the n-the covariant derivative of the given tensor table T \n");
  printf("        sw=0 => no simplification \n");
  printf("        sw=1 => simplification of the results \n");
  
  #% Laplacian
  printf("  + Laplacian(T::tensor table, [a[1],...,a[n]], Gdata) => (Laplacian T)[a[1],..a[n]] \n");
  #% DV
  printf("  + DV(X,Y,Gdata) =>  D_X Y ::list \n");
  printf("    calculates the covariant derivative of two vectors X and Y.\n");
  printf("    X,Y = vectors expressed as lists, Vectors, Matrices or Arrays.\n");
  printf("    Output D_X Y is given in the form of list. \n");
  #% Div
  printf("  + Div(T::tensor table,[mu,nu,..], Gdata) => add((DT)[alpha,alpha,mu, nu,..],alpha=1..dim) \n");
  printf("  + DivT(T::tensor table, Gdata[, sw])=> DivTtbl::tensor table:\n");
  printf("    DivTtbl[mu,nu,...]=Div(T, [mu,nu,...], Gdata) \n");
  printf("        sw=0 => no simplification \n");
  printf("        sw=1 => simplification of the results \n");
  ##% Algebraic operations
  printf("* Algebraic procedures\n");
  #% ContractT
  printf("  + ContractT(T::tensor table,[a,b],Gdata[,sw]) => CT::tensor table \n");
  printf("    Constructing a tensor table CT with two rank down from a tensor table T by contracting the specified indecies.\n")+
  printf("    [a,b] = the index positions to contract \n");  
  #% IP
  printf("  + IP(V,W, Gdata) => val::scalar \n");
  printf("    calculates the inner product of two tensors/vectors/forms of the same rank.\n");
  printf("    V,W:: tensor tables/vectors/dforms\n");
  
end proc:


######>Applications
#####>Job control
Job_KN := false:
#####> Kerr-Newmann spacetime 
if  Job_KN then 
    print("Kerr-Newmann spacetime");
###> metric & curvature 
     m:= 'm';  a:= 'a': e:='e':
     Delta := r^2 - 2*m*r+a^2+ e^2:
     Sigma2:= r^2+a^2*cos(theta)^2:
     Gamm:= (r^2+a^2)^2 -a^2*Delta*sin(theta)^2:
     Omega := a*(2*m*r-e^2)/Gamm:
    ds2_KN :=  -Delta*Sigma2/Gamm*dt^2 + Gamm*sin(theta)^2/Sigma2 *(dphi-Omega*dt)^2
                       + Sigma2*(dr^2/Delta + dtheta^2):
print("ds2"=ds2_KN):
    varlist_S := [ t, r, theta, phi]:
    Gdata_KN := Riemann(varlist_S, ds2_KN,1):
    dvarlist_S:= Gdata_KN["dvarlist"] :
    Ric11_KN := Gdata_KN["Ric11"]:
    Ric11_KN["indextype"]=[_u,_d]:
    Riem13_KN := Gdata_KN["Riem13"]:  #* tensordict *#
    Riem13_KN["indextype"]:=[_u,_d, _d,_d]:
    Weyl04_KN := Gdata_KN["Weyl04"]:  #* tensordict *#
    Weyl04_KN["indextype"]:=[_d,_d, _d,_d]:
    print("Calculation of curvatures:= done"):
###> Orthonormal basis 
    ONFB_KN:=  [ FB[1]= sqrt(Delta*Sigma2/Gamma)*dt, 
                          FB[2]=sqrt(Gamma/Sigma2)*sin(theta)*(dphi-Omega*dt), 
                          FB[3]=sqrt(Sigma2/Delta)*dr,
                          FB[4]= sqrt(Sigma2)*dtheta 
                        ]:
    Riem13_KN_ONB := CBtoFBT(Riem13_KN, varlist_S, ONFB_KN):
###> Null basis 
    NullFB_KN:= [ NFB[2] = -dt + a*sin(theta)^2*dphi+ Sigma2/Delta*dr, 
                           NFB[1]= Delta/(2*Sigma2)*(-dt+a*sin(theta)^2*dphi)-dr/2, 
                           NFB[4]=  1/sqrt(2*Sigma2)*(Sigma2*dtheta-I*a*sin(theta)*dt + I*(r^2+a^2)*sin(theta)*dphi),
                           NFB[3]= 1/sqrt(2*Sigma2)*(Sigma2*dtheta + I*a*sin(theta)*dt - I*(r^2+a^2)*sin(theta)*dphi)
                        ]:
    Weyl04_KN_ONB := CBtoFBT(Weyl04_KN, varlist_S, NullFB_KN):
###> NP coefficients 
    NPC_KN := NPC(Gdata_KN, NullFB_KN,1):
end if: 
#end:= Job_KN#

######>EOF