psi:=Sqrt[R[r]/r]

beta1:=psi^4* x*C/R[r]^3

beta2:=psi^4* y*C/R[r]^3

beta3:=psi^4 *z*C/R[r]^3

alpha:=Sqrt[1-2*M/R[r]+C^2/R[r]^4]

gamma:=1/Sqrt[1-beta^2]

gorig:={{-alpha^2+psi^(-4)*(beta1^2+beta2^2+beta3^2), beta1, beta2, beta3},{beta1,psi^4,0,0},{beta2,0,psi^4,0},{beta3,0,0,psi^4}}

Lambda:={{gamma,-gamma*beta,0,0},{-gamma*beta,gamma,0,0},{0,0,1,0},{0,0,0,1}}

g:=FullSimplify[Lambda.gorig.Lambda]

g3d:=g[[2;;4,2;;4]]

g3u:=FullSimplify[Inverse[g3d]]

g3dy:=FullSimplify[D[g3d,y,NonConstants->{r}]/.{R'[r]->R[r]*alpha/(r),D[r,y,NonConstants->{r}]->y/r}]

g3dz:=FullSimplify[D[g3d,z,NonConstants->{r}]/.{R'[r]->R[r]*alpha/(r),D[r,z,NonConstants->{r}]->z/r}]

g3dx:=FullSimplify[gamma*D[g3d,x,NonConstants->{r}]/.{R'[r]->R[r]*alpha/(r),D[r,x,NonConstants->{r}]->x/r}]

g3diff:={g3dx,g3dy,g3dz}

bbeta:=g[[1,2;;4]]

aalpha:=FullSimplify[Sqrt[-1/Inverse[g][[1,1]]]]

G := Table[Sum[1/2*g3u[[k,l]]*(g3diff[[i]][[j,l]]+g3diff[[j]][[i,l]]-g3diff[[l]][[i,j]])*bbeta[[k]],{k,3},{l,3}],{i,3},{j,3}]

Dbeta = FullSimplify[Table[D[bbeta[[i]],j,NonConstants->{r}]/.{R'[r]->R[r]*alpha/(r),D[r,z,NonConstants->{r}]->z/r,D[r,x,NonConstants->{r}]->x/r,D[r,y,NonConstants->{r}]->y/r},{i,3},{j,{x,y,z}}]]

Dbeta[[All,1]]=Dbeta[[All,1]]*gamma

K:=Table[1/(2*aalpha)*(beta*g3dx[[i,j]]+Dbeta[[i,j]]-G[[i,j]]+Dbeta[[j,i]]-G[[j,i]]),{i,3},{j,3}]

KK=FullSimplify[K]/.{r^2->x^2+y^2+z^2}

FullSimplify[KK[[1,3]]/.{beta->0}]

CForm[FullSimplify[KK[[2,3]]]]/.{R[r]->R,M->MASS}

Sum[KK[[i,j]]g3u[[i,j]],{i,3},{j,3}]-Sum[(g3u.KK.g3u)[[i,j]]KK[[i,j]],{i,3},{j,3}]