First commitHEADmaster

4 files changed, 91 insertions, 105 deletions

diff --git a/README.md b/README.md
new file mode 100644
index 0000000..6daae89
--- /dev/null
+++ b/README.md
@@ -0,0 +1,10 @@
+The present source code was written to reproduce the original results from
+
+"Gravitational waves in general relativity, VII. Waves from axi-symmetric isolated system"
+1962 - Bondi, Hermann  and Van der Burg, M. G. J.  and Metzner, A. W. K.
+
+and is intended to be used with maxima. The scripts require a common lisp
+implementation supporting UTF-8 due to the use of greek characters.
+
+Maxima can be downloaded from:
+https://maxima.sourceforge.io/
diff --git a/bondi.mac b/bondi.mac
index 87b8a53..6712043 100644
--- a/bondi.mac
+++ b/bondi.mac
@@ -26,7 +26,7 @@ lg[1, 1] : V * r^-1 * exp(2 * β) - U^2 * r^2 * exp(2 * γ);
 lg[1, 2] : exp(2 * β);
 lg[2, 1] : exp(2 * β);
 
-/* dr dθ */
+/* du dθ */
 lg[1, 3] : U * r^2 * exp(2 * γ);
 lg[3, 1] : U * r^2 * exp(2 * γ);
 
@@ -53,7 +53,7 @@ uricci(false);
 sum(sum(ug[α, ε] * mcs[α, ε, 1], ε, 1, 4), α, 1, 4),expand;
 
 /** General utilities **/
-exp_taylor: taylor(exp(x), x, 0, 3);
+exp_taylor(x) := subst(t=x, taylor(exp(t), t, 0, 3));
 
 /** Main equations **/
 
@@ -69,12 +69,10 @@ 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];
-
+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));
@@ -83,82 +81,52 @@ 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 */
-sub: cons(exp(2*γ - 2*β)=exp(-2*H)*subst(x=2*(γ - β + H), exp_taylor), ansatz)$
-tmp: combine(ev(subst(sub,eq2),derivlist(r,θ),expand))$
+tmp: combine(ev(apply1(subst(ansatz, eq2), r1),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)))$
+append(ansatz, [U1=2*diff(H,θ)*exp(2*H)]);
+tmp: combine(ev(apply1(subst(%, eq3), r1),derivlist(r,θ),expand))$
 /* Solve for V_{-1} */
-solve(trigsimp(ev(subst([U1=2*diff(H,θ)*exp(2*H)], coeff(tmp2, r, 0)),derivlist(θ))), VV1);
+solve(coeff(tmp, r, 0), VV1),expand;
 
 
 /* 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;
+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;
 
-ansatz: [γ=γ_ansatz, β=β_ansatz, U=U_ansatz, V=V_ansatz];
+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))$
-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;
+combine(expand(apply1(tmp, r2)))$
+sup_cond2: -solve(coeff(%, 1/r, 2)=0, 3*diff(N,u));
diff --git a/sphere.mac b/sphere.mac
index b58cd57..f65c932 100644
--- a/sphere.mac
+++ b/sphere.mac
@@ -7,8 +7,8 @@ linel: 150;
 
 derivabbrev: true;
 
+/**  **/
 load("ctensor");
-
 dim: 2;
 cframe_flag: false;
 ct_coords: [θ, ϕ];
@@ -24,34 +24,71 @@ lg[2, 2] : sin(θ)^2;
 cmetric();
 
 /* mcs[i, j, k]: Christoffel symbols of the second kind (one upper index - k) */
-christof(false);
+christof(true);
+ricci(true);
+
+/**  **/
+load("itensor");
+
+imetric(g);
+idim(2);
+metricconvert: true;
+
 
+/**  **/
+depends([c, M, N, P], [u, θ]);
 
-depends(c, [u, θ]);
+Cdown: zeromatrix(2, 2);
+Cdown[1, 1] : 2*c;
+Cdown[2, 2] : -2*c*sin(θ)^2;
 
 Ndown: zeromatrix(2, 2);
 Ndown[1, 1] : 2*diff(c, u);
 Ndown[2, 2] : -2*diff(c, u)*sin(θ)^2;
 
-f(m, n) := sum(ug[m, a] * sum(ug[n, b] * Ndown[a, b], b, 1, 2), a, 1, 2);
-
+f(m, n, T) := sum(ug[m, a] * sum(ug[n, b] * T[a, b], b, 1, 2), a, 1, 2);
 
 Nup: zeromatrix(2, 2);
-for i: 1 thru 2 do for j: 1 thru 2 do Nup[i, j]: f(i, j);
+for i: 1 thru 2 do for j: 1 thru 2 do Nup[i, j]: f(i, j, Ndown);
+
+Cup: zeromatrix(2, 2);
+for i: 1 thru 2 do for j: 1 thru 2 do Cup[i, j]: f(i, j, Cdown);
+
 
+/** Mass-loss formula **/
+mass_aspect_eqn: Mu = (1/4)*covdiff(Nu([],[A,B]), B, A)
+-(1/8)*Nd([A,B], [])*Nu([], [A,B]);
 
 
+Mu: block([Cu,Cd,Nu,Nd,expr],
+  expr: ic_convert(mass_aspect_eqn),
+  Cu: Cup, Cd: Cdown, Nu: Nup, Nd: Ndown,
+  expand(ev(expr)));
 
 
-X: sum(sum(Nup[A,B]*Ndown[A,B], A, 1, 2), B, 1, 2);
+/** Momentum-loss formula **/
+momentum_aspect_eqn: Pu([A], []) = covdiff(M([],[]), A)
++ (1/8)*covdiff(Cu([],[B,C])*Nd([C,B],[]), A)
++ (1/4)*covdiff(covdiff(Cu([], [B, C]), B, A)
+  - g([],[C,E])*g([],[B,D])*covdiff(Cd([A,B],[]), D, E), C)
++ (1/4)*covdiff(Nu([], [B,C])*Cd([A,C],[]) - Cu([],[B,C])*Nd([A,C],[]), B)
+- (1/4)*Nu([], [B,C])*covdiff(Cd([B,C],[]), A);
 
-Q: zeromatrix(2, 1);
-for i: 1 thru 2 do Q[i]: sum(diff(Nup[i,B], θ)
-  + sum(mcs[B, D, i] * Nup[D, B], D, 1, 2)
-  + sum(mcs[B, D, B] * Nup[i, D], D, 1, 2), B, 1, 2);
 
-Y: diff(Q[1], θ) + sum(sum(mcs[A, B, A], A, 1, 2) * Q[B], B, 1, 2),expand;
+block([Cu,Cd,Nu,Nd,expr],
+  expr: ic_convert(eqn),
+  Cu: Cup, Cd: Cdown, Nu: Nup, Nd: Ndown,
+  ev(expr),
+  Pu[1]),expand;
 
+sachs_coeff: -(2/3)*g([],[A,B])*(P([B], []) + Cd([B,C],[])*Ubu([],[C]) + βb([],[],B));
+bondi_coeff: 2*N + 3*c*diff(c,θ) + 4*c^2*cot(θ);
 
-/* M_u */
-- (1/8)*X + (1/4)*Y,expand;
+Nu: block([expr, X, Cu, Cd, Nu, Nd, Ubu, βb, N],
+  expr: ic_convert(sachs_coeff=X([], [A])),
+  expr: subst([Cu=Cup, Cd=Cdown, Nu=Nup, Nd=Ndown], expr),
+  Ubu[1]: -(diff(c, θ) + 2*c*cot(θ)),
+  βb: -(1/4)*c^2,
+  expand(ev(expr)),
+  N: rhs(first(solve(X[1]=bondi_coeff, N))),
+  subst(diff(P[1],u)=Pu[1], diff(-3*N, u)));
diff --git a/tmp.mac b/tmp.mac
deleted file mode 100644
index 4bd6efe..0000000
--- a/tmp.mac
+++ /dev/null
@@ -1,29 +0,0 @@
-part('(Ndown[a, b]), );
-raise(tensor, )
-
-
-args('(Ndown[a, b]));
-
-
-raise(tensor, index) := sum(ug[''index, m] * ''tensor, m, 1, 2);
-
-raise('(Ndown[a, b]), 'a);
---> sum(ug[])
-
-:lisp (format t "~A" #$'(Ndown[a,b])$)
-
-:lisp (mfuncall '$mnctimes $ug $ug)
-
-:lisp (displa (mfuncall 'mnctimes $ug $ug));
-
-:lisp #$mcs[1,1,1]$;
-
-mcs[1,1,1]
-
-ug . ug;
-
-:lisp (format t "~A" #$c . d$)
-
-'Ndown[a, b]
-
-('Ndown[a,b])