1 files changed, 165 insertions, 0 deletions
diff --git a/01-bondi.mac b/01-bondi.mac
new file mode 100644
index 0000000..dda6bf8
--- /dev/null
+++ b/01-bondi.mac
@@ -0,0 +1,165 @@
+load("setup.mac")$
+load("tensors.mac")$
+load("series.mac")$
+
+
+/* Equations */
+
+/* Stress-energy tensor of the theory */
+T_ab : zeromatrix(4,4)$
+
+/* Einstein equations tensor */
+for i : 1 thru 4 do for j : 1 thru 4 do (
+  E_ab[i,j] : expand(ric[i,j] - ric_s*lg[i,j]/2 - 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])))$
+
+/* solve(uE_ab[1,2]=0, diff(β,r))[1], */
+
+/* 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);
+
+
+/* Series ansatz */
+
+expansion_order : 6;
+
+/* Sequence elements inherit dependance on u and θ */
+depends([c, H, M, N, γc, βc, Ac, Bc], [u, θ])$
+
+/* Simplify ahead to speed up calculation */
+/* γβ_facts : [H=0, γc[2]=0]$ */
+
+/* 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)])$
+
+main_ansatz : build_main_ansatz(expansion_order)$
+
+funcs : [A, B, β, γ, β_expr, γ_expr];
+
+
+/**/
+
+/* Simplify ahead to speed up calculation */
+γβ_facts : [H=0, γc[2]=0];
+
+coeffs : γβ_facts;
+
+/* Equation 1 : determine β */
+e_hs_eqs[1] : expand_equation(solve(hs_eqs[1], diff(β, r))[1], coeffs);
+ce_hs_eqs[1] : coeff_equations(e_hs_eqs[1]);
+β_coeffs : block([len, lst],
+  len : length(ce_hs_eqs[1]),
+  lst : makelist(solve(ce_hs_eqs[1][i][2], βc[i]), i, 1, len),
+  apply(append, lst));
+
+coeffs : append(coeffs, β_coeffs);
+
+/* Equation 2 : determine A */
+e_hs_eqs[2] : expand_equation(expand(exp(2*β)*hs_eqs[2]), coeffs);
+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 */
+e_hs_eqs[3] : expand_equation(hs_eqs[3], coeffs);
+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)$
+
+/* Equation 4 : determine γ (evolution equation) */
+
+/*
+We wish M,N,C to be related to monopole, dipole and quadrupole moments
+respectively, also when converted from the static case. Thus, Ac[3] as to be of
+following form.
+*/
+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);
+
+e_ev_eqs[2] : expand_equation(ev_eqs[2], coeffs);
+ce_ev_eqs[2] : coeff_equations(e_ev_eqs[2]);
+
+eq_quad : expand(4*solve(ce_ev_eqs[2][1][2], diff(C,u))[1]);
+
+/* Convoluted way to show that the results for C_u are consistent */
+e_ev_eqs[1] : expand_equation(ev_eqs[1], coeffs);
+ce_ev_eqs[1] : coeff_equations(e_ev_eqs[1]);
+
+eq_quad_alt : trigreduce(expand(trigsimp(4*solve(ce_ev_eqs[2][2][2], diff(C,u))[1])));
+gammac4 : trigreduce(expand(trigsimp(trigreduce(expand(solve(ce_ev_eqs[1][2][2], diff(γc[4],u))))[1])));
+eq_quad_alt_res : expand(subst(gammac4, eq_quad_alt))
+
+
+/* Supplementary conditions */
+
+/* Supplementary Condition 1 */
+e_sc_eqs[1] : expand_equation(sc_eqs[1], coeffs);
+eq_mono : block(
+  [eq : coeff(e_sc_eqs[1], r, -2)=0],
+  trigreduce(expand(solve(eq, diff(M,u)))));
+
+/* Supplementary Condition 2 */
+e_sc_eqs[2] : expand_equation(sc_eqs[2], coeffs);
+eq_dipo : block(
+  [eq : coeff(e_sc_eqs[2], r, -2)=0],
+  trigreduce(expand(-3*solve(eq, diff(N,u)))));
+
+
+
+/* BUG: should not generate terms with r^-3 and so on ... */
+block([expansion_order : 2],
+  expand_equation(sc_eqs[1], coeffs))
+
+block([expansion_order : 2],
+  expand_equation(sc_eqs[2], coeffs))