(* ::Package:: *) (* Author: Jaime Varela Organization: University of California Berkeley *) (* This package is used to calculate various expressions in general relativity in coordinate form. The conventions used here follow S. Caroll. NOTE: in the following all indeces go from 0 to 3.*) BeginPackage["GRQUICK`"] (* This provides the usage of the functions *) Metin::usage = "For a dxd matrix g and d dimensional vector vec, Metin[g,vec] sets the working Metric to g and the coordinates to vec (default {t,x,y,z}). One can also use default metrics by inputing certain strings. For a list of metrics type ?DefaultMet" DefaultMet::usage = "The load decault metrics use Metin[string] with one of the following strings\n\n Schwar = Schwarchild metric with COORDINATES={t,r,\[Theta],\[Phi]} for a central mass of mass M;\n FRW = Friedman walker metric with curvature k and coordinates {t,r,\[Theta],\[Phi]} and expansion function a[t];\n Reissner = Charged black hole with total charge q;\n Kerr = rotating black hole with total angular momentum J;\n Kerr-Newman = charged rotating black hole with charge q and angular momentum J;\n DS static = Static desitter coordinates {t,r,\[Theta],\[Phi]} with the cosmological horizon at r=\[Alpha] = \!\(\*SqrtBox[\(3/\[CapitalLambda]\)]\); \n DS universe = Desitter universe with a[t]=\!\(\*SuperscriptBox[\(e\), \(H\\\ t\)]\) where H is the hubble constant;\n DS-Schwar = De-sitter Schwarzchild metric where \[CapitalLambda] is the cosmological constant and M is the mass of the black hole.\n\n Note: we use geometrical units in which G=c=1. The metrics are taken from Wikipedia and I make no claim to their accuracy." GINV::usage = "GINV, provides the inverse Metric in an array" LoadMetric::usage = "Loads defaul metrics (internal use only)" Metric::usage = "Metric[i,j] gives the ith jth component of the Metric" IMetric::usage = "IMetric[i,j] gives the ith jth inverse component of the Metric" Christoffel::usage = "Christoffel[{{i}, {j,k}}] provides (\!\(\*SuperscriptBox[\(\[CapitalGamma]\), \(i\)]\)\!\(\*SubscriptBox[\()\), \(jk\)]\)" Riemann::usage = "Riemann[{{l},{i,j,k}}] = (\!\(\*SuperscriptBox[\(R\), \(l\)]\)\!\(\*SubscriptBox[\()\), \(ijk\)]\) or Riemann[i,j,k,l]=Riemann[{i,j,k,l}]=\!\(\*SubscriptBox[\(R\), \(ijkl\)]\)" Ricci::usage = "The Ricci tensor Ricci[i,j] = \!\(\*SubscriptBox[\(R\), \(\(ij\)\(\\\\n\)\)]\)" Rscalar::usage = "R = \!\(\*SubscriptBox[SuperscriptBox[\(R\), \(u\)], \(u\)]\)" Einstein::usage = "Einstein[i,j] = \!\(\*SubscriptBox[\(R\), \(ij\)]\)" (* The follwing functions do not yield expressions which are stored in memory *) Geodesics::usage = "After the christofell symbols have been calculated, Geodesics[a], will display the geodesic equations x''[a\!\(\*SuperscriptBox[\(]\), \(i\)]\)+\!\(\*SubscriptBox[SuperscriptBox[\(\[CapitalGamma]\), \(i\)], \(jk\)]\) x'[a\!\(\*SuperscriptBox[\(]\), \(j\)]\) x'[a\!\(\*SuperscriptBox[\(]\), \(j\)]\) as a vector. " Tenstr::usage = "Tenstr[TENSOR] will give the trace of a rank two lower index tensor. (Input is a string) For example Tenstr[Einstein] gives \!\(\*SubscriptBox[SuperscriptBox[\(G\), \(i\)], \(i\)]\) . Note Tenstr can only intake GRQUICK rank 2 tensors, not user defined tensors." helpGRQUICK::usage = "To begin one must input a metric using Metin. Type in ?Metin for usage. For the usage of any function use ?function Standard Differential Geometry objects are: Metric, IMetric, Christoffel,Riemann, Weyl,Ricci,Rscalar,Einstein, Kretschmann, Metprod, Metnorm, Normalizevec, RaiseVector, LowerVector, and Setlight. Functions usefull for geodesic analysis are: Geodesics, Geoplot, Issol, ReplaceCoords Functions which operate on tensore are: Coderv, and Tenstr. Raychaudhuri functions are: SpatialMetric, BExpansion, ThetaExpansion, Shear, and Twist. The vector inputs to these function must be properly timelike normalized (I have not implemented a null formulation as of yet). For definitions see Wald's General Relativity text. Important objects/tensors can be displayed using the commands: (GINV; GMET; DIM; COORDINATES) ,CHARRAY, GARRAY, RARRAY, RSARRAY, SPEEDOFLIGHT, and GEOSOL. These object are populated once the functions (Metin), Christoffel, Einstein, Ricci, Rscalar, Setlight, and Geoplot have been computed respectively." GMET::usage = "the Metric tensor" CHARRAY::usage = "an array of Christofell symbols" GARRAY::usage = "array of the einstein tensor" RARRAY::usage = "array of the ricci tensor" RSARRAY::usage = "array of the ricci tensor" COORDINATES::usage = "the coordinate vector" DIM::usage = "space-time dimension. DIM-1 is the spatial dim" SPEEDOFLIGHT::usage = "value of the speed of light (default c=1)" setzero::usage = "to reset values" constructnondummy::usage = "used for indeces" Coderv::usage = "Coderv[string,{{indtop},{indbot}},a] results in the ath covariant derivative of the tensor string that has the index structure {{indtop},{indbot}}. For example Coderv[Ricci,{1,2},0] the 0th covariant derivative of R_{12}"; Setlight::usage = "Setlight[val] will set the numerical speed of light to be c=val (default c=1). This is needed for numerical functions" Geoplot::usage = "Geoplot[t,init,vinit,range,plotvars] uses NDSOLVE to plot the timelike (or null) Geodesics for a system with parameter t, an initial position vector (init), an initial four velocity (vinit), and a two-vector range={ti,tf} specifying the range. The four velocity must be properly normalized and the initial conditions are for t=ti. (remember to change the value to the speed of light to reflect the units in your Metric)." GEOSOL::usage = "after Geoplot has been called once, GEOSOL contains the interpolation solutions to the previous Geoplot call." Metnorm::usage = "given a CONTRAVARIANT vector, Metnorm[vec] gives its norm, namely g_{i,j} vec^i vec^j" Metprod::usage = "product of two CONTRAVARIANT vectors, Metprod[vec1,vec2] gives the scalar product in curved spacetime, namely vec1^i g_{i,j} vec2^j." Normalizevec::usage = "Normalizevec[X], for any D-dimensional array it returns a normalize vector X/Sqrt[-Metprod[X]]" Ricciint::usage = "for programs use only, used as an intermediate Ricci tensor" Einsint::usage = "for programs use only, used as an intermediate Einstein" Issol::usage = "Issol[x[\[Tau]],\[Tau]] will give True if the four vector x[\[Tau]] is a solution to the geodesic equations and False otherwise. Issol may also give an expression which must be zero for \!\(\*SuperscriptBox[\(x\), \(\[Mu]\)]\) to be a geodesic" ReplaceCoords::usage = "ReplaceCoords[exp,coords] replaces the default coordinates (COORDINATES) to thos specified in coords." RaiseVector::usage = "Give an vector X={x0,x1,...}, RaiseVector[X] will treat X as covariant vector and return the contravariant vector" LowerVector::usage = "Give an vector X={x0,x1,...}, RaiseVector[X] will treat X as contravariant vector and return the covariant vector" SpatialMetric::usage = "SpatialMetric[{\[Mu],\[Nu]},{\!\(\*SubscriptBox[\(\[Zeta]\), \(0\)]\),\!\(\*SubscriptBox[\(\[Zeta]\), \(1\)]\),\!\(\*SubscriptBox[\(\[Zeta]\), \(2\)]\),...}] gives \!\(\*SubscriptBox[\(h\), \(\[Mu]\[Nu]\)]\)=\!\(\*SubscriptBox[\(g\), \(\[Mu]\[Nu]\)]\)+\!\(\*SubscriptBox[\(\[Zeta]\), \(\[Mu]\)]\)\!\(\*SubscriptBox[\(\[Zeta]\), \(\[Nu]\)]\)" BExpansion::usage = "BExpansion[{\[Mu],\[Nu]},{\!\(\*SubscriptBox[\(\[Zeta]\), \(0\)]\)(x),\!\(\*SubscriptBox[\(\[Zeta]\), \(1\)]\)(x),\!\(\*SubscriptBox[\(\[Zeta]\), \(2\)]\)(x),\!\(\*SubscriptBox[\(\[Zeta]\), \(3\)]\)(x)}] give \!\(\*SubscriptBox[\(B\), \(\[Mu]\[Nu]\)]\)=\!\(\*SubscriptBox[\(\[Del]\), \(\[Nu]\)]\)\!\(\*SubscriptBox[\(\[Zeta]\), \(\[Mu]\)]\)" ThetaExpansion::usage = "ThetaExpansion[{\!\(\*SubscriptBox[\(\[Zeta]\), \(0\)]\)(x),\!\(\*SubscriptBox[\(\[Zeta]\), \(1\)]\)(x),\!\(\*SubscriptBox[\(\[Zeta]\), \(2\)]\)(x),\!\(\*SubscriptBox[\(\[Zeta]\), \(3\)]\)(x)}] gives the Raychauduri scalar \[Theta]" Shear::usage = "Shear[{\[Mu],\[Nu]},{\!\(\*SubscriptBox[\(\[Zeta]\), \(0\)]\)(x),\!\(\*SubscriptBox[\(\[Zeta]\), \(1\)]\)(x),\!\(\*SubscriptBox[\(\[Zeta]\), \(2\)]\)(x),\!\(\*SubscriptBox[\(\[Zeta]\), \(3\)]\)(x)}] gives \!\(\*SubscriptBox[\(\[Sigma]\), \(\[Mu]\[Nu]\)]\)=\!\(\*SubscriptBox[\(B\), \([a\\\ b]\)]\)- \!\(\*FractionBox[\(1\), \(3\)]\)\[Theta] \!\(\*SubscriptBox[\(h\), \(\[Mu]\[Nu]\)]\)" Twist::usage = "Twist[{\[Mu],\[Nu]},{\!\(\*SubscriptBox[\(\[Zeta]\), \(0\)]\)(x),\!\(\*SubscriptBox[\(\[Zeta]\), \(1\)]\)(x),\!\(\*SubscriptBox[\(\[Zeta]\), \(2\)]\)(x),\!\(\*SubscriptBox[\(\[Zeta]\), \(3\)]\)(x)}] gives \!\(\*SubscriptBox[\(\[Omega]\), \(\[Mu]\[Nu]\)]\)=\!\(\*SubscriptBox[\(B\), \([a\\\ b]\)]\)" Kretschmann::usage = "Kretschmann invariant" Weyl::usage = "The Weyl tensor given as defined in Carrol is given by Weyl[ijkl]=Weyl[{ijkl}]= \!\(\*SubscriptBox[\(C\), \(ijkl\)]\)" Begin["Global`"] LoadMetric[met_]:=Block[{g,v}, g=101; v=0; Switch[met, "Schwar", {g=DiagonalMatrix[{-(1-(2 M/r)),1/(1-(2 M/r)),r*r, r^2 Sin[\[Theta]]^2}] ;v={t,r,\[Theta],\[Phi]}} ,"FRW", {g=DiagonalMatrix[{-1,a[t]^2/(1-k r^2),a[t]^2 r^2,a[t]^2 r^2 Sin[\[Theta]]^2}] ;v={t,r,\[Theta],\[Phi]}} ,"DS static", {g=DiagonalMatrix[{-(1-r^2/\[Alpha]^2),1/(1-r^2/\[Alpha]^2),r^2,r^2 Sin[\[Theta]]^2}] ;v={t,r,\[Theta],\[Phi]}} ,"DS universe", {g=DiagonalMatrix[{-1,Exp[H t]^2/(1-k r^2),Exp[H t]^2 r^2,Exp[H t]^2 r^2 Sin[\[Theta]]^2}] ;v={t,r,\[Theta],\[Phi]}} ,"DS-Schwar", {g=DiagonalMatrix[{-(1-( \[CapitalLambda] r^2)/3-(2 M)/r),1/(1-( \[CapitalLambda] r^2)/3-(2 M)/r),r^2,r^2 Sin[\[Theta]]^2}] ;v={t,r,\[Theta],\[Phi]}} ,"Reissner", {g=DiagonalMatrix[{-(1-(2 M)/r+q^2/r^2),1/(1-(2 M)/r+q^2/r^2),r^2,r^2 Sin[\[Theta]]^2}] ;v={t,r,\[Theta],\[Phi]}} ,"Kerr", {g= {{-1+(2 M r)/(r^2+(J^2 Cos[\[Theta]]^2)/M^2),0,0,-((2 J r Sin[\[Theta]]^2)/(r^2+(J^2 Cos[\[Theta]]^2)/M^2))},{0,(M^2 r^2+J^2 Cos[\[Theta]]^2)/(J^2+M^2 r (-2 M+r)),0,0}, {0,0,r^2+(J^2 Cos[\[Theta]]^2)/M^2,0}, {-((2 J r Sin[\[Theta]]^2)/(r^2+(J^2 Cos[\[Theta]]^2)/M^2)),0,0,Sin[\[Theta]]^2 (J^2/M^2+r^2+(2 J^2 M r Sin[\[Theta]]^2)/(M^2 r^2+J^2 Cos[\[Theta]]^2))}} ;v={t,r,\[Theta],\[Phi]}} ,"Kerr-Newman", {g= {{-1+(2 M^2 (-q^2+2 M r))/(J^2+2 M^2 r^2+J^2 Cos[2 \[Theta]]),0,0,(J M (q^2-2 M r) Sin[\[Theta]]^2)/(M^2 r^2+J^2 Cos[\[Theta]]^2)}, {0,(r^2+(J^2 Cos[\[Theta]]^2)/M^2)/(J^2/M^2+q^2+r (-2 M+r)),0,0},{0,0,r^2+(J^2 Cos[\[Theta]]^2)/M^2,0}, {(J M (q^2-2 M r) Sin[\[Theta]]^2)/(M^2 r^2+J^2 Cos[\[Theta]]^2),0,0,(Sin[\[Theta]]^2 ((J^2+M^2 r^2)^2-J^2 (J^2+M^2 (q^2-2 M r+r^2)) Sin[\[Theta]]^2))/(M^4 r^2+J^2 M^2 Cos[\[Theta]]^2)}} ;v={t,r,\[Theta],\[Phi]}} ]; If[g==101,Return[{"Not a default metric","Not a default metric"}]]; Return[{g,v}]; ] End[] Begin["`Private`"] (*The program sets GMET as a global matrix *) (*Default GMET = diag(-1,1,1,1) *) GMET = {{-1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}}; GINV= {{-1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}}; (*Default tensor arrays which will store expressions for later use*) CHARRAY = 0; RARRAY = 0; RSARRAY = Null; GARRAY = 0; COORDINATES = {t,x,y,z}; DIM=4; (*Sets the global coordinates*) (** The functions **) (*Resets default values*) setzero[a_]:=Block[{}, CHARRAY = 0; RARRAY = 0; RSARRAY = Null; GARRAY = 0; SPEEDOFLIGHT=1; Return[Print["Metric Success"]]; ] (*Inputs the Metric and coordinates*) Metin[gin__]:= Block[{g,V,vec}, vec={gin}; If[StringQ[vec[[1]]]==False, g=vec[[1]]; V=vec[[2]];, g=LoadMetric[vec[[1]]][[1]]; V=LoadMetric[vec[[1]]][[2]]; ]; DIM = Length[V]; (*If[StringQ[g]==True,Return[g]];*) If[Transpose[g]-g != IdentityMatrix[DIM]-IdentityMatrix[DIM], Return[Print["Invalid Metric"]], GMET = g; GINV=Simplify[Inverse[GMET]]; COORDINATES = V; DIM = Length[COORDINATES]; setzero[0]; ]] Metric[i_,j_]:=Return[GMET[[i+1,j+1]]] IMetric[i_,j_]:= Return[GINV[[i+1,j+1]]] Metnorm[vec_]:= Sum[vec[[ii]] Metric[ii-1,jj-1] vec[[jj]], {ii,1,DIM},{jj,1,DIM}] Metprod[vec1_,vec2_]:= Sum[vec1[[ii]] Metric[ii-1,jj-1] vec2[[jj]], {ii,1,DIM},{jj,1,DIM}] Normalizevec[vec_]:= Block[{a}, If[Length[vec] != Length[COORDINATES], Print["Vector is of incorrect size"]]; Return[vec/Sqrt[-Metnorm[vec]]]; ] Christoffel[{{i_},{k_,l_}}]:=Block[{exp}, If[Length[CHARRAY] !=0, exp = CHARRAY[[i+1]][[k+1,l+1]], (*populates the CHARRAY if this has not been done*) CHARRAY = Table[Table[Simplify[(1/2) Sum[(GINV[[is,var]]) (D[GMET[[var,ks]],COORDINATES[[ls]]] + D[GMET[[var,ls]],COORDINATES[[ks]]] +-D[GMET[[ks,ls]],COORDINATES[[var]]]), {var,1,Length[COORDINATES]}]], {ks,1,Length[COORDINATES]}, {ls,1,Length[COORDINATES]}], {is,1, Length[COORDINATES]}]; exp = CHARRAY[[i+1]][[k+1,l+1]]; ]; Return[exp] ] Riemann[ind__]:= Block[{vec,l,i,j,k}, vec= Quiet[Check[Length[ind]*ind,4{ind}]]; If[Length[CHARRAY]==0, Christoffel[{{0},{0,0}}]]; If[Length[vec]==2, vec=vec/2; {l=vec[[1]][[1]],i=vec[[2]][[1]],j=vec[[2]][[2]],k=vec[[2]][[3]]}; Return[ Simplify[ D[(CHARRAY[[l+1]][[i+1,k+1]]),COORDINATES[[j+1]] ] -D[(CHARRAY[[l+1]][[i+1,j+1]]),COORDINATES[[k+1]] ] +Sum[(CHARRAY[[l+1]][[j+1,ss]]) (CHARRAY[[ss]][[i+1,k+1]]) ,{ss,1 Length[COORDINATES]}] -Sum[(CHARRAY[[l+1]][[k+1,ss]]) (CHARRAY[[ss]][[i+1,j+1]]) ,{ss,1 Length[COORDINATES]}] ] ] , vec=vec/4; {l=vec[[1]],i=vec[[2]],j=vec[[3]],k=vec[[4]]}; Return[ Simplify[1/2 (D[Metric[l,k],COORDINATES[[i+1]],COORDINATES[[j+1]]]+D[Metric[i,j],COORDINATES[[l+1]],COORDINATES[[k+1]]] -D[Metric[l,j],COORDINATES[[i+1]],COORDINATES[[k+1]]]-D[Metric[i,k],COORDINATES[[l+1]],COORDINATES[[j+1]]])+ Sum[Metric[nn,pp](Christoffel[{{nn},{i,j}}]Christoffel[{{pp},{l,k}}]-Christoffel[{{nn},{i,k}}]Christoffel[{{pp},{l,j}}]) ,{nn,0,Length[COORDINATES]-1},{pp,0, Length[COORDINATES]-1}] ] ] ]] Ricciint[i_,j_]:= Block[{sexp}, If[Length[RARRAY]==0, RARRAY = Table[Simplify[ Sum[Riemann[{{ls}, {is, ls, js}}],{ls,0,DIM-1}]] ,{is,0,DIM-1},{js,0,DIM-1}]; ]; sexp = RARRAY[[i+1,j+1]]; Return[sexp] ] Ricci[m__]:=Block[{check}, check = Quiet[Check[0*(m[[1]]+m[[2]]),1]]; If[check==0, Return[Ricciint[m[[1]],m[[2]]] ] ,Return[Ricciint[m]] ] ] Rscalar:=Block[{}, If[Length[RARRAY]==0, Ricci[0,0]]; If[RSARRAY==Null ,RSARRAY = Simplify[Tr[GINV.RARRAY]]]; Return[RSARRAY]] Einsint[i_,j_]:=Block[{}, If[Length[RARRAY]==0,Ricci[0,0]]; If[Length[GARRAY]==0, GARRAY = Simplify[RARRAY - (Rscalar/2) GMET ]]; Return[GARRAY[[i+1,j+1]]] ] Einstein[m__]:=Block[{check}, check = Quiet[Check[0*(m[[1]]+m[[2]]),1]]; If[check==0, Return[Einsint[m[[1]],m[[2]]] ] ,Return[Einsint[m]] ] ] Weyl[ind__]:= Block[{vec,i,k,l,m}, vec= Quiet[Check[Length[ind]*ind,4{ind}]]; vec=vec/4; If[Length[CHARRAY]==0, Christoffel[{{0},{0,0}}]]; {i=vec[[1]],k=vec[[2]],l=vec[[3]],m=vec[[4]]}; Return[Riemann[{i,k,l,m}]+ (1/(DIM-2))(-Ricci[{i,l}] Metric[k,m]+Ricci[i,m] Metric[k,l]+Ricci[k,l] Metric[i,m]-Ricci[k,m] Metric[i,l]) +Rscalar/((DIM-1)(DIM-2)) (Metric[i,l]Metric[k,m]-Metric[i,m]Metric[k,l])] ] Coderv[expi_,ind_,a_]:=Block[ {topi,boti,p,q,rank,dummy,codv,int,indvec,inter,interind,vector,exp,derv,codvp,codvq}, If[Length[CHARRAY]==0, Christoffel[{{0},{0,0}}]]; If[StringQ[expi]==False && Length[expi]==0, exp=ToString[expi],exp=expi]; If[Length[ind[[1]] ]==0,topi=0,topi=ind[[1]]+1]; If[Length[topi]==0,boti=ind+1,boti=ind[[2]]+1]; If[Length[topi]==0,p=0,p=Length[topi]]; q= Length[boti]; rank = p+q; (*vectors*) If[rank==1 && StringQ[exp]==True, vector= Table[ToExpression[exp][ii],{ii,0,DIM-1}]; derv = ToExpression[exp][Switch[Length[ind[[1]]],0,ind[[1]],1,ind[[1]][[1]] ] ] ]; If[rank==1 && Length[exp]>0, vector=exp; derv = exp[[Switch[Length[ind[[1]]],0,ind[[1]]+1,1,ind[[1]][[1]]+1]]] ]; If[rank==1 && q>0, Return[D[derv,COORDINATES[[a+1]]] -Sum[CHARRAY[[dummy+1]][[a+1,ind[[1]]+1]] (vector[[dummy+1]]),{dummy,0,DIM-1}]]]; If[rank==1 && p>0, interind=ind[[1]][[1]]; Return[D[derv,COORDINATES[[a+1]]] +Sum[CHARRAY[[interind+1]][[a+1,dummy+1]] (vector[[dummy+1]]),{dummy,0,DIM-1}]] ]; (*tensors*) indvec=ind; codvp=If[p>0,Table[0,{k,1,p}],{0}]; codvq=If[q>0,Table[0,{k,1,q}],{0}]; If[p>0.1, For[int=1,int
0,
For[int=1,int0,indvec[[2]][[int]]=dummy, indvec[[int]]=dummy];
codvq[[int]] = Evaluate[-(Sum[CHARRAY[[ dummy+1 ]][[a+1,boti[[int]] ]]
(ToExpression[exp][indvec])
,{dummy,0,DIM-1}])];
]];
Return[Sum[codvp[[kk]],{kk,1,If[p>0,p,1]}] +Sum[codvq[[kk]],{kk,1,If[q>0,q,1]}] +
D[ToExpression[exp][ind],COORDINATES[[a+1]]] ]
]
Geodesics[a_]:= Block[{v,int},
If[Length[CHARRAY]==0, Christoffel[{{0},{0,0}}]];
v=Table[0,{int,1,DIM}];
Clear[int];
For[int=1,int < Length[COORDINATES]+1, int++,
v[[int]]=COORDINATES[[int]][a] ];
Table[D[v[[ss]],a,a]
+ Sum[
(CHARRAY[[ss]][[kk,tt]]/.Table[COORDINATES[[yy]]->v[[yy]],{yy,1,Length[COORDINATES]}])
D[v[[kk]],a] D[v[[tt]],a],{kk,1,Length[COORDINATES]},{tt,1,Length[COORDINATES]}],{ss,1,Length[COORDINATES]}]
]
Issol[vi_,s_]:=Block[{v,vec,val,dvec,ddvec,eqs,int,epss},
(*TODO: Make more efficient*)
(*We add a variable in the chanche there is an 1/r infinity*)
If[Length[CHARRAY]==0, Christoffel[{{0},{0,0}}]];
v=vi+epss;
vec=Geodesics[s];
dvec=D[v,s];
ddvec=D[v,s,s];
eqs=Table[10,{kk,1,DIM}];
For[int=1,int