/* 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 * β);

/* dr 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: taylor(exp(x), x, 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;
β_ansatz: H + β1*r^-1 + β2*r^-2 + β3*r^-3 + β4*r^-4;
U_ansatz: U1*r^-1 + U2*r^-2 + U3*r^-3 + U4*r^-4;
V_ansatz: VV1*r + V0 + V1*r^-1 + V2*r^-2 + V3*r^-3 + V4*r^-4;
ansatz: [γ=γ_ansatz, β=β_ansatz, U=U_ansatz, V=V_ansatz];


/* 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);

/* Derive U from eq2 */
sub: cons(exp(2*γ - 2*β)=exp(-2*H)*subst(x=2*(γ - β + H), exp_taylor), ansatz)$
tmp: combine(ev(subst(sub,eq2),derivlist(r,θ),expand))$
/* Solve for U_1 */
solve(coeff(tmp, r, 1), U1);




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

combine(ev(apply1(subst(ansatz, eq2), r1),derivlist(r,θ),expand));
solve(coeff(%, r, 1), U1);


/* Derive V from eq2 */
tmp: subst(ansatz, eq3),derivlist(r,θ),expand$
ex: part(coeff(tmp, r, 0), 1, 6, 2, 2, allbut(9))$
sub: [exp(ex + 2*H)=exp(2*H)*subst([x=ex], taylor(exp(x), x, 0, 3)),
      exp(- ex - 2*H)=exp(-2*H)*subst([x=-ex], taylor(exp(x), x, 0, 3))]$
tmp2: combine(expand(subst(sub, tmp)))$
/* Solve for V_{-1} */
solve(trigsimp(ev(subst([U1=2*diff(H,θ)*exp(2*H)], coeff(tmp2, r, 0)),derivlist(θ))), VV1);


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

/*
U_ansatz: 2 * diff(H, θ) * exp(2 * H) / r
V_ansatz: r * exp(2 * H) * (1 + 2 * diff(H, θ) * cot(θ) + 4 * diff(H, θ) ^ 2 + 2 * diff(H, θ, 2))
*/


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

γ_ansatz: c*r^-1 + (C - c^3/6)*r^-3;
/* Assumption given H=0, not on paper*/
β_ansatz: -(1/4)*c^2*r^-2;
U_ansatz: -(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_ansatz: 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;

ansatz: [γ=γ_ansatz, β=β_ansatz, U=U_ansatz, V=V_ansatz];

/* 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))$
sup_cond2: distrib(trigsimp(-3*first(solve(coeff(%, 1/r, 2)=0, diff(N,u))))),expand;

taylor(rhs(sup_cond2), x, 0, 3)

rhs(sup_cond2)

ex: exp(- 2*c*r^-1 - c^2*r^-2/2 + c^3*r^-3/3 - 2*C*r^-3)

ex: exp(2*c*r^-1 + c^2*r^-2/2 - c^3*r^-3/3 + 2*C*r^-3)

8*c*diff(c,u)*cot(θ) - (3/2)*diff(c, u)*diff(c,θ) * ex + c*diff(diff(c,u),θ)*ex/2 + diff(M, θ)*ex + 5*diff(c,u)*diff(c,θ)/2 + 9*c*diff(diff(c,u),θ)/2;

ser: subst([x=2*c*r^-1 + c^2*r^-2/2 - c^3*r^-3/3 + 2*C*r^-3], taylor(exp(x), x, 0, 3));
tmp2: subst([ex=ser], tmp),expand;

tmp3: combine(tmp2);

coeff(tmp3, 1/r, 2)=0

distrib(trigsimp(-3*first(solve(coeff(%, 1/r, 2)=0, diff(N,u))))),expand;