/* 1920x1080 - HD screen */
linel: 230;
/* 3840x1600 - Ultrawide screen */
linel: 470;
/* 1366x768 - Poor man's thinkpad */
linel: 150;

/* 3840x1600 - Ultrawide screen with work on side */
linel: 350;

derivabbrev: true;

load("ctensor");

dim: 4;
cframe_flag: false;
ct_coords: [u, r, θ, ϕ];

depends([U, V, β, γ], [u, r, θ]);

lg : zeromatrix(4, 4);
/* du du */
lg[1, 1] : V * r^-1 * exp(2 * β) - U^2 * r^2 * exp(2 * γ);

/* du dr */
lg[1, 2] : exp(2 * β);
lg[2, 1] : exp(2 * β);

/* du dθ */
lg[1, 3] : U * r^2 * exp(2 * γ);
lg[3, 1] : U * r^2 * exp(2 * γ);

/* dθ dθ */
lg[3, 3] : - r^2 * exp(2 * γ);

/* dϕ dϕ */
lg[4, 4] : - r^2 * exp(- 2 * γ) * sin(θ)^2;

/* Load metric and compute inverse */
cmetric();

/* Force simplification of some expressions */
ug[2,2] : expand(ug[2,2]),simp;

/* mcs[i, j, k]: Christoffel symbols of the second kind (one upper index - k) */
christof(false);
/* ric[i, j]: R_ab Ricci tensor*/
ricci(false);
/* uric[i, j]: R^a_b Ricci tensor */
uricci(false);

/* Bondi - Eq. 17 */
sum(sum(ug[α, ε] * mcs[α, ε, 1], ε, 1, 4), α, 1, 4),expand;

/** General utilities **/
exp_taylor(x) := subst(t=x, taylor(exp(t), t, 0, 3));

/** Main equations **/

/* Main equations, Bondi - Eqs. 22 - 25 */
eq1 : ric[2, 2] = 0, expand;
eq2 : - 2 * r^2 * ric[2, 3] = 0, expand;
eq3 : ric[3, 3] * exp(2 * (β - γ)) - r^2 * uric[4, 4] * exp(2 * β) = 0, expand;
eq4 : - uric[4, 4] * exp(2 * β) * r^2 = 0, expand;

depends([H], [u, θ]);
depends([γ1, γ2, γ3, γ4], [u, θ]);
depends([β1, β2, β3, β4], [u, θ]);
depends([U1, U2, U3, U4], [u, θ]);
depends([VV1, V0, V1, V2, V3, V4], [u, θ]);

ansatz: [γ = γ1*r^-1 + γ2*r^-2 + γ3*r^-3 + γ4*r^-4,
         β = H + β1*r^-1 + β2*r^-2 + β3*r^-3 + β4*r^-4,
         U = U1*r^-1 + U2*r^-2 + U3*r^-3 + U4*r^-4,
         V = VV1*r + V0 + V1*r^-1 + V2*r^-2 + V3*r^-3 + V4*r^-4];

/* Derive β from eq1; Replace in eq1 -> β2 = - γ1^2 / 4, β1 = 0*/
tmp: combine(ev(subst(ansatz, eq1),derivlist(r),expand));
/* Solve for β_1 */
solve(coeff(tmp, 1/r, 3), β1);
/* Solve for β_2 */
solve(coeff(tmp, 1/r, 4), β2);

matchdeclare (nn, lambda ([e], not is(equal(listofvars(e), [H]))));
defrule (r1, %e^nn, (partition(nn, H), exp_taylor(first(%%)) * exp(second(%%))));

/* Derive U from eq2 */
tmp: combine(ev(apply1(subst(ansatz, eq2), r1),derivlist(r,θ),expand))$
/* Solve for U_1 */
solve(coeff(tmp, r, 1), U1);


/* Derive V from eq2 */
append(ansatz, [U1=2*diff(H,θ)*exp(2*H)]);
tmp: combine(ev(apply1(subst(%, eq3), r1),derivlist(r,θ),expand))$
/* Solve for V_{-1} */
solve(coeff(tmp, r, 0), VV1),expand;


/* subst(ansatz, eq4),derivlist(r,θ),expand; */


/** Supplementary conditions **/
depends([c, C, N, M], [u, θ]);

ansatz: [γ = c*r^-1 + (C - c^3/6)*r^-3,
         β = -(1/4)*c^2*r^-2,
         U = -(diff(c,θ) + 2*c*cot(θ))*r^-2
             + (2*N + 3*c*diff(c,θ) + 4*c^2*cot(θ))*r^-3
             + (1/2)*(3*diff(C,θ) + 6*C*cot(θ) - 6*c*N - 8*c^2*diff(c,θ)
               - 8*c^3*cot(θ))*r^-4,
         V = r - 2*M
             - (diff(N,θ) + N*cot(θ) - diff(c,θ)^2 - 4*c*diff(c,θ)*cot(θ)
               - c^2*(1 + 8*cot(θ)^2)/2)*r^-1
             - (1/2)*(diff(C,θ,2) + 3*diff(C,θ)*cot(θ) - 2*C
               + 6*N*(diff(c,θ) + 2*c*cot(θ)) + 8*c*(diff(c,θ)^2 + 3*c*diff(c,θ)
               + 2*c^2*cot(θ)^2))*r^-2];

C_u: (2*c^2 + 2*c*M + N*cot(θ) - diff(N,θ)) / 4;

matchdeclare (mm, lambda ([e], true));
defrule (r2, %e^mm, exp_taylor(mm));

/* Mass loss formula, Bondi - Eq. 35 */
tmp: combine(ev(subst(ansatz, ric[1,1]),derivlist(u,r,θ),expand))$
distrib(trigsimp(first(solve(coeff(%, 1/r, 2)=0, diff(M,u))))),expand$
sup_cond1: substpart(trigreduce(part(%, 2, 1)), %, 2, 1);

/* Bondi - Eq. 36  */
tmp: combine(ev(subst(ansatz, ric[1,3]),derivlist(u,r,θ),expand))$
combine(expand(apply1(tmp, r2)))$
sup_cond2: -solve(coeff(%, 1/r, 2)=0, 3*diff(N,u));