2024-08-04

3 files changed, 229 insertions, 16 deletions

diff --git a/bondi.mac b/bondi.mac
index 3577589..87b8a53 100644
--- a/bondi.mac
+++ b/bondi.mac
@@ -1,37 +1,164 @@
-debugmode: true;
-
+/* 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, theta, phi];
+ct_coords: [u, r, θ, ϕ];
 
-depends([U, V, _beta, _gamma], [u, r, theta]);
+depends([U, V, β, γ], [u, r, θ]);
 
 lg : zeromatrix(4, 4);
-lg[1, 1] : V * r^-1 * exp(2 * _beta) - U^2 * r^2 * exp(2 * _gamma);
+/* du du */
+lg[1, 1] : V * r^-1 * exp(2 * β) - U^2 * r^2 * exp(2 * γ);
 
-lg[1, 2] : exp(2 * _beta);
-lg[2, 1] : exp(2 * _beta);
+/* du dr */
+lg[1, 2] : exp(2 * β);
+lg[2, 1] : exp(2 * β);
 
-lg[1, 3] : U * r^2 * exp(2 * _gamma);
-lg[3, 1] : U * r^2 * exp(2 * _gamma);
+/* dr dθ */
+lg[1, 3] : U * r^2 * exp(2 * γ);
+lg[3, 1] : U * r^2 * exp(2 * γ);
 
-lg[3, 3] : - r^2 * exp(2 * _gamma);
-lg[4, 4] : - r^2 * exp(- 2 * _gamma) * sin(theta)^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();
 
-/* For some reason some components must be simplified by hand ... */
+/* 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; */
+
 /*
-    Calculate Christoffel symbols
-    mcs[i, j, k]: Christoffel symbols of the second kind (one upper index - k)
+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))
 */
 
-christof(false);
-ricci(false);
+
+/** 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;
diff --git a/sphere.mac b/sphere.mac
new file mode 100644
index 0000000..b58cd57
--- /dev/null
+++ b/sphere.mac
@@ -0,0 +1,57 @@
+/* 1920x1080 - HD screen */
+linel: 230;
+/* 3840x1600 - Ultrawide screen */
+linel: 470;
+/* 1366x768 - Poor man's thinkpad */
+linel: 150;
+
+derivabbrev: true;
+
+load("ctensor");
+
+dim: 2;
+cframe_flag: false;
+ct_coords: [θ, ϕ];
+
+lg : zeromatrix(2, 2);
+/* dθ dθ */
+lg[1, 1] : 1;
+
+/* dϕ dϕ */
+lg[2, 2] : sin(θ)^2;
+
+/* Load metric and compute inverse */
+cmetric();
+
+/* mcs[i, j, k]: Christoffel symbols of the second kind (one upper index - k) */
+christof(false);
+
+
+depends(c, [u, θ]);
+
+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);
+
+
+Nup: zeromatrix(2, 2);
+for i: 1 thru 2 do for j: 1 thru 2 do Nup[i, j]: f(i, j);
+
+
+
+
+
+X: sum(sum(Nup[A,B]*Ndown[A,B], A, 1, 2), B, 1, 2);
+
+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;
+
+
+/* M_u */
+- (1/8)*X + (1/4)*Y,expand;
diff --git a/tmp.mac b/tmp.mac
new file mode 100644
index 0000000..4bd6efe
--- /dev/null
+++ b/tmp.mac
@@ -0,0 +1,29 @@
+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])