###############################################################################
#
#  bl.ms:
# -------
#   -- This Maple script verifies that the "shift_phi", "beta_phi" and
#      "lapse" variables are derived correctly in the Einstein Toolkit file:
#
#          EinsteinInitialData/IDAnalyticBH/src/Kerr.c
#
#
#   Written by Scott C. Noble  Thu Oct  7 11:06:34 EDT 2010
#
#
###############################################################################

with(linalg);

################################################################################
# EinsteinInitialData/IDAnalyticBH/src/Kerr.c Notation:
#  -- I copied lines of the code here and then converted its C syntax
#     to Maple.
#  -- Also, I commented out some definitions that were not necessary
#     for the comparison.
#  -- Note that  "rK" here is the typical "r" radius of
#     Boyer-Lindquist coordinates,  and "R" is the radius variable of
#     isotropic coordinates.
#
################################################################################

#  CCTK_REAL m=mass,a=a_Kerr,a_2=a*a,m2_a2=m*m-a_2;
   m:=M;
  rK:=r;

    a_2:=a*a;
    m2_a2:=m*m-a_2;


#    xx=x[i];  yy=y[i];  zz=z[i];
#    rho_2=xx*xx+yy*yy;
#    rho=sqrt(rho_2);
#    R_2=rho_2+zz*zz;
    R_2 := R*R;
#    R=sqrt(R_2);
    R_3:=R*R_2;
#    cth=zz/R;
    cth := cos(th);
    cth_2:=cth*cth;
#    sth_2=rho_2/R_2;
#    sth=rho/R;
    sth := sin(th);
    sth_2 := sth*sth;

#    rK:=R+m+m2_a2/4/R;
   solR:=solve(rK=R+m+m2_a2/4/R,R);
   R:=solR[1];
#   R:=solR[2];
    r_2:=rK*rK;
    sqrt_Delta:=R-m2_a2/4/R;
    Sigma:=r_2+a_2*cth_2;
    # \beta_\phi
    beta_phi:=-2*m*rK*a*sth_2/Sigma;
    p2:=a_2+r_2-a*beta_phi;
   # /* drdR=sqrt_Delta/R; */
    lapse:=sqrt_Delta/sqrt(p2);
   # /* shift_phi=-2*m*rK*a/p2; */
    # \beta^\phi

    shift_phi:=-2*m*rK*a/(p2*Sigma);  # Proposed "correct" expression
    shift_phi_old:=-2*m*rK*a/(p2);    # "Wrong" expression from original version of Kerr.c


################################################################################
# Following  Shapiro & Teukolsky, "Black Holes, White Dwarfs....", (1983).
################################################################################

Delta:=r^2-2*M*r+a^2;
Sigma2:=r^2+a^2*cos(th)^2;
br := 2*M*r/Sigma2;
A:= r^2 + a^2 * ( 1 + br*sin(th)^2 );

gcon := array(1..4,1..4);
gcov := array(1..4,1..4);
betaup := array(1..4);
betadn := array(1..4);

##############################
#### Boyer-Lindquist:

gcov[1,1] := -(1-br);
gcov[1,2] := 0;
gcov[1,3] := 0;
gcov[1,4] := -br*a*sin(th)^2;
gcov[2,1] := gcov[1,2];
gcov[2,2] := Sigma2/Delta;
gcov[2,3] := 0;
gcov[2,4] := 0;
gcov[3,1] := gcov[1,3];
gcov[3,2] := gcov[2,3];
gcov[3,3] := Sigma2;
gcov[3,4] := 0;
gcov[4,1] := gcov[1,4];
gcov[4,2] := gcov[2,4];
gcov[4,3] := gcov[3,4];
gcov[4,4] := sin(th)^2*A;

# This was derived by hand by myself:
gcon[1,1] := - A / Delta;
gcon[1,2] := 0 ;
gcon[1,3] := 0 ;
gcon[1,4] := -br*a/Delta;
gcon[2,1] := 0;
gcon[2,2] := Delta/Sigma2;
gcon[2,3] := 0;
gcon[2,4] := 0;
gcon[3,1] := 0;
gcon[3,2] := 0;
gcon[3,3] := 1/Sigma2;
gcon[3,4] := 0;
gcon[4,1] := gcon[1,4];
gcon[4,2] := gcon[2,4];
gcon[4,3] := gcon[3,4];
gcon[4,4] := (1-br)/(Delta*sin(th)^2);

#### TESTS:
#verify my gcon calculation:
gcontest := map(simplify, evalm(inverse(gcov) - gcon));
gcontest2 := map( simplify, evalm(gcon &* gcov));


##############################
# ADM functions:
##############################

# lapse :
alpha := sqrt(simplify(-1/gcon[1,1]));

# \beta^i
for i from 1 by 1 to 4 do
    betadn[i]:=0;
    betaup[i]:=0;
end do;

for i from 2 by 1 to 4 do
    betaup[i] := alpha^2*gcon[1,i];
end do;

for i from 2 by 1 to 4 do
    for j from 2 by 1 to 4 do
        betadn[i] := betadn[i]  + gcov[i,j]*betaup[j];
    end do;
end do;

# \beta_\phi :
beta_phi2 := -2*M*r*a*(sin(th))^2 / Sigma2;
print("This is the difference between my closed form expression for beta_\phi and the Maple-derived quantity:");
check_beta_phi_loc := simplify(beta_phi2 -betadn[4]);
print("This is the difference between my derived expression for beta_\phi and the expression in Kerr.c:");
check_beta_phi_code := simplify(betadn[4] - beta_phi);

# \beta^\phi :
shift_phi2 := -2*M*r*a/(Sigma2*A);
print("This is the difference between my closed form expression for beta^\phi and the Maple-derived quantity:");
check_shift_phi_loc := simplify(shift_phi2 - betaup[4]);
print("This is the difference between my derived expression for beta^\phi and the corrected expression in Kerr.c:");
check_shift_phi_code := simplify(betaup[4] - shift_phi);
print("This is the difference between my derived expression for beta^\phi and the original (wrong) expression in Kerr.c:");
check_shift_phi_old_code := simplify(betaup[4] - shift_phi_old);

# lapse:
print("This is the difference between my derived expression for the lapse and the one in Kerr.c:");
check_lapse := simplify(alpha^2 - lapse^2);