/* 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;