From 34b945681740c320c3a011326d493d7bd09019f3 Mon Sep 17 00:00:00 2001
From: Thomas Guillermo Albers Raviola <thomas@thomaslabs.org>
Date: Thu, 16 Jul 2026 23:10:26 +0200
Subject: Initial commit

---
 01-bondi.mac | 165 +++++++++++++++++++++++++++++++++
 series.mac   | 296 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 setup.mac    | 105 +++++++++++++++++++++
 tensors.mac  |  21 +++++
 4 files changed, 587 insertions(+)
 create mode 100644 01-bondi.mac
 create mode 100644 series.mac
 create mode 100644 setup.mac
 create mode 100644 tensors.mac

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))
diff --git a/series.mac b/series.mac
new file mode 100644
index 0000000..5617a1e
--- /dev/null
+++ b/series.mac
@@ -0,0 +1,296 @@
+/* Term generation */
+exp_taylor(expr, n) := expand(subst([ξ=expr], taylor(exp(ξ), ξ, 0, n)));
+neg_pow_series(c, r, a, b) := sum(c[i] / r^i, i, a, b);
+
+
+/* Term filtering */
+
+filter_order(expr, x, a, b) := block([acc, len],
+  acc : 0,
+  len : if atom(expr) or not(op(expr) = "+") then 1 else length(expr),
+
+  for i thru len do block([deg, p],
+    p : if len=1 then expr else part(expr, i),
+    deg : hipow(partition([p], x)[2], x),
+    if (deg >= a and deg <= b) then acc : acc + p),
+
+  acc);
+
+/* Speed up by expanding and mapping only once over the list (check operator is +) */
+keep_upto_order(expr, var, start, end) := block([parts],
+  expr : expand(expr),
+  parts : makelist(coeff(expr, var, i) * r^i, i, start, end),
+  apply("+", parts));
+
+
+/* Term generation with filtering */
+
+/*
+Calculate the exponential series expansion evaluated at `expr' up to `n'th order
+in `x'
+*/
+exp_series(expr, x, n, a, b) := block([acc, xn, len],
+  expr : expand(expr),
+  xn : 1,
+  acc : 1,
+  len : if atom(expr) or not(is(op(expr)="+")) then 1 else length(expr),
+  /* Assume a, b < 0 for now */
+  for i thru n do block([subexpr],
+    /* Discard terms that would end up creating terms outside range */
+    /*
+    subexpr : 0,
+    for j thru len do block([term, deg],
+      term : if len=1 then expr else part(expr, j),
+      deg : hipow(term, x),
+      if not((deg + lopow(xn, x)) > b or (deg + hipow(xn, x)) < a)
+      then
+      subexpr : subexpr + term),
+    */
+    /**/
+    xn : filter_order(ratexpand(xn * expr / i), x, a, b),
+    acc : acc + xn),
+  acc);
+
+horner_exp_series(expr, x, n, a, b) := block([acc],
+  /* expand(1 + x*(1 + x*(1 + x*(1+x/4)/3)/2)) */
+  expr : expand(expr),
+  acc : 1,
+  for i thru n do acc : keep_upto_order(1 + expr*acc/(n-i+1), x, a, b),
+  acc);
+
+
+/* Multiply series expansions upto nth order */
+
+matchdeclare(m_freeof_r, lambda([e], freeof(r, e)))$
+matchdeclare(m_small, lambda([e], e < -expansion_order))$
+defrule(r_simp_small, m_freeof_r*r^m_small, 0)$
+
+__simple_filter(series, max_order) := block([ret, tmp],
+  tmp : expansion_order,
+  expansion_order : max_order,
+  ret : apply1(series, r_simp_small),
+  expansion_order : tmp,
+  ret);
+
+__table_size(factors) := block(
+  [nfactors, ranges, t_min, t_max],
+  nfactors : length(factors),
+  ranges : make_array(any, nfactors, 2),
+  /* Fill table with maximal and minimal powers */
+  for i : 0 thru nfactors-1 do block(
+    ranges[i, 0] : lopow(factors[i+1], r),
+    ranges[i, 1] : hipow(factors[i+1], r)),
+  /* Find greatest and smallest r powers */
+  t_min : ranges[0, 0],
+  t_max : ranges[0, 1],
+  for i : 0 thru nfactors-1 do (
+    t_min : min(t_min, ranges[i, 0]),
+    t_max : max(t_max, ranges[i, 1])),
+  [t_min, t_max]);
+
+__fill_term_table(table, factors, t_min, t_max, len) := block(
+  [nfactors],
+  nfactors : length(factors),
+  /* Split factors into term table */
+  for i : 0 thru nfactors-1 do (
+    /* Initialize table row with 0s */
+    for j : 0 thru len-1 do table[i, j] : 0,
+    /**/
+    if is(op(factors[i+1])="+") then (
+      for term in args(factors[i+1]) do block(
+        [j, n],
+        n : lopow(term, r),
+        j : n - t_min,
+        table[i, j] : table[i, j] + coeff(term, r, n)))
+    else block(
+      [j, n, term],
+      term : factors[i+1],
+      n : lopow(term, r),
+      j : n  - t_min,
+      table[i, j] : coeff(term, r, n))));
+
+__mul_series_level(table, max_order, i, stack, nfactors, len, t_min) := block(
+  [],
+  if is(i = nfactors) then block(
+    [order, acc],
+    order : 0,
+    for j : 0 thru nfactors-1 do order : order + stack[j] + t_min,
+    if order >= -max_order then (
+      acc : 1,
+      for j : 0 thru nfactors-1 do acc : acc * table[j, stack[j]],
+      acc * r^order)
+    else 0)
+  else block([acc],
+    /* TODO: omit external 0s */
+    acc : 0,
+    for j : 0 thru len-1 do (
+      stack[i] : j,
+      acc : acc + __mul_series_level(
+        table, max_order, i+1, stack, nfactors, len, t_min)),
+    acc));
+
+/*
+Multiply list of factors, each expressed as terms of series, only upto max_order
+*/
+mul_series(factors, max_order) := block(
+  [term_table, len, t_min, t_max, stack, nfactors],
+  nfactors : length(factors),
+
+  /* If there is nothing to multiple, still filter higher terms */
+  if is(nfactors=1) then (
+    return(__simple_filter(factors[1], max_order))),
+
+  /* Allocate stack and table for ranges */
+  stack : make_array(any, nfactors),
+
+  /* Compute table size */
+  block([pair],
+    pair : __table_size(factors),
+    t_min : pair[1],
+    t_max : pair[2]),
+  len : abs(t_max - t_min) + 1,
+
+  /* Allocate table for terms */
+  term_table : make_array(any, nfactors, len),
+  __fill_term_table(term_table, factors, t_min, t_max, len),
+
+  __mul_series_level(term_table, max_order, 0, stack, nfactors, len, t_min))$
+
+
+/**/
+
+β_expr_val(x, n) := block(
+  if x < 0 then (
+    exp(x*H)*horner_exp_series(abs(x)*(H-rhs(main_ansatz[2])), r, n, -n, 2))
+  else (
+    exp(x*H)*horner_exp_series(x*(rhs(main_ansatz[2])-H), r, n, -n, 2)));
+
+γ_expr_val(x, n) := horner_exp_series(x*rhs(main_ansatz[1]), r, n, -n, 2);
+
+matchdeclare(m_any1, lambda ([e], freeof(β, γ, e)));
+matchdeclare(m_any2, lambda ([e], freeof(β, γ, e)));
+matchdeclare(m_any3, lambda ([e], freeof(β, γ, e)));
+matchdeclare(m_beta, lambda ([e], not freeof(β, e) and freeof(γ, e)));
+matchdeclare(m_gamma, lambda ([e], not freeof(γ, e) and freeof(β, e)));
+
+/* Replace β and γ expressions into a symbolic representation */
+defrule(r_wrap_exp_β, exp(m_any1*m_beta), β_expr(m_any1));
+defrule(r_wrap_exp_γ, exp(m_any1*m_gamma), γ_expr(m_any1));
+defrule(r_wrap_exp_βγ, exp(m_any1*m_beta + m_any2*m_gamma),
+  β_expr(m_any1)*γ_expr(m_any2));
+defrule(r_wrap_exp_βγ2, exp(m_any1*m_beta + m_any2*m_gamma + m_any3),
+  exp(m_any3)*β_expr(m_any1)*γ_expr(m_any2));
+
+/* Apply wrap rules to expression */
+wrap_exp(expr) := (
+  apply1(expr, r_wrap_exp_β, r_wrap_exp_γ, r_wrap_exp_βγ, r_wrap_exp_βγ2));
+
+/* Replace symbolic representations of β and γ expressions into their their
+expansions */
+defrule(r_exp_β_val, β_expr(m_any1), β_expr_val(m_any1, expansion_order));
+defrule(r_exp_γ_val, γ_expr(m_any1), γ_expr_val(m_any1, expansion_order));
+
+
+/**/
+
+freeof_many(expr, funcs) := apply(freeof, append(funcs, [expr]));
+
+/* Split product in factors containing any funcs and factors free of them */
+split_prod(expr, funcs) := block(
+  [deps_on_funcs : []],
+  /* If not a product, assume product with single factor */
+  for f in funcs do block(
+    if atom(expr) then (
+      if member(expr, funcs) then (
+        push(expr, deps_on_funcs),
+        expr : 1),
+      return())
+    elseif op(expr)="*" then block(
+      [split : partition(expr, f)],
+      expr : split[1],
+      push(split[2], deps_on_funcs))
+    elseif not freeof_many(expr, funcs) then (
+      push(expr, deps_on_funcs),
+      expr : 1,
+      return())
+    else (
+      return())),
+  [expr, apply("*", deps_on_funcs)]);
+
+/* Substitute functions by their series expansions and apply substitutions in
+coeffs */
+__expand_term_subst(expr, coeffs) := block(
+  expr : apply1(expr, r_exp_β_val, r_exp_γ_val),
+  expr : subst(append(main_ansatz, coeffs), expr),
+  ev(expr, eval, diff, expand));
+
+/* Substitute product of functions by result of mul_series */
+__expand_term_prod(expr, coeffs) := block(
+  [free, factors],
+  expr : wrap_exp(expr),
+  /* Split product into list of factors containing funcs and those independent
+  of funcs*/
+  block([split : split_prod(expr, funcs)],
+    free : split[1],
+    if not atom(split[2]) and op(split[2])="*" then (
+      factors : args(split[2]))
+    else (
+      factors : [split[2]])),
+  /* Substitute funcs factors by expansion */
+  factors : __expand_term_subst(factors, coeffs),
+  /* Total product using series optimization */
+  free * mul_series(factors, expansion_order + lopow(free, r)));
+
+/*
+Split product factors into list substitute series expansions, substitute
+coefficients and finally multiply only upto expansion_order.
+*/
+expand_term(expr, coeffs) := block(
+  [tmp],
+  /* TODO: Recursive call instead of error */
+  if op(expr)="+" then (
+    error("Term is itself a sum!"))
+  /* For single negation just call recursively */
+  elseif op(expr)="-" then (
+    if length(args(expr)) # 1 then (
+      error("Unexpected '-' operator with more than 1 argument")),
+    expand(- expand_term(part(expr, 1), coeffs)))
+  elseif op(expr)="*" then (
+    expand(__expand_term_prod(expr, coeffs)))
+  elseif op(expr)="/" then block(
+    [n, d, expansion_order:expansion_order - lopow(d, r)],
+    n : part(expr, 1),
+    d : part(expr, 2),
+    expand(expand_term(n, coeffs) / d))
+  /* TODO: better handle exponents */
+  else (
+    apply1(__expand_term_subst(wrap_exp(expr), coeffs), r_simp_small)));
+
+
+/* Solve main equations term by term */
+
+/* TODO: Divide and conquer on term list for long expansions */
+/* coeffs: facts already known about series expansions */
+
+/* Series expansion inside each instance of a value in the terms of an equation
+*/
+expand_equation(expr, coeffs) := block(
+  [len, terms],
+  expr : rhs(expr) - lhs(expr),
+  0 = if atom(expr) then (
+    expr)
+  elseif op(expr)="+" then (
+    apply("+", map(lambda([term], expand_term(term, coeffs)), args(expr))))
+  else (
+    expand_term(expr, coeffs)))$
+
+
+/* Misc */
+
+/* For an equation of polynomials, returns the term by term list of equations of
+the coefficients */
+coeff_equations(eq) := block([start, end],
+  eq : expand(eq),
+  start : lopow(eq, r),
+  end : hipow(eq, r),
+  reverse(makelist([r^i, coeff(eq, r, i)], i, start, end)));
diff --git a/setup.mac b/setup.mac
new file mode 100644
index 0000000..becfdc3
--- /dev/null
+++ b/setup.mac
@@ -0,0 +1,105 @@
+/*
+ * Setup
+ */
+linel: 120;
+
+derivabbrev: true;
+
+load("ctensor");
+
+dim: 4;
+cframe_flag: false;
+ct_coords: [u, r, θ, ϕ];
+
+/* Original Bondi U and V are replaced by A and B; V is reserved for potential */
+depends([A, B, β, γ], [u, r, θ]);
+
+/*
+Bondi metric ansatz
+NOTE: Signature is different than original Bondi paper
+*/
+bondi_ansatz() := block([m],
+  m: zeromatrix(4, 4),
+  /* du du */
+  m[1, 1] : -B * r^-1 * exp(2 * β) + A^2 * r^2 * exp(2 * γ),
+  /* du dr */
+  m[1, 2] : -exp(2 * β),
+  m[2, 1] : -exp(2 * β),
+  /* du dθ */
+  m[1, 3] : -A * r^2 * exp(2 * γ),
+  m[3, 1] : -A * r^2 * exp(2 * γ),
+  /* dθ dθ */
+  m[3, 3] : r^2 * exp(2 * γ),
+  /* dϕ dϕ */
+  m[4, 4] : r^2 * exp(- 2 * γ) * sin(θ)^2,
+  m)$
+
+lg: bondi_ansatz()$
+
+/* 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);
+
+/* Ricci scalar */
+ric_s : expand(scurvature());
+
+
+/* 2-sphere submanifold */
+
+/* 2-sphere submanifold metric */
+h_AB[3,3] : exp(2*γ)$
+h_AB[4,4] : exp(-2*γ)*sin(θ)^2$
+h_AB[3,4] : 0$
+h_AB[4,3] : 0$
+
+/* 2-sphere submanifold inverse metric */
+uh_AB[3,3] : exp(-2*γ)$
+uh_AB[4,4] : exp(2*γ)*sin(θ)^-2$
+uh_AB[3,4] : 0$
+uh_AB[4,3] : 0$
+
+/* 2-sphere metric */
+q_AB[3,3] : 1$
+q_AB[4,4] : sin(θ)^2$
+q_AB[3,4] : 0$
+q_AB[4,3] : 0$
+
+/* 2-sphere inverse metric */
+uq_AB[3,3] : 1$
+uq_AB[4,4] : sin(θ)^-2$
+uq_AB[3,4] : 0$
+uq_AB[4,3] : 0$
+
+/* Christoffel symbols of 2-sphere */
+for i : 3 thru 4 do for j : 3 thru 4 do for k : 3 thru 4 do (
+  h_mcs[i, j, k] : sum(uh_AB[k,C]*(
+      cdiff(h_AB[j, C], i)
+      + cdiff(h_AB[i, C], j)
+      + cdiff(h_AB[i, j], C))/2, C, 3, 4));
+
+/* Riemann tensor of 2-sphere */
+for i : 3 thru 4 do for j : 3 thru 4 do for k : 3 thru 4 do for l : 3 thru 4 do (
+  h_riem[i,j,k,l] : expand(
+    cdiff(h_mcs[i,k,l], j)
+    - cdiff(h_mcs[j,k,l], i)
+    + sum(h_mcs[i,k,n]*h_mcs[l,n,j] - h_mcs[j,k,n]*h_mcs[n,i,l], n, 3, 4)));
+
+/* Ricci tensor of 2-sphere */
+for i : 3 thru 4 do for j : 3 thru 4 do (
+  h_ric[i,j] : expand(sum(h_riem[i,k,j,k], k, 3, 4)));
+
+/* Ricci scalar of 2-sphere */
+h_ric_s : expand(sum(sum(uh_AB[i,j]*h_ric[i,j], i, 3, 4), j, 3, 4));
+
+print("Setup done!");
diff --git a/tensors.mac b/tensors.mac
new file mode 100644
index 0000000..dea722a
--- /dev/null
+++ b/tensors.mac
@@ -0,0 +1,21 @@
+/* Contract indices in expression */
+
+/* NOTE: buildq works like splicing a list with quasiquote in lisp; escaped
+items are defined at start of quote */
+contr(expr, var) ::= (
+  if listp(var) and length(var) > 1 then
+  block([head, tail],
+    head : first(var),
+    tail : rest(var),
+    buildq([expr, head, tail], sum(contr(expr, tail), head, 1, 4)))
+  elseif listp(var) and length(var) = 1 then
+  block([head],
+    head : first(var),
+    buildq([expr, head], sum(expr, head, 1, 4)))
+  elseif listp(var) then
+  expr
+  else
+  buildq([expr, var], sum(expr, var, 1, 4)))$
+
+/* Differentiate with respect to nth coordinate */
+cdiff(func, coord) := diff(func, ct_coords[coord])$
-- 
cgit v1.3