# James Shifflett
#
# version: 27Feb2011
#
# This program checks that the LRES electromagnetic plane-wave solution
# satisfies the LRES field equations. The solution is checked in pp-wave
# x,y,u,v coordinates as in Stephani,Kramer,MacCallum,Hoenselaers,
# but the solution is also transformed to x,y,z,t coordinates
# and to the alternative xyuv coordinates in Misner,Thorne & Wheeler, p.961.
# The program also contains a routine to calculate connections for null fields.
# 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:
gam:=array(1..DIM,1..DIM,1..DIM):
k:=array(1..DIM):
j:=array(1..DIM):
gradA3:=array(1..DIM):
uh:=array(1..DIM,1..DIM,1..DIM):
ub:=array(1..DIM,1..DIM,1..DIM):
cuh:=array(1..DIM):
cuh1:=array(1..DIM):
uh1:=array(1..DIM,1..DIM,1..DIM):
gamt:=array(1..DIM,1..DIM,1..DIM):

print("x=x, y=y, u=(-z+t)/sqrt(2), v=( z+t)/sqrt(2)");
coord:=[x,y,u,v];

HP:=hp(u):
HC:=hc(u):
FC:=fc(u):
FCX:=fcx(u):
FCY:=fcy(u):
L:=l(u):
HH:=H(x,y,u):

star:=proc(arg) local asdf;
asdf:=subs(I=(-I),arg);
asdf
end:

re:=proc(arg) local asdf;
asdf:=subs(I=0,arg);
asdf
end:

print("wavenumber");
k:=[0,0,-1,0];

A3:=-sqrt(2)*(FCX*x+FCY*y);

print("the electromagnetic potential");
A:=create([-1],array(1..DIM,[0,0,A3,0]));

print("covariant metric g in xyuv coordinates");
gll:=create([-1,-1],array(1..DIM,1..DIM, \
        [[-1 , 0 ,  0  , 0 ] \
        ,[ 0 ,-1 ,  0  , 0 ] \
        ,[ 0 , 0 ,  HH , 1 ] \
        ,[ 0 , 0 ,  1  , 0 ]]));

print("contravariant metric g in xyuv coordinates");
grr:=invert(gll,'detgll');
rmg:=sqrt(-detgll);

FFll:=create([-1,-1],array(1..DIM,1..DIM)):
for JB from 1 to DIM do \
for JC from 1 to DIM do \
FFll[compts][JB,JC]:=diff(A[compts][JC],coord[JB])-diff(A[compts][JB],coord[JC]); \
od: \
od;

for JB from 1 to DIM do \
gradA3[JB]:=diff(A3,coord[JB]): \
od:

fll:=create([-1,-1],array(1..DIM,1..DIM)):
for JB from 1 to DIM do \
for JC from 1 to DIM do \
fll[compts][JB,JC]:=k[JB]*gradA3[JC]-k[JC]*gradA3[JB]; \
od: \
od;

print("electromagnetic field F calculated from curl of A");
eval(FFll);
print("electromagnetic field F calculated from A3 and k");
eval(fll);
FFlr:=prod(FFll,grr,[2,1]):
print("raised electromagnetic field density F");
FFrr:=lin_com(rmg,prod(grr,FFlr,[2,1]));

print("transformation to xyzt coordinates");
cxyzt:=[x,y,(v-u)/sqrt(2),(v+u)/sqrt(2)];
XYZT:=create([1,-1],array(1..DIM,1..DIM)):
for JB from 1 to DIM do \
for JC from 1 to DIM do \
XYZT[compts][JB,JC]:=diff(cxyzt[JB],coord[JC]); \
od: \
od;
print("The transformation matrix is");
eval(XYZT);
print("The MAPLE invert function really gives inverse transpose");
XYZTI:=invert(XYZT,'detXYZT');
print("covariant g, A and F in xyzt coordinates");
gllb:=prod(prod(XYZTI,gll,[2,1]),XYZTI,[2,2]);
Ab:=prod(XYZTI,A,[2,1]);
FFllb:=prod(prod(XYZTI,FFll,[2,1]),XYZTI,[2,2]);
print("from Landau-Lifshitz p61 covariant F in xyzt coordinates is");
print("( 0  -Hz  Hy -Ex )");
print("( Hz  0  -Hx -Ey )");
print("(-Hy  Hx  0  -Ez )");
print("( Ex  Ey  Ez  0  )");

print("transformation to Misner,Thorne,Wheeler x,y,u,v coordinates (MTW p.961)");
cMTW:=[x/L,y/L,u*sqrt(2),v*sqrt(2)-(x^2+y^2)*diff(L,u)/(L*sqrt(2))];
MTW:=create([1,-1],array(1..DIM,1..DIM)):
for JB from 1 to DIM do \
for JC from 1 to DIM do \
MTW[compts][JB,JC]:=diff(cMTW[JB],coord[JC]); \
od: \
od;
print("The transformation matrix is");
eval(MTW);
print("The MAPLE invert function really gives inverse transpose");
MTWI:=invert(MTW,'detMTW');
print("covariant g, A and F in Misner,Thorne,Wheeler x,y,u,v coordinates");
gllp:=prod(prod(MTWI,gll,[2,1]),MTWI,[2,2]);
Ap:=prod(MTWI,A,[2,1]);
FFllp:=prod(prod(MTWI,FFll,[2,1]),MTWI,[2,2]);


# change f into ftilde
fll:=lin_com(I*sqrt(2)/sqrt(Lb),fll):
frl:=prod(grr,fll,[2,1]):
frr:=prod(frl,grr,[2,1]):
NI:=lin_com(1,grr,-1,frr):
N:=invert(NI,'dNI'):
rmN:=sqrt(-dNI):
print("fundamental tensor, N");
N:=lin_com(rmN,N);

print("Christoffel connections, gam");
for JA from 1 to DIM do \
for JB from 1 to DIM do \
for JC from 1 to DIM do \
gam[JA,JB,JC]:=0; \
for JD from 1 to DIM do \
gam[JA,JB,JC]:=gam[JA,JB,JC]+grr[compts][JA,JD]*( diff(gll[compts][JC,JD],coord[JB]) \
                                                 +diff(gll[compts][JD,JB],coord[JC]) \
                                                 -diff(gll[compts][JB,JC],coord[JD]))/2; \
od: \
if not(gam[JA,JB,JC] = 0) then print("gam",JA,JB,JC,"=",gam[JA,JB,JC]); fi; \
od: \
od: \
od;

print("calculate connections for null fields");

# covariant uh1, the order f solution
for JA from 1 to DIM do \
for JS from 1 to DIM do \
for JM from 1 to DIM do \
uh1[JA,JS,JM]:=diff(fll[compts][JS,JM],coord[JA]) \
              +diff(fll[compts][JA,JM],coord[JS]) \
              -diff(fll[compts][JA,JS],coord[JM]); \
for JN from 1 to DIM do \
uh1[JA,JS,JM]:=uh1[JA,JS,JM]-gam[JN,JS,JA]*fll[compts][JN,JM]-gam[JN,JM,JA]*fll[compts][JS,JN] \
                            -gam[JN,JA,JS]*fll[compts][JN,JM]-gam[JN,JM,JS]*fll[compts][JA,JN] \
                            +gam[JN,JA,JM]*fll[compts][JN,JS]+gam[JN,JS,JM]*fll[compts][JA,JN]; \
od: \
uh1[JA,JS,JM]:=uh1[JA,JS,JM]/2; \
od: \
od: \
od;

for JA from 1 to DIM do \
cuh1[JA]:=0: \
for JT from 1 to DIM do \
for JN from 1 to DIM do \
cuh1[JA]:=cuh1[JA]+uh1[JN,JT,JA]*frr[compts][JT,JN]; \
od: \
od: \
od:

# covariant uh
for JA from 1 to DIM do \
for JS from 1 to DIM do \
for JM from 1 to DIM do \
uh[JA,JS,JM]:=uh1[JA,JS,JM]+cuh1[JA]*fll[compts][JM,JS]/2; \
for JN from 1 to DIM do \
uh[JA,JS,JM]:=uh[JA,JS,JM]+( cuh1[JN]*frl[compts][JN,JM]*gll[compts][JS,JA] \
                            -cuh1[JN]*frl[compts][JN,JS]*gll[compts][JM,JA])/2; \
for JT from 1 to DIM do \
uh[JA,JS,JM]:=uh[JA,JS,JM]+uh1[JA,JN,JT]*frl[compts][JN,JM]*frl[compts][JT,JS] \
                         +(uh1[JM,JN,JT]+uh1[JN,JM,JT])*frl[compts][JN,JA]*frl[compts][JT,JS]/2 \
                         -(uh1[JS,JN,JT]+uh1[JN,JS,JT])*frl[compts][JN,JA]*frl[compts][JT,JM]/2 \
                         +(uh1[JS,JM,JN]-uh1[JM,JS,JN])*frl[compts][JN,JT]*frl[compts][JT,JA]/2; \
od: \
od: \
od: \
od: \
od;

for JA from 1 to DIM do \
cuh[JA]:=0: \
for JT from 1 to DIM do \
for JN from 1 to DIM do \
cuh[JA]:=cuh[JA]+uh[JN,JT,JA]*frr[compts][JT,JN]; \
od: \
od: \
od:

# covariant ub
for JA from 1 to DIM do \
for JS from 1 to DIM do \
for JM from 1 to DIM do \
ub[JA,JS,JM]:=-cuh[JA]*gll[compts][JM,JS]/2+(cuh[JM]*gll[compts][JS,JA]+cuh[JS]*gll[compts][JM,JA])/2; \
for JT from 1 to DIM do \
ub[JA,JS,JM]:=ub[JA,JS,JM]+(uh[JA,JM,JT]-uh[JM,JA,JT])*frl[compts][JT,JS]/2 \
                          +(uh[JA,JS,JT]-uh[JS,JA,JT])*frl[compts][JT,JM]/2 \
                          +(uh[JS,JM,JT]+uh[JM,JS,JT])*frl[compts][JT,JA]/2; \
od: \
od: \
od: \
od;

print("nonsymmetric connections, gamt");
for JA from 1 to DIM do \
for JS from 1 to DIM do \
for JM from 1 to DIM do \
gamt[JA,JS,JM]:=gam[JA,JS,JM]; \
for JT from 1 to DIM do \
gamt[JA,JS,JM]:=gamt[JA,JS,JM]+grr[compts][JA,JT]*(ub[JT,JS,JM]+uh[JT,JS,JM]); \
od: \
if not(gamt[JA,JS,JM] = 0) then print("gamt",JA,JS,JM,"=",gamt[JA,JS,JM]); fi; \
od: \
od: \
od;

print("symmetric Ricci tensor, R");
R:=create([-1,-1],array(1..DIM,1..DIM)):
for JS from 1 to DIM do \
for JM from 1 to DIM do \
R[compts][JS,JM]:=0; \
for JA from 1 to DIM do \
R[compts][JS,JM]:=R[compts][JS,JM]+diff(gam[JA,JS,JM],coord[JA])-diff(gam[JA,JS,JA],coord[JM]); \
for JN from 1 to DIM do \
R[compts][JS,JM]:=R[compts][JS,JM]+gam[JN,JS,JM]*gam[JA,JN,JA]-gam[JN,JS,JA]*gam[JA,JN,JM]; \
od: \
od: \
od: \
od;
eval(R);

print("nonsymmetric Ricci tensor, Rt");
Rt:=create([-1,-1],array(1..DIM,1..DIM)):
for JS from 1 to DIM do \
for JM from 1 to DIM do \
Rt[compts][JS,JM]:=0; \
for JA from 1 to DIM do \
Rt[compts][JS,JM]:=Rt[compts][JS,JM]+diff(gamt[JA,JS,JM],coord[JA])-diff(gamt[JA,JS,JA],coord[JM]); \
for JN from 1 to DIM do \
Rt[compts][JS,JM]:=Rt[compts][JS,JM]+gamt[JN,JS,JM]*gamt[JA,JN,JA]-gamt[JN,JS,JA]*gamt[JA,JN,JM]; \
od: \
od: \
od: \
od;
eval(Rt);

for JB from 1 to DIM do \
j[JB]:=0: \
for JC from 1 to DIM do \
j[JB]:=j[JB]+diff(rmg*frr[compts][JB,JC],coord[JC]); \
od: \
od:

print("Einstein-Maxwell Einstein equations with unspecified H");
lin_com(1,R,-2,prod(FFlr,FFll,[2,1]),prod(FFlr,FFlr,[2,1],[1,2])[compts][1]/2,gll);

print("LRES Einstein equations with unspecified H");
lin_com(1,Rt,-sqrt(2*Lb)*I,FFll,Lb,N,-Lb,gll);

print("check that gamt satisfies connection equations");
oo:=array(1..DIM):
for JA from 1 to DIM do \
for JS from 1 to DIM do \
for JM from 1 to DIM do \
oo[JM]:=diff(N[compts][JA,JS],coord[JM]); \
for JT from 1 to DIM do \
oo[JM]:=oo[JM]-gamt[JT,JA,JM]*N[compts][JT,JS]-gamt[JT,JM,JS]*N[compts][JA,JT]; \
od: \
oo[JM]:=simplify(oo[JM]); \
od: \
print(JA,JS," ",oo[1],oo[2],oo[3],oo[4]);
od: \
od;

H(x,y,u):=2*(HP*x^2+HC*x*y-HP*y^2)+2*(FCX^2+FCY^2)*(x^2+y^2);

print("check Einstein-Maxwell Einstein equations");
lin_com(1,R,-2,prod(FFlr,FFll,[2,1]),prod(FFlr,FFlr,[2,1],[1,2])[compts][1]/2,gll);

print("check LRES Einstein equations");
lin_com(1,Rt,-sqrt(2*Lb)*I,FFll,Lb,N,-Lb,gll);

print("check Ampere's law");
eval(j);