##---------------------------------------
# In this Maple sheet, we demonstrate the linear dependence of the components in vector p (Eq. 21) on the particle velocity v, see Section 4.5.
##---------------------------------------

# Restart maple
restart;

# Include dependent libraries
with(plots):
with(LinearAlgebra):

# Initialize the random number generator, with a free rolls.
roll := rand(1..20):
roll();
roll();
roll();
roll();
roll();
roll();

##---------------------------------------
# Rather than demonstrating the linear dependence symbolically, we show it for a random vector field.
# The symbolic expressions quickly become very complicated, since a linear 18x18 system must be solved along the way.
##---------------------------------------

# Create the u component and all its first and second-order derivatives randomly.
uc := (roll()-10)/10:
u_x := (roll()-10)/10:
u_y := (roll()-10)/10:
u_z := (roll()-10)/10:
u_t := (roll()-10)/10:
u_xx := (roll()-10)/10:
u_yy := (roll()-10)/10:
u_zz := (roll()-10)/10:
u_tt := (roll()-10)/10:
u_xy := (roll()-10)/10:
u_xz := (roll()-10)/10:
u_xt := (roll()-10)/10:
u_yz := (roll()-10)/10:
u_yt := (roll()-10)/10:
u_zt := (roll()-10)/10:

# Create the v component and all its first and second-order derivatives randomly.
vc := (roll()-10)/10:
v_x := (roll()-10)/10:
v_y := (roll()-10)/10:
v_z := (roll()-10)/10:
v_t := (roll()-10)/10:
v_xx := (roll()-10)/10:
v_yy := (roll()-10)/10:
v_zz := (roll()-10)/10:
v_tt := (roll()-10)/10:
v_xy := (roll()-10)/10:
v_xz := (roll()-10)/10:
v_xt := (roll()-10)/10:
v_yz := (roll()-10)/10:
v_yt := (roll()-10)/10:
v_zt := (roll()-10)/10:

# Create the w component and all its first and second-order derivatives randomly.
wc := (roll()-10)/10:
w_x := (roll()-10)/10:
w_y := (roll()-10)/10:
w_z := (roll()-10)/10:
w_t := (roll()-10)/10:
w_xx := (roll()-10)/10:
w_yy := (roll()-10)/10:
w_zz := (roll()-10)/10:
w_tt := (roll()-10)/10:
w_xy := (roll()-10)/10:
w_xz := (roll()-10)/10:
w_xt := (roll()-10)/10:
w_yz := (roll()-10)/10:
w_yt := (roll()-10)/10:
w_zt := (roll()-10)/10:

# We create the test vector field in Taylor expansion around the space-time point (0,0,0,0).
t0 := 0;
x0 := 0;
y0 := 0;
z0 := 0;

# We form the Taylor expansion of our vector field (uf, vf, wf), using the derivatives that we just defined above.
uf := uc + u_x*(x-x0) + u_y*(y-y0) + u_z*(z-z0) + u_t*(t-t0)
+ (1/2)*u_xx*(x-x0)*(x-x0)
+ (1/2)*u_yy*(y-y0)*(y-y0)
+ (1/2)*u_zz*(z-z0)*(z-z0)
+ (1/2)*u_tt*(t-t0)*(t-t0)
+ (1/1)*u_xy*(x-x0)*(y-y0)
+ (1/1)*u_xz*(x-x0)*(z-z0)
+ (1/1)*u_xt*(x-x0)*(t-t0)
+ (1/1)*u_yz*(y-y0)*(z-z0)
+ (1/1)*u_yt*(y-y0)*(t-t0)
+ (1/1)*u_zt*(z-z0)*(t-t0);
vf := vc + v_x*(x-x0) + v_y*(y-y0) + v_z*(z-z0) + v_t*(t-t0)
+ (1/2)*v_xx*(x-x0)*(x-x0)
+ (1/2)*v_yy*(y-y0)*(y-y0)
+ (1/2)*v_zz*(z-z0)*(z-z0)
+ (1/2)*v_tt*(t-t0)*(t-t0)
+ (1/1)*v_xy*(x-x0)*(y-y0)
+ (1/1)*v_xz*(x-x0)*(z-z0)
+ (1/1)*v_xt*(x-x0)*(t-t0)
+ (1/1)*v_yz*(y-y0)*(z-z0)
+ (1/1)*v_yt*(y-y0)*(t-t0)
+ (1/1)*v_zt*(z-z0)*(t-t0);
wf := wc + w_x*(x-x0) + w_y*(y-y0) + w_z*(z-z0) + w_t*(t-t0)
+ (1/2)*w_xx*(x-x0)*(x-x0)
+ (1/2)*w_yy*(y-y0)*(y-y0)
+ (1/2)*w_zz*(z-z0)*(z-z0)
+ (1/2)*w_tt*(t-t0)*(t-t0)
+ (1/1)*w_xy*(x-x0)*(y-y0)
+ (1/1)*w_xz*(x-x0)*(z-z0)
+ (1/1)*w_xt*(x-x0)*(t-t0)
+ (1/1)*w_yz*(y-y0)*(z-z0)
+ (1/1)*w_yt*(y-y0)*(t-t0)
+ (1/1)*w_zt*(z-z0)*(t-t0);

# Gravity (ug,vg,wg) and particle response time r are also chosen randomly.
ug := (roll()-10)/10;
vg := (roll()-10)/10;
wg := (roll()-10)/10;
r := roll()/20;


##---------------------------------------
# Next, we compute the Jacobian and the temporal derivative of our vector field.
##---------------------------------------

# First, we put the three vector field components (uf,vf,wf) into a vector Uf.
Uf := Vector([uf,vf,wf]);

# The Jacobian is computed by arranging the spatial derivatives as column vectors in a matrix.
Jf := 
Matrix([
[diff(Uf[1],x), diff(Uf[1],y), diff(Uf[1],z)],
[diff(Uf[2],x), diff(Uf[2],y), diff(Uf[2],z)],
[diff(Uf[3],x), diff(Uf[3],y), diff(Uf[3],z)]]);

# Finally, the temporal derivative is computed.
Uf_t := Vector([diff(Uf[1],t),diff(Uf[2],t),diff(Uf[3],t)]);


##---------------------------------------
# Given the vector field and its derivatives, we can not setup the matrix M, see Eq. (22).
##---------------------------------------

# First, we introduce a few abbreviations for the location X, the particle velocity V and the gravity vector G.
X := Vector([x,y,z]);
V := Vector([u,v,w]);
G := Vector([ug,vg,wg]);

# Next, we apply the operator _p of Eq. (1) to location X, flow velocity U, inertial particle velocity V and gravity G.
# Define x_p
XM := Matrix([
[  0   ,-X[3],  X[2]],
[ X[3],  0   , -X[1]],
[-X[2], X[1],   0   ]]);

# Define u_p
UM := Matrix([
[  0   ,-Uf[3],  Uf[2]],
[ Uf[3],  0   , -Uf[1]],
[-Uf[2], Uf[1],   0   ]]);

# Define g_p
GM := Matrix([
[  0   ,-G[3],  G[2]],
[ G[3],  0   , -G[1]],
[-G[2], G[1],   0   ]]);

# Define v_p
VM := Matrix([
[  0   ,-V[3],  V[2]],
[ V[3],  0   , -V[1]],
[-V[2], V[1],   0   ]]);


# Given all symbols above, we can now define matrix M according to Eq. (22).
# Using a small helper variable:
H1 := Multiply(-Jf,XM) + UM;

M := (1/r)*Matrix([
[H1[1,1],H1[1,2],H1[1,3], Jf[1,1],Jf[1,2],Jf[1,3], XM[1,1],XM[1,2],XM[1,3], -1,0,0,  0,0,0, 0,0,0 ],
[H1[2,1],H1[2,2],H1[2,3], Jf[2,1],Jf[2,2],Jf[2,3], XM[2,1],XM[2,2],XM[2,3], 0,-1,0,  0,0,0, 0,0,0 ],
[H1[3,1],H1[3,2],H1[3,3], Jf[3,1],Jf[3,2],Jf[3,3], XM[3,1],XM[3,2],XM[3,3], 0,0,-1,  0,0,0, 0,0,0 ]])
+
Matrix([
[GM[1,1],GM[1,2],GM[1,3], 0,0,0, 2*VM[1,1],2*VM[1,2],2*VM[1,3], 0,0,0,  XM[1,1],XM[1,2],XM[1,3],  -1,0,0 ],
[GM[2,1],GM[2,2],GM[2,3], 0,0,0, 2*VM[2,1],2*VM[2,2],2*VM[2,3], 0,0,0,  XM[2,1],XM[2,2],XM[2,3],  0,-1,0 ],
[GM[3,1],GM[3,2],GM[3,3], 0,0,0, 2*VM[3,1],2*VM[3,2],2*VM[3,3], 0,0,0,  XM[3,1],XM[3,2],XM[3,3],  0,0,-1 ]
]);


##---------------------------------------
# Now that we have an analytic expression for matrix M, we compute the system matrix \hat{M} in Eq. (25) and the right hand side Eq. (26).
##---------------------------------------

# We sample the neighborhood U in the spatial domain n times
n := 50;

# The random 6d sample points all have the same particle velocity V, since we sample in space only.
for i from 1 by 1 to n do 
xx0[i] := (roll()-10)/10:
yy0[i] := (roll()-10)/10:
zz0[i] := (roll()-10)/10:
uu0[i] := u0:
vv0[i] := v0:
ww0[i] := w0:
end do:

# Evaluate M and the temporal derivative U_t of the vector at the sample points
for i from 1 by 1 to n do
M0[i] := eval(M,{x=xx0[i],y=yy0[i],z=zz0[i],t=t0, u=uu0[i], v=vv0[i],w=ww0[i]  }):
U0_t[i] := eval(Uf_t,{x=xx0[i],y=yy0[i],z=zz0[i],t=t0, u=uu0[i], v=vv0[i],w=ww0[i]   }):
end do:

# Compute the system matrix  M^T M and the right hand side M^T U_t by summing up all the samples (Eqs. 25 and 26), which gives an 18x18 system matrix and a 18x1 right hand side vector.
MM := Matrix(18,18);
YY := Vector(18);
for i from 1 by 1 to n do 
MM := MM + Multiply(Transpose(M0[i]),M0[i]):
YY := YY + (1/r)*Multiply(Transpose(M0[i]),U0_t[i]):
end do:

# Compute the rank of the system matrix, which gives a rank of 12.
Rank(MM);

# Following Section 4.4, we introduce a regularization matrix E.
EE := Matrix(18,18);
EE[13,13] := 1:
EE[14,14] := 1:
EE[15,15] := 1:
EE[16,16] := 1:
EE[17,17] := 1:
EE[18,18] := 1:

# Apply the regularization as in Eq. (30), which gives a full rank of 18.
MM := MM + EE:
Rank(MM);

##---------------------------------------
# Now comes the interesting part. How does the solution depend on particle velocity v?
##---------------------------------------

# In our right hand side of the linear system only the components 7-9 depend on v.
YY[1];
YY[2];
YY[3];
YY[4];
YY[5];
YY[6];
YY[7];
YY[8];
YY[9];
YY[10];
YY[11];
YY[12];
YY[13];
YY[14];
YY[15];
YY[16];
YY[17];
YY[18];

# Now, we solve the linear system in Eq. (24), which produces the optimal reference frame parameters in vector p:
P := Multiply(MatrixInverse(MM),YY);

# In vector p, only components 10-12 depend on the particle velocity.
simplify(P[1]);
simplify(P[2]);
simplify(P[3]);
simplify(P[4]);
simplify(P[5]);
simplify(P[6]);
simplify(P[7]);
simplify(P[8]);
simplify(P[9]);
simplify(P[10]);
simplify(P[11]);
simplify(P[12]);
simplify(P[13]);
simplify(P[14]);
simplify(P[15]);
simplify(P[16]);
simplify(P[17]);
simplify(P[18]);

# Looking at Eq. (21), we see that these components correspond to \ddot{d}.
# In other words, only \ddot{d} depends on the particle velocity.
# Further, we can see that \dddot{d} and \dddot{S} both become zero as expected from the regularization.