Refactor Bondi

1 files changed, 41 insertions, 43 deletions

diff --git a/01-bondi.mac b/01-bondi.mac
index 9240820..e9133ab 100644
--- a/01-bondi.mac
+++ b/01-bondi.mac
@@ -14,11 +14,11 @@ for i : 1 thru 4 do for j : 1 thru 4 do (
 
 /* 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]))$
+  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])))$
+  uE_ab[i,j] : expand(contr(ug[i,k]*E_ab[k,j], [k])));
 
 /* solve(uE_ab[1,2]=0, diff(β,r))[1], */
 
@@ -31,7 +31,7 @@ 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 : [
@@ -42,7 +42,7 @@ E_ab[4,4]=0
 /* Fulfilled trivially */
 /* E_ab[3,4]=0 */
 /* E_ab[4,3]=0 */
-];
+]$
 
 /* Suplementary conditions */
 sc_eqs : [
@@ -50,9 +50,9 @@ uE_ab[2,1]=0,
 uE_ab[2,3]=0
 /* Fulfilled trivially */
 /* uE_ab[2,4]=0 */
-];
+]$
 
-hs_eqs : trigreduce(hs_eqs);
+hs_eqs : trigreduce(hs_eqs)$
 
 
 /* Series ansatz */
@@ -62,9 +62,6 @@ 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),
@@ -85,38 +82,41 @@ funcs : [A, B, β, γ, β_expr, γ_expr];
 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]);
+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));
+  apply(append, lst))$
 
-coeffs : append(coeffs, β_coeffs);
+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]);
+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));
+  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]);
+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 (
+  for i : 1 thru expansion_order-1 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));
+  ev(vals, eval, diff, expand))$
 
 coeffs : append(coeffs, B_coeffs)$
 
@@ -129,42 +129,40 @@ 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);
+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]);
+e_ev_eqs[2] : expand_equation(expand(exp(2*γ)*exp(2*β)*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))
+eq_quad : block(
+  [eq : 4*solve(ce_ev_eqs[2][1][2], diff(C,u))[1]],
+  trigreduce(expand(trigsimp(eq))));
 
 
 /* Supplementary conditions */
 
 /* Supplementary Condition 1 */
-e_sc_eqs[1] : expand_equation(sc_eqs[1], coeffs);
+e_sc_eqs[1] : expand_equation(expand(exp(2*β)*sc_eqs[1]), coeffs)$
+
 eq_mono : block(
   [eq : coeff(e_sc_eqs[1], r, -2)=0],
-  trigreduce(expand(solve(eq, diff(M,u)))));
+  eq : solve(eq, diff(M,u))[1],
+  trigreduce(expand(trigsimp(eq))));
 
 /* Supplementary Condition 2 */
-e_sc_eqs[2] : expand_equation(sc_eqs[2], coeffs);
+e_sc_eqs[2] : expand_equation(expand(exp(2*β)*sc_eqs[2]), coeffs)$
+
 eq_dipo : block(
   [eq : coeff(e_sc_eqs[2], r, -2)=0],
-  trigreduce(expand(-3*solve(eq, diff(N,u)))));
-
+  eq : -3*solve(eq, diff(N,u))[1],
+  trigreduce(expand(trigsimp(eq))));
 
 
-/* BUG: should not generate terms with r^-3 and so on ... */
-/* Check if bug is related to expansion of powers */
-block([expansion_order : 2],
-  expand_equation(sc_eqs[1], coeffs))
+/* Convoluted way to show that the results for C_u are consistent */
 
-block([expansion_order : 2],
-  expand_equation(sc_eqs[2], coeffs))
+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))