load("setup.mac")$
load("tensors.mac")$
load("series.mac")$


/* Equations */

depends(Φ, [u,r,θ,ϕ]);

/* Einstein tensor */
for i : 1 thru 4 do for j : 1 thru 4 do (
  G_ab[i,j] : expand(ric[i,j] - ric_s*lg[i,j]/2));

/* Equation (8) */

/*
(1-Φ^2) * ric[i,j]
= 4*cdiff(Φ,i)*cdiff(Φ,j)
- lg[i,j]*contr(ug[a,b]*cdiff(Φ,a)*cdiff(Φ,b))
- 2*Φ*(cdiff(cdiff(Φ, j), i) - contr(mcs[i,j,k]*cdiff(Φ,k),[k]))

Add -(1/2)*lg[i,j]*ric_s = 0 on both sides.
Yields G_ab on lhs and does nothing on rhs
*/

/* Stress-energy tensor of the theory */
for i : 1 thru 4 do for j : 1 thru 4 do (
  T_ab[i,j] : expand(
    (1/8/%pi)*(
      4*cdiff(Φ,i)*cdiff(Φ,j)
      - lg[i,j]*contr(ug[a,b]*cdiff(Φ,a)*cdiff(Φ,b),[a,b])
      - 2*Φ*(cdiff(cdiff(Φ, j), i) - contr(mcs[i,j,k]*cdiff(Φ,k),[k])))))$

/* Einstein equations tensor */
for i : 1 thru 4 do for j : 1 thru 4 do (
  E_ab[i,j] : expand(G_ab[i,j] - 8*%pi*T_ab[i,j]))$

/* Einstein equations tensor with mixed indices */
for i : 1 thru 4 do for j : 1 thru 4 do (
  uE_ab[i,j] : expand(contr(ug[i,k]*E_ab[k,j], [k])))$

/* Hypersurface equations */
hs_eqs : [
/* Determine β */
uE_ab[1,2]=0,
/* Determine A */
uE_ab[1,3]=0,
/* uE_ab[1,4]=0 is trivial, would determine A' otherwise */
/* Determine B */
uE_ab[1,1]=0
]$

/* Evolution equations */
ev_eqs : [
/* Determine γ */
E_ab[3,3]=0,
/* Should be fulfilled automatically (?) */
E_ab[4,4]=0
/* Fulfilled trivially */
/* E_ab[3,4]=0 */
/* E_ab[4,3]=0 */
]$

/* Suplementary conditions */
sc_eqs : [
uE_ab[2,1]=0,
uE_ab[2,3]=0
/* Fulfilled trivially */
/* uE_ab[2,4]=0 */
]$

hs_eqs : trigreduce(hs_eqs)$
sc_eqs: trigreduce(sc_eqs)$

hs_eqs : subst([diff(Φ,ϕ)=0], hs_eqs)$
ev_eqs : subst([diff(Φ,ϕ)=0], ev_eqs)$
sc_eqs : subst([diff(Φ,ϕ)=0], sc_eqs)$

/* Scalar wave equation */
sw_eq : trigreduce(expand(box_op(Φ)=0))$


/* Series ansatz */

expansion_order : 6;

/* Sequence elements inherit dependance on u and θ */
depends([c, H, M, N, γc, βc, Ac, Bc, Φc], [u, θ])$

/* Main series ansatz */
build_main_ansatz(n) := block([], [
  γ = c/r + neg_pow_series(γc, r, 2, n),
  β = H + neg_pow_series(βc, r, 1, n),
  A = neg_pow_series(Ac, r, 1, n),
  B = neg_pow_series(Bc, r, -1, n),
  Φ = neg_pow_series(Φc, r, 0, n)])$

main_ansatz : build_main_ansatz(expansion_order)$

funcs : [A, B, β, γ, β_expr, γ_expr, Φ];

/**/

/* Simplify ahead to speed up calculation */
γβ_facts : [H=0, γc[2]=0];
sim_facts : ev([
  /* From Equation 4 (1/r term) */
  diff(Φc[0],u)=0,
  diff(Φc[0],u,2)=0,
  /* From Equation 4 (1/r^2 term) */
  diff(Φc[0],θ)=0,
  diff(Φc[0],θ, 2)=0,
  diff(diff(Φc[0],θ),u)=0,
  diff(diff(Φc[0],u),θ)=0,
  diff(diff(Φc[0],u),θ,2)=0],
  eval,diff);

coeffs : γβ_facts;
coeffs : append(γβ_facts, sim_facts);

/* Equation 1 : determine β */
e_hs_eqs[1] : expand_equation(expand(exp(2*β)*hs_eqs[1]), coeffs)$
ce_hs_eqs[1] : coeff_equations(e_hs_eqs[1])$
β_coeffs : block(
  [vals : []],
  for i : 1 thru 5 do (
    vals : append(vals, solve(subst(vals, ce_hs_eqs[1][i][2]), βc[i]))),
  vals)$

coeffs : append(coeffs, β_coeffs)$

/* Equation 2 : determine A */
e_hs_eqs[2] : expand_equation(expand(exp(2*β)*hs_eqs[2]), coeffs)$
e_hs_eqs[2] : subst(coeffs, e_hs_eqs[2])$
ce_hs_eqs[2] : coeff_equations(e_hs_eqs[2])$

A_coeffs : block(
  [vals : []],
  for i : 1 thru expansion_order do (
    if not(i=3) then (
      vals : append(vals, solve(subst(vals, ce_hs_eqs[2][i][2]), Ac[i])))),
  ev(vals, eval, diff, expand))$

coeffs : append(coeffs, A_coeffs)$

/* Equation 3 : determine B */
print_progress : true;
e_hs_eqs[3] : expand_equation(expand(exp(2*β)*hs_eqs[3]), coeffs)$
/* TODO: 17 31 34 slow due to squares? */
/* BUG: diff(λc[0],θ)=0 does not get substituted; due to missing ev? */
e_hs_eqs[3] : subst(coeffs, e_hs_eqs[3])$
ce_hs_eqs[3] : coeff_equations(e_hs_eqs[3])$
B_coeffs : block(
  [vals : []],
  for i : 1 thru expansion_order do (
    if i=2 then (
      vals : append(vals, [Bc[0]=-2*M]))
    else (
      vals : append(vals, solve(subst(vals, ce_hs_eqs[3][i][2]), Bc[i-2])))),
  ev(vals, eval, diff, expand))$

coeffs : append(coeffs, B_coeffs)$


/* Evolution equations */

depends(D, [u, θ]);

/* Equation 4 : determine λ */
e_sw_eq : expand_equation(sw_eq, coeffs)$
ce_sw_eq : coeff_equations(e_sw_eq)$
Φu_coeffs : block(
  [vals : [Φc[1]=D]],
  for i : 3 thru 6 do (
    vals : append(vals, solve(subst(vals, ce_sw_eq[i][2]), diff(Φc[i-1],u)))),
  vals);

coeffs : append(coeffs, Φu_coeffs)$

depends(C, [u, θ]);
extra_coeffs : [γc[3]=C-c^3/6, Ac[3]=2*N + 3*c*diff(c,θ) + 4*c^2*cot(θ)]$
coeffs : append(coeffs, extra_coeffs)$

/* Equation 5 : determine γ */
e_ev_eqs[2] : expand_equation(trigreduce(expand(exp(2*β)*exp(2*γ)*ev_eqs[2]/sin(θ)^2)), coeffs)$
e_ev_eqs[2] : subst(coeffs, e_ev_eqs[2])$
ce_ev_eqs[2] : coeff_equations(e_ev_eqs[2])$

eq_quad : trigreduce(expand(trigsimp(subst(coeffs,
        4*solve(ce_ev_eqs[2][2][2], diff(C, u))[1]))))$

subst(coeffs,ce_ev_eqs[2][1]);

/* Relation of c to c_AB */
ce_ev_eqs[2][1]; /* TODO: D_u * c = c_u * D ??? */


/* Supplementary conditions */

/* Supplementary Condition 1 */
e_sc_eqs[1] : expand_equation(expand(exp(2*β)*sc_eqs[1]), coeffs)$
/* BUG */
e_sc_eqs[1] : subst(sim_facts, e_sc_eqs[1])$
eq_mono : block(
  [eq : coeff(e_sc_eqs[1], r, -2)=0, sim_facts:sim_facts],
  push(subst(sim_facts, expand(diff(Φu_coeffs[2],u))), sim_facts),
  eq : solve(eq, diff(M,u))[1],
  eq : subst(sim_facts, eq),
  trigreduce(expand(trigsimp(eq))));

/* Supplementary Condition 2 */
e_sc_eqs[2] : block(
  /* BUG: Single A with no further factors crashes */
  [tmp : lhs(expand(exp(2*β)*sc_eqs[2]+A))=0],
  tmp : expand_equation(tmp, coeffs),
  expand(tmp - expand(subst(append(main_ansatz, coeffs), A))));

eq_dipo : block(
  [eq : coeff(e_sc_eqs[2], r, -2)=0],
  subst(sim_facts, trigreduce(expand(-3*solve(eq, diff(N,u))))));


/* Cache */

path(name) := printf(false, "cache/04-bbmd-jor/~a.lisp", name);

save(path("e_hs_eqs"), e_hs_eqs);
save(path("ce_hs_eqs"), ce_hs_eqs);
save(path("e_sw_eq"), e_sw_eq);
save(path("ce_sw_eq"), ce_sw_eq);
save(path("e_ev_eqs"), e_ev_eqs);
save(path("ce_ev_eqs"), ce_ev_eqs);

save(path("beta_coeffs"), β_coeffs);
save(path("a_coeffs"), A_coeffs);
save(path("b_coeffs"), B_coeffs);
save(path("phiu_coeffs"), Φu_coeffs);

save(path("e_sc_eqs"), e_sc_eqs);
save(path("coeffs"), coeffs)

save(path("eq_quad"), eq_quad);
save(path("eq_mono"), eq_mono);
save(path("eq_dipo"), eq_dipo);

load_cache() := block(
  [cache : [
  "e_hs_eqs", "ce_hs_eqs",
  "e_sw_eq", "ce_sw_eq",
  "e_ev_eqs", "ce_ev_eqs",
  "beta_coeffs", "a_coeffs", "b_coeffs", "phiu_coeffs",
  "e_sc_eqs",
  "coeffs",
  "eq_quad", "eq_mono", "eq_dipo"]],
  for file in cache do load(path(file)));

load_cache();