/*
Some basics ...
x_s(t) = position of the object at time t ... I'm guessing this is coordinate time, and not local time, which means x_s is only a function of coordinate time (no two positions for one point in time) which means we're not considering things going backwards in time
v_s(t) = d/dt x_s(t) = the velocity wrt coordinate time
r_s(t,x,y,z) = sqrt( (x - x_s(t))^2 + y^2 + z^2) = distance from position of object at time t
f(r_s) = ( tanh(sigma*(r_s + R)) - tanh(sigma*(r_s - R)) )/( 2*tanh(sigma*R) )
R = warp bubble radius
sigma = warp bubble thickness

ADM:
alpha = 1
gamma_ij = delta_ij
beta^x = -v_s(t) * f(r_s(t))
beta^y = beta^z = 0

u = v_s * f(r_s)
... v_s(t), r_s(t,x,y,z)

u,_u = v_s_,u f(r_s) + v_s f'(r_s) r_s_,u
v_s_,t = x_s_,tt
r_s_,t = -v_s*(x - x_s)/r_s
r_s_,x = (x - x_s)/r_s
r_s_,y = y/r_s
r_s_,z = z/r_s
*/;

load(ctensor);

dim : 4;

ct_coords : [t,x,y,z];

cframe_flag : false;

depends(u, [t,x,y,z]);

depends(x_s, t);

v_s_def : diff(x_s, t);

f_def(r_s) := ( tanh(sigma*(r_s + R)) - tanh(sigma*(r_s - R)) )/( 2*tanh(sigma*R) );

depends(v_s, t);

depends(r_s, [t,x,y,z]);

depends(f, [r_s]);

u_from_v_s_f : u = v_s * f;

r_s_def(t,x,y,z) := sqrt( (x - x_s(t))^2 + y^2 + z^2  );

u_def(t,x,y,z) := v_s(t) * f_def(r_s_def(t,x,y,z));

diff(u_def(t,x,y,z), t)$
ratsimp(%);

diff(u_def(t,x,y,z), x)$
ratsimp(%);

lg : matrix(
    [u^2-1, -u, 0, 0],
    [-u, 1-u^2, 0, 0],
    [0, 0, 1, 0],
    [0, 0, 0, 1]
);

cmetric()$

depends(g, [t,x,y,z]);

gdet;

ug$
ratsimp(%)$
subst([-gdet = -g], %)$
ug : %;

christof(false)$

for i:1 while i<=4 do (
    for j:1 while j<=4 do (
        for k:1 while k<=4 do (
            mcs[i,j,k] : ratsimp(mcs[i,j,k]),
            lcs[i,j,k] : ratsimp(lcs[i,j,k])
        )
    )
)$

cdisplay(lcs)$

cdisplay(mcs)$

ricci(false)$

for i:1 while i<=4 do (
    for j:1 while j<=4 do (
        ric[i,j] : ratsimp(ric[i,j])
    )
)$

cdisplay(ric)$

leinstein(false)$

for i:1 while i<=4 do (
    for j:1 while j<=4 do (
        lein[i,j] : ratsimp(lein[i,j])
    )
)$

cdisplay(lein)$

diff(f,x);

diff(f*v_s,z,2);

genmatrix(lein, 4, 4)$
subst([g=gdet],%)$
subst([u_from_v_s_f],%)$
ev(%,diff)$
ratsimp(%)$
subst([diff(f,r_s)=df],%)$
subst([diff(f,r_s,2)=d2f],%)$
/* "moving with bubble" means first derivatives are zero ... but not 2nd derivatives, right */ 
subst([
    diff(diff(u,t),x)=0,
    diff(diff(u,t),y)=0,
    diff(diff(u,t),z)=0,
    diff(diff(u,x),y)=0,
    diff(diff(u,x),z)=0,
    diff(diff(u,y),z)=0,
    diff(u,t,2)=0,
    diff(u,x,2)=0,
    diff(u,y,2)=0,
    diff(u,z,2)=0
], %)$
ratsimp(%)$
/*
r_s_,t = -v_s*(x - x_s)/r_s
r_s_,x = (x - x_s)/r_s
r_s_,y = y/r_s
r_s_,z = z/r_s
*/
subst([
    diff(r_s,t,2) = diff( -v_s * (x - x_s) / r_s, t),
    diff(r_s,x,2) = diff( (x - x_s) / r_s, x),
    diff(r_s,y,2) = diff( y / r_s, y),
    diff(r_s,z,2) = diff( z / r_s, z)
],%)$
subst([
    diff(r_s,t) = -v_s * (x - x_s) / r_s,
    diff(r_s,x) = (x - x_s) / r_s,
    diff(r_s,y) = y / r_s,
    diff(r_s,z) = z / r_s
],%)$
ratsimp(%)$
lein_mat : %
/*
subst([
    diff(u,t)=0,
    diff(u,x)=0,
    diff(u,y)=0,
    diff(u,z)=0,
    diff(diff(u,t),x)=0,
    diff(diff(u,t),y)=0,
    diff(diff(u,t),z)=0,
    diff(diff(u,x),y)=0,
    diff(diff(u,x),z)=0,
    diff(diff(u,y),z)=0,
    diff(u,t,2)=0,
    diff(u,x,2)=0,
    diff(u,y,2)=0,
    diff(u,z,2)=0
], %);
*/;

lein_mat[1,1];