# James Shifflett
#
# version: 17Feb2007
#
# This program checks that the LRES charged solution satisfies the
# LRES field equations. It also gives the tetrad representation
# and calculates spin coefficients and some other things.
# It uses geometrized units where F_ik has units of 1/cm rather than
# the special units in some of my papers where F_ik is unitless.

with(tensor):
DIM:=4:
coord:=[t,r,theta,phi];

`tensor/simp`:=proc(x) local ANS,D0; global theta,r;
 ANS:=radsimp(expand(normal(x)));
end:

CH:=ch(r):
SH:=sh(r):
AA:=aa(r):
VH:=vh(r):
VHAT:=vhat(r):
F01:=f01(r):
A0:=a0(r):

print("various functions used in the solution");
sh(r):=(Q/r^2)*2^(1/2)*I/L^(1/2);
ch(r):=(1+SH^2)^(1/2);
vh(r):=(r*L/Q^2)*(Int(r^2*CH,r)-r^3/3);
vhat(r):=(r*L/Q^2)*(int(r^2*CH,r)-r^3/3);
aa(r):=1-2*M/r-Lambda*r^2/3+(Q^2*VH/r^2)*(1-Lambda/L);
f01(r):=(Q/r^2)*(1+2*diff(AA,r)/(L*r));
A0(r):=-int(f01,r);
re:=Q^(1/2)*2^(1/4)/L^(1/4);

print("the tetrads");
h:=create([1,-1],array(1..DIM,1..DIM, \
[[AA*CH/2,1/2,0,0], \
 [1,-1/(AA*CH),0,0], \
 [0,0,r*(CH/2)^(1/2),-I*r*sin(theta)*(CH/2)^(1/2)],\
 [0,0,r*(CH/2)^(1/2), I*r*sin(theta)*(CH/2)^(1/2)]]));

print("inverse tetrads");
hi:=invert(h,'deth');

# help-npspin specifies the form of h, help-connexF doesn't make sense
# for a contravariant or covariant tensor, invert(h)=inverse(h)
# for a mixed tensor, invert(h)=transpose(inverse(h)), see help-invert
#
#  (a)        j               j
# h       hi       dhi   =hi   ,
#    j      (a)       (a)   (a) j
#  
# dhi can be used to take derivatives easily
#
#   (a)(b)        (a)(b)         i     (a)i      (a)(b)
#  N      ,    = N      ,  * hi    =  N    ,  - N      dhi
#          (b)           i     (b)          i             (b)
#
# This method is required in some cases because the student version has a
# limit on the number of indices that tensors can have (3?)

# calculate spin coefficients
SPN_CF:=npspin(coord,h,'G','any'):

print("calculate curvature and Weyl tensor components");
CURVE:=npcurve(SPN_CF,any);

print("covariant metric in tetrad form");
gtetrad:=create([-1,-1],array(symmetric,1..DIM,1..DIM, \
                                        [[ 0, 1, 0, 0] \
                                        ,[ 1, 0, 0, 0] \
                                        ,[ 0, 0, 0,-1] \
                                        ,[ 0, 0,-1, 0]])) ;
gitetrad:=invert(gtetrad,'detgtetrad'):

print("fundamental tensor in tetrad form");
Ntetrad:=create([-1,-1],array(1..DIM,1..DIM, \
          [[  0  ,CH-SH, 0    , 0    ] \
          ,[CH+SH,  0  , 0    , 0    ] \
          ,[  0  ,  0  , 0    ,-1/CH ] \
          ,[  0  ,  0  ,-1/CH , 0    ]]));
print("inverse fundamental tensor in tetrad form");
NtetradI:=invert(Ntetrad,'detNtetrad');

# calculate spin coeffiecent array

print("check that Maple's two methods of calculating spin coefficients agree");

# for some bizarre reason the two do not check if I substitute Lb for L !!!

gamtetrad:=connexF(coord,gtetrad,h):
#Rm:=RiemannF(coord,gitetrad,hi,gamtetrad);
# connexF IS OPPOSITE SIGN FROM CHANDRASECHAR, SPIN COEFFICIENTS ARE THE SAME
gamtetrad:=lin_com(-1,gamtetrad):
expand(gamtetrad[compts][3,1,1]-SPN_CF[kappa]);
expand(gamtetrad[compts][2,4,4]-SPN_CF[lambda]);
expand(gamtetrad[compts][3,1,3]-SPN_CF[sigma]);
expand(gamtetrad[compts][2,4,2]-SPN_CF[nu]);
expand(gamtetrad[compts][3,1,4]-SPN_CF[rho]);
expand(gamtetrad[compts][2,4,3]-SPN_CF[mu]);
expand(gamtetrad[compts][2,4,1]-SPN_CF[pi]);
expand(gamtetrad[compts][3,1,2]-SPN_CF[tau]);
expand((gamtetrad[compts][2,1,1]+gamtetrad[compts][3,4,1])/2-SPN_CF[epsilon]);
expand((gamtetrad[compts][2,1,2]+gamtetrad[compts][3,4,2])/2-SPN_CF[gamma]);
expand((gamtetrad[compts][2,1,4]+gamtetrad[compts][3,4,4])/2-SPN_CF[alpha]);
expand((gamtetrad[compts][2,1,3]+gamtetrad[compts][3,4,3])/2-SPN_CF[beta]);

print("covariant metric g");
g:=create([-1,-1],array(1..DIM,1..DIM, \
    [[ CH*AA ,     0      ,    0    ,             0         ] \
    ,[   0   , -1/(CH*AA) ,    0    ,             0         ] \
    ,[   0   ,     0      , -CH*r^2 ,             0         ] \
    ,[   0   ,     0      ,    0    ,-CH*r^2*(sin(theta))^2 ]]));
print("contravariant metric, gI");
gI:=invert(g,'detg');

print("covariant fundamental tensor, N");
N:=create([-1,-1],array(1..DIM,1..DIM, \
    [[ AA*CH^2 ,   SH  ,   0  ,          0         ] \
    ,[   -SH   , -1/AA ,   0  ,          0         ] \
    ,[    0    ,   0   , -r^2 ,          0         ] \
    ,[    0    ,   0   ,   0  ,-r^2*(sin(theta))^2 ]]));
print("inverse fundamental tensor, NI");
NI:=invert(N,'detN');

print("check that g=h^T*gtetrad*h and N=h^T*Ntetrad*h");
for js from 1 to DIM do \
for jm from 1 to DIM do \
og:=g[compts][js,jm];
oN:=N[compts][js,jm];
for jb from 1 to DIM do \
for ja from 1 to DIM do \
og:=og-h[compts][ja,js]*gtetrad[compts][ja,jb]*h[compts][jb,jm];
oN:=oN-h[compts][ja,js]*Ntetrad[compts][ja,jb]*h[compts][jb,jm];
od: \
od: \
og:=simplify(og);
oN:=simplify(oN);
print('og and oN',js,jm,"=",og,oN);
od: \
od:

print("calculate the symmetric connections, gam");
gam:=create([1,-1,-1],array(1..DIM,1..DIM,1..DIM)):
for ja from 1 to DIM do \
for js from 1 to DIM do \
for jm from 1 to DIM do \
gam[compts][ja,js,jm]:=0;
for jn from 1 to DIM do \
gam[compts][ja,js,jm]:=gam[compts][ja,js,jm] \
                      +gI[compts][ja,jn]*(diff(g[compts][jm,jn],coord[js])  \
                                         +diff(g[compts][jn,js],coord[jm])  \
                                         -diff(g[compts][js,jm],coord[jn]));
od: \
od: \
od: \
od:

print("calculate the symmetric Ricci tensor, Rs");
Rs:=create([-1,-1],array(1..DIM,1..DIM)):
for jn from 1 to DIM do \
for jm from 1 to DIM do \
Rs[compts][jn,jm]:=0;
for ja from 1 to DIM do \
Rs[compts][jn,jm]:=Rs[compts][jn,jm]+diff(gam[compts][ja,jn,jm],coord[ja]) \
                                    -diff(gam[compts][ja,jn,ja],coord[jm]);
for js from 1 to DIM do \
Rs[compts][jn,jm]:=Rs[compts][jn,jm]+gam[compts][js,jn,jm]*gam[compts][ja,js,ja] \
                                    -gam[compts][js,jn,ja]*gam[compts][ja,js,jm]; 
od: \
od: \
Rs[compts][jn,jm]:=simplify(Rs[compts][jn,jm]);
Rs[compts][jn,jm]:=radsimp(Rs[compts][jn,jm]);
# print("Rs",jn,jm,"=",Rs[compts][jn,jm]);
od: \
od:

# x:=simplify(radsimp(CH^2*subs(Lambda=0,Rs[compts][2,2]))):
# x:=subs(diff(CH,r)=diff(SH,r)*SH/CH,x):
# x:=subs(diff(CH,r)=diff(SH,r)*SH/CH,x):
# x:=radsimp(subs(Int(r^2*CH,r)=y*Q^(3/2)*(2/L)^(3/4),x)):
# x:=simplify(radsimp(subs(ch(r)=0,x))):
# x:=simplify(radsimp(subs(r=sqrt(Q)*2^(1/4)/L^(1/4),x))):
# M1:=(Q^(1/2)*2^(1/4)*1+Q^(3/2)*2^(3/4)*L^(1/2)*(3*y+1))/(6*L^(1/4)):
# M2:=(Q^(1/2)*2^(1/4)*3+Q^(3/2)*2^(3/4)*L^(1/2)*(3*y-1))/(6*L^(1/4)):
# z:=expand(numer(x)-9*36*sqrt(2)*L*(M-M1)*(M-M2)):

gamt:=create([1,-1,-1],array(1..DIM,1..DIM,1..DIM,
      [[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]] \
      ,[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]] \
      ,[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]] \
      ,[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]])):

print("the nonsymmetric connections, gamt");
# note that t,r,theta,phi=1,2,3,4 not 0,1,2,3
gamt[compts][2,1,1]:=diff(AA,r)*AA*CH^2/2+2*AA^2*SH^2/r;
gamt[compts][1,2,1]:=diff(AA,r)/(2*AA);
gamt[compts][1,1,2]:=diff(AA,r)/(2*AA);
gamt[compts][2,2,2]:=-diff(AA,r)/(2*AA);
gamt[compts][3,2,3]:=1/r;
gamt[compts][3,3,2]:=1/r;
gamt[compts][4,2,4]:=1/r;
gamt[compts][4,4,2]:=1/r;
gamt[compts][2,3,3]:=-AA*r;
gamt[compts][2,4,4]:=-AA*r*(sin(theta))^2;
gamt[compts][4,3,4]:=cot(theta);
gamt[compts][4,4,3]:=cot(theta);
gamt[compts][3,4,4]:=-sin(theta)*cos(theta);
gamt[compts][3,1,3]:=-AA*SH/r;
gamt[compts][3,3,1]:= AA*SH/r;
gamt[compts][4,1,4]:=-AA*SH/r;
gamt[compts][4,4,1]:= AA*SH/r;
gamt[compts][2,2,1]:=-2*AA*SH/r;
gamt[compts][2,1,2]:= 2*AA*SH/r;

F:=create([-1,-1],array(1..DIM,1..DIM, \
           [[  0 ,F01,0,0] \
           ,[-F01, 0 ,0,0] \
           ,[  0 , 0 ,0,0] \
           ,[  0 , 0 ,0,0]])):

print("calculate the nonsymmetric Ricci tensor, Rt");
Rt:=create([-1,-1],array(1..DIM,1..DIM)):
for jn from 1 to DIM do \
for jm from 1 to DIM do \
Rt[compts][jn,jm]:=0;
for ja from 1 to DIM do \
Rt[compts][jn,jm]:=Rt[compts][jn,jm]+diff(gamt[compts][ja,jn,jm],coord[ja]) \
                                    -diff(gamt[compts][ja,ja,jn],coord[jm])/2 \
                                    -diff(gamt[compts][ja,ja,jm],coord[jn])/2;
for js from 1 to DIM do \
Rt[compts][jn,jm]:=Rt[compts][jn,jm]+gamt[compts][js,jn,jm]*gamt[compts][ja,js,ja] \
                                    -gamt[compts][js,jn,ja]*gamt[compts][ja,js,jm]; 
od: \
od: \
Rt[compts][jn,jm]:=simplify(Rt[compts][jn,jm]);
print("Rt",jn,jm,"=",Rt[compts][jn,jm]);
od: \
od:

print("check that Ampere's law is satisfied");
for jn from 1 to DIM do \
oj:=0;
for jm from 1 to DIM do \
oj:=oj+diff(sqrt(-detN)*(NI[compts][jn,jm]-NI[compts][jm,jn]),coord[jm]);
od: \
oj:=simplify(oj);
print(jn,"=",oj);
od:

print("check that connection equations are satisfied");
ogam:=create([-1],array(1..DIM)):
for jn from 1 to DIM do \
for jm from 1 to DIM do \
for jb from 1 to DIM do \
ogam[jb]:=diff(N[compts][jn,jm],coord[jb]);
for ja from 1 to DIM do \
ogam[jb]:=ogam[jb]-gamt[compts][ja,jn,jb]*N[compts][ja,jm] \
                  -gamt[compts][ja,jb,jm]*N[compts][jn,ja]; 
od: \
ogam[jb]:=simplify(ogam[jb]);
od: \
print(jn,jm,"=",ogam[1],ogam[2],ogam[3],ogam[4]);
od: \
od:

print("check that Einstein equation are satisfied");
for jn from 1 to DIM do \
for jm from 1 to DIM do \
oRt:=Rt[compts][jn,jm] \
    -F[compts][jn,jm]*I*sqrt(2*L) \
    +L*N[compts][jn,jm] \
    +(Lambda-L)*g[compts][jn,jm];
oRt:=simplify(oRt);
print(jn,jm,"=",oRt);
od: \
od: