Select Git revision
GRTensor.m 71.24 KiB
(* ::Package:: *)
(************************************************************************)
(* This file was generated automatically by the Mathematica front end. *)
(* It contains Initialization cells from a Notebook file, which *)
(* typically will have the same name as this file except ending in *)
(* ".nb" instead of ".m". *)
(* *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs. Doing so is equivalent *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end. *)
(* *)
(* DO NOT EDIT THIS FILE. This entire file is regenerated *)
(* automatically each time the parent Notebook file is saved in the *)
(* Mathematica front end. Any changes you make to this file will be *)
(* overwritten. *)
(************************************************************************)
(* Probably not all the extra packages are really needed *)
BeginPackage["GRTensor`"];
MetDet::usage="MetDet[g_]. Compute the determinant of the metric";
InverseMetric::usage="InverseMetric[g_]. Compute the inverse of the metric";
ChristoffelSymbol::usage="ChristoffelSymbol[coords_,g_,pert_:0]. Compute Christoffel symbols. Default for perturbation variabel pert is 0.";
WeylTensor::usage="WeylTensor[coords_,g_,pert_:0]. Compute Weyl tensor following the convention of Misner et al., that is, [S2] = 1, [S3] = 1. Default for perturbation variabel pert is 0. ";
RiemannTensor::usage="RiemannTensor[coords_,g_,pert_:0]. Compute Riemann tensor following the convention of Misner et al., that is, [S2] = 1, [S3] = 1 .Default for perturbation variabel pert is 0. ";
RicciTensor::usage="RicciTensor[coords_,g_,pert_:0]. Compute Riemann tensor following the convention of Misner et al., that is, [S2] = 1, [S3] = 1. Default for perturbation variabel pert is 0. ";
RicciScalar::usage="RicciScalar[coords_,g_,pert_:0]. Compute RicciScalar scalar. Default for perturbation variabel pert is 0. ";
KrScalar::usage="KrScalar[coords_,g_]. Compute Kretschmann scalar";
WeylTrace::usage"WeylTrace[coords_,g_]: Compute Weyl Tensor trace";
Einstein::usage="Einstein[coords_,g_,\[Epsilon]p:]. Compute Einstein tensor following the convention of Misner et al., that is, [S2] = 1, [S3] = 1 with perturbation index \[Epsilon]p";
ETensor::usage="ETensor[coords_,g_,\[Epsilon]p_]. Compute the energy momentum tensor with perturbation index \[Epsilon]p.";
DAlembert::usage="DAlembert[coords_,g_,func_]. Compute D'Alembert operator for func[coords]";
CovDer::usage="CovDer[coords_,metric_,tensor_,comps_]. Compute the covariant derivative (default covariant version) for scalar and 1-2 forms. The components are given in a list as: {a},{a,b},{a,b,c}";
NonZeroChristoffel::usage="NonZeroChristoffel[\[CapitalGamma]]. Show the nonzero Christoffel components.";
NonZeroMetricComp::usage="NonZeroMetricComp[g]. Show the nonzero metric components.";
NonZeroTensorComp::usage="NonZeroTensorComp[T]. Show the nonzero Tensor components. It works with any symmetric m xmxmxmx... tensor";
LeviCivitaTensorCurv::usage="LeviCivitaTensorCurv[coords_,g_]. Compute the Levi-Civita antisymmetric tensor for curvilinear coordinates. For cartesian xx recovers the usual \[Epsilon]_(abc).";
CheckTetrad::usage="[gab_,nullv_]. Check whether the 4 null tetrad vectors satisfy orthonormality conditions.";
EinsteinfR::usage="EinsteinfR[coords_,g_,fR_]. Compute fR Einstein equations following the convention of Misner et al., that is, [S2] = 1, [S3] = 1 ";
EinsteinST::usage="EinsteinST[coords_,g_,v\[Phi]_]. Compute scalar-tensor Einstein equations following the convention of Misner et al., that is, [S2] = 1, [S3] = 1 ";
STensorT\[Psi]::usage="STensorT\[Psi][coords_,g_,v\[Phi]_]. Compute scalar-tensor Energy momentum tensor following the convention of Misner et al., that is, [S2] = 1, [S3] = 1 ";
TeffFR::usage="TeffFR[coords_,g_,fR_]. Compute fR Teff tensor such Gab=8\[Pi]/f'[R](Tab + Teff) following the convention of Misner et al., that is, [S2] = 1, [S3] = 1";
TeffST::usage="TeffST[coords_,g_,{V\[CurlyPhi],\[CurlyPhi]}]. Compute ST-EF/JF Teff tensor such Gab=8\[Pi] (Tab + Teff) following the convention of Misner et al., that is, [S2] = 1, [S3] = 1. Allowed options for the Frame are Einstein, Jordan";
FRTOV::usage="FRTOV[coords_,g_,fR_,vars_]. Compute fR TOV eqs such Gab=8\[Pi]/f'[R](Tab + Teff) following the convention of Misner et al., that is, [S2] = 1, [S3] = 1";
STTOV::usage="STTOV[coords_,g_,{V\[CurlyPhi]_,var\[CurlyPhi]_},vars_]. Compute ST-EF/JF TOV eqs such Gab=8\[Pi](Tab + Teff) following the convention of Misner et al., that is, [S2] = 1, [S3] = 1";
fR2Pot::usage="fR2UJF[fR_]. From fR model to the ST potential U(\[Phi])-V(\[Phi])";
fR2\[Phi]::usage="fR2\[Phi][fR_]. From fR model to the ST potential";
CurlCurvilinear::usage"CurlCurvilinear[xx,g,vec]. It computes the curl tensor in curvilinear coordinates";
ElectricTensor3p1Dev::usage="ElectricTensor3p1Dev[xx_,g_,Kt_,pert_:0,OptionsPattern[]]";
MagneticTensor3p1Dev::usage="MagneticTensor3p1Dev[xx_,g_,Kt_,pert_:0,OptionsPattern[]]"
EoSCallParsBSks::usage="EoSCallParsBSks[model_]: Parameters for BSK1-3 EOS for the matter density.";
EoSCallParsSly::usage="EoSCallParsSly[model_]: Parameters for SLy1 EOS for the matter density.";
EoSBSks::usage="EoSBSks[model_]: Analytic EOS for BSK1-3for BSK1-3 EOS for the matter density.";
EoSSly::usage="EoSSly[model_]: Analytic EOS for SLy1 EOS for the matter density.";
SlyInner::usage="";
SlyLCore::usage="";SlyLCoreAll::usage="";
SlyInnerAll::usage="";
EoSFitsPars::usage="EoSFitsPars[model_,verbose_:False]. Parameters for the JRead(arxiv:0812.2163) parameters.";
EoSSlyCrust::usage="EoSSlyCrust[\[Rho]_]. Crust model for the NS of JRead(arxiv:0812.2163)";
EoSFits::usage="EoSFits[model_,OptionsPattern[]]. NS EOS of JRead(arxiv:0812.2163)";
EoSPol::usage="EoSPol[model_]. NS EOS for a non-relativistic (NR) and relativistic (R, default) NS respectively";
EoSPol\[Epsilon]::usage="EoSPol\[Epsilon][model_]. NS EOS for a non-relativistic (NR) and relativistic (R, default) NS respectively";
From\[Rho]To\[Epsilon]Fits::usage"From\[Rho]To\[Epsilon]Fits[eos_]. NS EOS of JRead(arxiv:0812.2163) for the energy density.";
\[Omega]eff::usage="\[Omega]eff[fR_,OptionsPattern]. It computes the effective equation of state given the fR model. It also provides the same for scalar tensor theory in the JF.";
ShootingNStars::usage="ShootingNStars[eqs_,eqsRg_,rvar_,vars_,shtdInd_,varshtdRg_]. Shooting function of the index var shtdInd for a set of eqs integrated in eqsRg on the variables vars ;";
BracketingSTNStars::usage="BracketingSTNStars[eqs_,eqsRg_,rvar_,vars_,shtdInd_,varshtdRg_]. Shooting function of the index var shtdInd for a set of eqs integrated in eqsRg on the variables vars ;";
RK4::usage="RK4[func_?ListQ,vars_?ListQ,ivals_?ListQ,pars_?ListQ,step_]";
TestCode::usage"TestCode[eqs_,vars_,icond_,rlst_,drlst_]";
AtomsList::usage="Take the coefficients out";
InterpolationDomain::usage="InterpolationDomain[interpolatedfunction]. It outputs the domain in a format {tmin,tmax}";
TakeColumn::usage="TakeColumn[list1_,list2_] extracts columns list2 from list1, i.e. it gives the functionality of the pre-Mathematica 6 Column function.";
(* All here below are dev. versions. Replace them when they are sufficiently tested. *)
ChristoffelSymbolDev::usage="ChristoffelSymbolDev[coords_,g_,pert_:0]. Compute Christoffel symbols. Default for perturbation variabel pert is 0."
RiemannTensorDev::usage="RiemannTensor[coords_,g_,pert_:0]. Compute Riemann tensor following the convention of Misner et al., that is, [S2] = 1, [S3] = 1 .Default for perturbation variabel pert is 0. "
RicciTensorDev::usage="RicciTensorDev[coords_,g_,pert_:0]. Compute Riemann tensor following the convention of Misner et al., that is, [S2] = 1, [S3] = 1 .Default for perturbation variabel pert is 0. ";
Begin["`Private`"];
Options[MetDet]={"PerturbationIndex"->1,"SimplifyFunction"->Identity};
MetDet[g_,OptionsPattern[]]:=Block[{simpl},simpl=OptionValue["SimplifyFunction"]; simpl@Det[g]];
Options[InverseMetric]={"PerturbationIndex"->1,"SimplifyFunction"->Identity};
InverseMetric[g_,OptionsPattern[]]:=Block[{simpl},simpl=OptionValue["SimplifyFunction"]; simpl@Inverse[g]]
Options[ChristoffelSymbol]={"Verbose"->False,"PerturbationIndex"->1,"SimplifyFunction"->Identity};
ChristoffelSymbol[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{n,ig,res,perti,simpl,verbose},
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
verbose=OptionValue["Verbose"];
n=Length@xx;
ig=InverseMetric[g];
res=Table[(1/2)*If[NumericQ[pert],Sum[ig[[i,s]]*(-D[g[[j,k]],xx[[s]]]+D[g[[j,s]],xx[[k]]]+D[g[[s,k]],xx[[j]]]),{s,1,n}],
Normal@Series[Sum[ig[[i,s]]*(-D[g[[j,k]],xx[[s]]]+D[g[[j,s]],xx[[k]]]+D[g[[s,k]],xx[[j]]]),{s,1,n}],{pert,0,perti}]],{i,1,n},{j,1,n},{k,1,n}];
simpl@res]
Options[RiemannTensor]=Options[ChristoffelSymbol];
RiemannTensor[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{n,Chr,res,perti,simpl,verbose},
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
verbose=OptionValue["Verbose"];
n=Length@xx;
If[verbose,Print["Starting with Christoffel symbols..."]];
Chr=ChristoffelSymbol[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
If[verbose,Print["Christoffel symbols computed. Starting with Riemann..."]];
res=Table[If[NumericQ[pert],D[Chr[[i,k,m]],xx[[l]]]-D[Chr[[i,k,l]],xx[[m]]]+Sum[Chr[[i,s,l]]*Chr[[s,k,m]],{s,1,n}]-Sum[Chr[[i,s,m]]*Chr[[s,k,l]],{s,1,n}],
Normal@Series[D[Chr[[i,k,m]],xx[[l]]]-D[Chr[[i,k,l]],xx[[m]]]+Sum[Chr[[i,s,l]]*Chr[[s,k,m]],{s,1,n}]-Sum[Chr[[i,s,m]]*Chr[[s,k,l]],{s,1,n}],{pert,0,perti}]],{i,1,n},{k,1,n},{l,1,n},{m,1,n}];
If[verbose,Print["...Riemann computed"]];
simpl@res];
Options[WeylTensor]=Options[ChristoffelSymbol];WeylTensor[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{n,Chr,riemann,riemanndown,ricciS,ricciT,res,perti,simpl,verbose},
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
verbose=OptionValue["Verbose"];
n=Length@xx;
If[verbose,Print["Starting with RicciScalar..."]];
ricciS=RicciScalar[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
If[verbose,Print["Following with RicciTensor..."]];
ricciT=RicciTensor[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
If[verbose,Print["Following with Riemann..."]];
riemann=RiemannTensor[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
riemanndown=Table[Sum[g[[a,\[Alpha]]]riemann[[\[Alpha],b,c,d]],{\[Alpha],4}],{a,n},{b,n},{c,n},{d,n}];
res=Table[If[NumericQ[pert],riemanndown[[i,k,l,m]]+1/(n-2)(ricciT[[i,m]]g[[k,l]]-ricciT[[i,l]]g[[k,m]]+ricciT[[k,l]]g[[i,m]]-ricciT[[k,m]]g[[i,l]])+1/((n-1)(n-2))ricciS(g[[i,l]]g[[k,m]]-g[[i,m]]g[[k,l]]),
Normal@Series[riemanndown[[i,k,l,m]]+1/(n-2)(ricciT[[i,m]]g[[k,l]]-ricciT[[i,l]]g[[k,m]]+ricciT[[k,l]]g[[i,m]]-ricciT[[k,m]]g[[i,l]])+1/((n-1)(n-2))ricciS(g[[i,l]]g[[k,m]]-g[[i,m]]g[[k,l]]),{pert,0,perti}]],{i,1,n},{k,1,n},{l,1,n},{m,1,n}];
(*Simplify[res]*)
simpl@res];
Options[RicciTensor]=Options[ChristoffelSymbol];
RicciTensor[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{Rie,res,n,perti,simpl},
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
n=Length@xx;
Rie=RiemannTensor[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
res=Table[If[NumericQ[pert],Sum[Rie[[s,i,s,j]],{s,1,n}],
Normal@Series[Sum[Rie[[s,i,s,j]],{s,1,n}],{pert,0,perti}]],{i,1,n},{j,1,n}];
(*Simplify[res]*)
simpl@res]
Options[RicciScalar]=Options[ChristoffelSymbol];
RicciScalar[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{Ricc,ig,res,n,perti,simpl},
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
n=Length@xx;
Ricc=RicciTensor[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
ig=InverseMetric[g,"SimplifyFunction"->simpl];
res=If[NumericQ[pert],Sum[ig[[s,i]] Ricc[[s,i]],{s,1,n},{i,1,n}],
Normal@Series[Sum[ig[[s,i]] Ricc[[s,i]],{s,1,n},{i,1,n}],{pert,0,perti}]];
simpl@res]
Options[KrScalar]=Options[ChristoffelSymbol];
KrScalar[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{Rie,res,n,Ried,Rieup,gup,perti,simpl},
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
n=Length@xx;
Rie=RiemannTensor[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
gup=Inverse@g;
Ried=Table[Sum[g[[i,k]] Rie[[i,j,m,l]],{i,1,n}],{k,1,n},{j,1,n},{m,1,n},{l,1,n}];
Rieup=Table[Sum[gup[[j,k]]gup[[m,\[Mu]]] gup[[l,\[Nu]]] Rie[[i,j,m,l]],{j,1,n},{m,1,n},{l,1,n}],{i,1,n},{k,1,n},{\[Mu],1,n},{\[Nu],1,n}];
res=Sum[Ried[[i,j,m,l]] Rieup[[i,j,m,l]],{i,1,n},{j,1,n},{m,1,n},{l,1,n}];
If[NumericQ[pert],simpl@res,Normal@Series[simpl@res,{pert,0,perti}]]];
Options[WeylTrace]=Options[ChristoffelSymbol];
WeylTrace[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{Chr,ig,n,riemann,Rieup,ricciS,ricciT,res,perti,simpl,weylT,weylTup},
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
n=Length@xx;
ig=InverseMetric[g,"SimplifyFunction"->simpl];
weylT=WeylTensor[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
weylTup=Table[Sum[ig[[a,\[Alpha]]] ig[[b,\[Beta]]] ig[[c,\[Gamma]]] ig[[d,\[Delta]]] weylT[[\[Alpha],\[Beta],\[Gamma],\[Delta]]],{\[Alpha],n},{\[Beta],n},{\[Gamma],n},{\[Delta],n}],{a,n},{b,n},{c,n},{d,n}];
res=If[NumericQ[pert],Sum[weylT[[i,k,l,m]]weylTup[[i,k,l,m]],{i,n},{k,n},{l,n},{m,n}],
Normal@Series[Sum[weylT[[i,k,l,m]]weylTup[[i,k,l,m]],{i,n},{k,n},{l,n},{m,n}],{pert,0,perti}]];
simpl@res];
Options[Einstein]=Options[ChristoffelSymbol];
Einstein[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{res,perti,simpl},
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
res=RicciTensor[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl]-1/2 RicciScalar[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl] g;
If[NumericQ[pert],simpl@res,Normal@Series[simpl@res,{pert,0,perti}]]
]
DAlembert[xx_,metric_,scalfun_]:=Block[{metdet,dal,metup,sqrtdet},
metup=Inverse[metric];
metdet=MetDet[metric];
sqrtdet=Sqrt[-metdet];
dal=FullSimplify[(1/sqrtdet) Sum[
D[sqrtdet metup[[\[Nu],\[Mu]]] D[scalfun,xx[[\[Mu]]]],xx[[\[Nu]]]],{\[Nu],4},{\[Mu],4}]]
]
Options[ETensor]=Join[{"Signature"->1},Options[ChristoffelSymbol]];
ETensor[coors_,met_,pert_:0,OptionsPattern[]]:=Block[{g,gup,Global`\[Rho],\[Rho]c,Global`p,perti,pc,riscal,riscalvars,sign,simpl,T\[Mu]\[Nu],T\[Mu]\[Nu]up,u,udown,\[CapitalOmega]},
perti=OptionValue["PerturbationIndex"];
sign=OptionValue["Signature"];
simpl=OptionValue["SimplifyFunction"];
g=met;
riscal=RicciScalar[coors,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
riscalvars=Complement[coors,Complement[coors,AtomsList[riscal]]];
gup=Inverse[g];
pc=Global`p@@riscalvars;
\[Rho]c=Global`\[Rho]@@riscalvars;
u=If[sign==1,{1/Sqrt[-g[[1,1]]],0,0,pert \[CapitalOmega]/Sqrt[-g[[1,1]]]},{1/Sqrt[g[[1,1]]],0,0,pert \[CapitalOmega]/Sqrt[g[[1,1]]]}];
u=Normal@Series[u,{pert,0,perti}];
udown=Table[Sum[u[[i]]g[[i,l]],{i,Length@coors}],{l,Length@coors}];
udown=Normal@Series[udown,{pert,0,perti}];
If[Simplify[Normal@Series[u.udown,{pert,0,1}]]!=-1,Return["Wrong normalization"]];
T\[Mu]\[Nu]up=If[sign==1,Normal@Series[Table[(pc+\[Rho]c) u[[i]] u[[j]]+pc gup[[i,j]],{i,Length@coors},{j,Length@coors}],{pert,0,perti}],
Normal@Series[Table[(pc+\[Rho]c) u[[i]] u[[j]]-pc gup[[i,j]],{i,Length@coors},{j,Length@coors}],{pert,0,perti}]]//Simplify;
T\[Mu]\[Nu]=Normal@Series[g(T\[Mu]\[Nu]up.g),{pert,0,1}]
]
Options[NonZeroTensorComp]={"Verbose"->True,"TensorString"->"T"};
NonZeroTensorComp[gamma_,OptionsPattern[]]:=Module[{dimension,ii,nonzerpos,tstring,verbose,
allpos,zerpos},
tstring=OptionValue["TensorString"];
verbose=OptionValue["Verbose"];
(* set dimensions and total number or elements *)
dimension=Dimensions[gamma];
(*allpos=Flatten[Table[{i,j,k},{i,dimension[[1]]},{j,dimension[[1]]},{k,dimension[[1]]}],1];*)
allpos=Tuples[Table[i,{i,dimension[[1]]}],Length@dimension];
zerpos=Position[gamma,_?(#== 0 &)];
nonzerpos=Complement[allpos,zerpos];
If[verbose,Do[Print[ToString[Subscript[tstring, StringJoin[ToString/@(nonzerpos[[i,1;;Length@dimension]]-1)]],StandardForm]<>" = "<>ToString[Extract[gamma,nonzerpos[[i]]],StandardForm]],{i,Length@nonzerpos}]];
nonzerpos
]
Options[NonZeroMetricComp]={"Verbose"->True,"TensorString"->"T"};NonZeroMetricComp[gamma_,OptionsPattern[]]:=Module[{allpos,dimension,
n,nonzerpos,tstring,verbose,zerpos},
verbose=OptionValue["Verbose"];
tstring=OptionValue["TensorString"];
n=Length@gamma;
(* set dimensions and total number or elements *)
dimension=Dimensions[gamma][[1]];
allpos=Flatten[Table[{i,j},{i,dimension},{j,dimension}],1];
zerpos=Position[gamma,_?(#== 0 &)];
nonzerpos=Complement[allpos,zerpos];
If[verbose,Do[Print[ToString[tstring,StandardForm]<>" = "<>ToString[gamma[[nonzerpos[[i,1]],nonzerpos[[i,2]]]],StandardForm]],{i,Length@nonzerpos}]];
nonzerpos
]
CovDer[coords_,metric_,tensor_,index_,OptionsPattern[{"Valence"->"Covariant","Verbose"->False}]]:=Block[{Crhistoffel,g,xx,n,h\[Eta]\[Nu],rank,cov,valence,verbose,c,b,a},
(*Print[Style["It is wrong!",Red]];*)
valence=OptionValue["Valence"];
verbose=OptionValue["Verbose"];
n=Length@coords;
g=metric;
h\[Eta]\[Nu]=tensor;
xx=coords;
If[ListQ[h\[Eta]\[Nu]],
rank=Length@Dimensions@tensor;
Crhistoffel=ChristoffelSymbol[xx,g],
rank=0;
];
Which[rank==0, If[verbose,Print["Scalar , "];Print[valence];];
{a}=index;
cov=D[h\[Eta]\[Nu],xx[[index]]]; ,
rank==1, If[verbose,Print["Vector , "];Print[valence];];
{a,b}=index;
Which[valence=="Covariant",
cov=D[h\[Eta]\[Nu][[b]],xx[[a]]]-Sum[Crhistoffel[[\[Rho],a,b]]h\[Eta]\[Nu][[\[Rho]]],{\[Rho],Length@xx}];,
valence=="Contravariant",
cov=D[h\[Eta]\[Nu][[b]],xx[[a]]]+Sum[Crhistoffel[[\[Rho],a,b]]h\[Eta]\[Nu][[\[Rho]]],{\[Rho],Length@xx}];
];,
rank==2, If[verbose,Print["Tensor , "];Print[valence];];
{a,b,c}=index;
Which[valence=="Covariant",
cov=D[h\[Eta]\[Nu][[b,c]],xx[[a]]]-Sum[Crhistoffel[[d,a,b]]h\[Eta]\[Nu][[c,d]]+Crhistoffel[[d,a,c]]h\[Eta]\[Nu][[b,d]],{d,n}];,
valence=="Contravariant",
cov=D[h\[Eta]\[Nu][[b,c]],xx[[a]]]+Sum[Crhistoffel[[b,a,d]]h\[Eta]\[Nu][[c,d]]+Crhistoffel[[c,a,d]]h\[Eta]\[Nu][[d,b]],{d,n}];,
valence=="Mixed",
cov=D[h\[Eta]\[Nu][[b,c]],xx[[a]]]+Sum[Crhistoffel[[b,a,d]]h\[Eta]\[Nu][[c,d]]-Crhistoffel[[c,a,d]]h\[Eta]\[Nu][[d,b]],{d,n}];];,
True, If[verbose,Print["Rank not recognised"]];
Return[];
];
Simplify@cov
]
LeviCivitaTensorCurv[xx_,g_]:=Module[{dim},
dim=Length@xx;
Sqrt[Det[g]]LeviCivitaTensor[dim,List]
]
CurlCurvilinear[xx_,g_,vec_]:=Module[{lv,vecb,DCov,lvup,gup,vecn},
lv=LeviCivitaTensorCurv[xx,g];
vecb=Sqrt[Diagonal[g]];
vecn=vec*vecb;
DCov=Table[CovDer[xx,g,vecn,{a,b}],{a,Length@xx},{b,Length@xx}];gup=InverseMetric[g];
lvup=Simplify[Table[Sum[gup[[i,a]]gup[[j,b]] gup[[k,c]] lv[[a,b,c]],{a,Length@xx},{b,Length@xx},{c,Length@xx}],{i,Length@xx},{j,Length@xx},{k,Length@xx}]];
Table[Simplify[Sum[lvup[[i,e,d]]DCov[[e,d]],{e,Length@xx},{d,Length@xx}]]vecb[[i]],{i,Length@xx}]
]
Options[CheckTetrad]=Options[ChristoffelSymbol];
CheckTetrad[gab_,nullv_,OptionsPattern[]]:=Block[{l,n,m,mb,test,verbose},
{l,n,m,mb}=nullv;
verbose=OptionValue["Verbose"];
test=Chop@Simplify@{l.gab.l,
n.gab.n,
m.gab.m,
mb.gab.mb,
l.gab.m,
l.gab.mb,
n.gab.m,
n.gab.mb,
l.gab.n,
m.gab.mb};
If[verbose,Print["{l.l,n.n,m.m,mb.mb,l.m,l.mb,n.m,n.mb,l.n,m.mb} = ",test]];
If[Total@test[[1;;8]]==0&&test[[9]]*test[[10]]==-1,Print[Style["Your tetrad satisfies the properties of orthonormality",Blue]];,Print[Style["Your tetrad has some troubles",Blue]]]];
Options[ChristoffelSymbolDev]={"Verbose"->False,"PerturbationIndex"->1,"SimplifyFunction"->Identity,"Compile"->False,"CompileCoordinateIndex"->{1,2,3,4}};
ChristoffelSymbolDev[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{n,ig,compile,compilecind,coords,res,perti,simpl,verbose},
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
verbose=OptionValue["Verbose"];
compile=OptionValue["Compile"];
compilecind=OptionValue["CompileCoordinateIndex"];
n=Length@xx;
ig=InverseMetric[g];
res=ConstantArray[0,{n,n,n}];
If[NumericQ[pert], Do[res[[i,j,k]]=Sum[ig[[i,s]]*(-D[g[[j,k]],xx[[s]]]+D[g[[j,s]],xx[[k]]]+D[g[[s,k]],xx[[j]]]),{s,n}],{i,n},{j,n},{k,j,n}];,
Do[res[[i,j,n]]=Normal@Series[Sum[ig[[i,s]]*(-D[g[[j,k]],xx[[s]]]+D[g[[j,s]],xx[[k]]]+D[g[[s,k]],xx[[j]]]),{s,n}],{pert,0,perti}],{i,n},{j,n},{k,j,n}]];
res=simpl[1/2 res];
(* Compile *)
If[compile, Do[res[[i,j,k]]=If[NumberQ[res[[i,j,k]]],res[[i,j,k]],Compile[Evaluate@({#,_Real}&/@(xx[[compilecind]])),Evaluate[res[[i,j,k]]],CompilationTarget->"C",CompilationOptions->"InlineExternalDefinitions"->True]],{i,n},{j,n},{k,j,n}];];
(* Applying symmetries *)
Do[res[[i,j+1,k]]=res[[i,k,j+1]];,{i,n},{j,n-1},{k,j}];
res
]
Options[RiemannTensorDev]=Join[Options[ChristoffelSymbolDev],{"IndexDown"->False}];
RiemannTensorDev[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{caux,n,Chr,compile,compilecind,index,res,perti,simpl,time,verbose},
index=OptionValue["IndexDown"];
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
compile=OptionValue["Compile"];
compilecind=OptionValue["CompileCoordinateIndex"];
verbose=OptionValue["Verbose"];
n=Length@xx;
If[verbose,Print["Starting with Christoffel symbols..."]];
Chr=ChristoffelSymbolDev[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl,"Compile"->False];
If[verbose,Print["Christoffel symbols computed. Starting with Riemann..."]];
res=ConstantArray[0,{n,n,n,n}];
If[index,
If[NumericQ[pert], Do[res[[i,k,l,m]]=Sum[g[[i,p]](D[Chr[[p,k,m]],xx[[l]]]-D[Chr[[p,k,l]],xx[[m]]]+Sum[Chr[[p,s,l]]*Chr[[s,k,m]]-Chr[[p,s,m]]*Chr[[s,k,l]],{s,n}]),{p,n}],{i,n},{k,i,n},{l,n},{m,l,n}], Do[res[[i,k,l,m]]=Sum[g[[i,p]](Normal@Series[D[Chr[[p,k,m]],xx[[l]]]-D[Chr[[p,k,l]],xx[[m]]]+Sum[Chr[[p,s,l]]*Chr[[s,k,m]],{s,n}]-Sum[Chr[[p,s,m]]*Chr[[s,k,l]],{s,n}],{pert,0,perti}]),{p,n}],{i,n},{k,n},{l,n},{m,n}]];
(* Applying simmetries *)
(*Do[res[[i+1,k,l,m]]=res[[l,k,i+1,m]],{i,n-1},{k,n},{l,i},{m,n}];*)
Do[res[[i,k,l+1,m]]=-res[[i,k,m,l+1]],{i,n},{k,n},{l,n-1},{m,l}];
Do[res[[i+1,k,l,m]]=-res[[k,i+1,l,m]],{i,n-1},{k,i},{l,n},{m,n}];
If[compile,caux=0; Do[time=Timing[res[[i,k,l,m]]=Compile[Evaluate@({#,_Real}&/@(xx[[compilecind]])),Evaluate[res[[i,k,l,m]]],RuntimeOptions->"Speed"]][[1]];If[verbose,Print["compiling: ",{i,k,l,m}," ",time," s"]];,{i,n},{k,n},{l,n},{m,n}];]
,
If[NumericQ[pert], Do[res[[i,k,l,m]]=(D[Chr[[i,k,m]],xx[[l]]]-D[Chr[[i,k,l]],xx[[m]]]+Sum[Chr[[i,s,l]]*Chr[[s,k,m]]-Chr[[i,s,m]]*Chr[[s,k,l]],{s,n}]),{i,n},{k,n},{l,n},{m,n}],
Do[res[[i,k,l,m]]=(Normal@Series[D[Chr[[i,k,m]],xx[[l]]]-D[Chr[[i,k,l]],xx[[m]]]+Sum[Chr[[i,s,l]]*Chr[[s,k,m]],{s,n}]-Sum[Chr[[i,s,m]]*Chr[[s,k,l]],{s,n}],{pert,0,perti}]),{i,n},{k,n},{l,n},{m,n}]];
];
If[verbose,Print["...Riemann computed"]];
simpl@res];
Options[RicciTensorDev]=Options[ChristoffelSymbolDev];
RicciTensorDev[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{compile,Rie,res,n,perti,simpl},
perti=OptionValue["PerturbationIndex"];
compile=OptionValue["Compile"];
simpl=OptionValue["SimplifyFunction"];
n=Length@xx;
Rie=RiemannTensorDev[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
res=ConstantArray[0,{n,n}];
If[NumericQ[pert], Do[res[[i,j]]=Sum[Rie[[s,i,s,j]],{s,n}],{i,n},{j,n}],
Do[res[[i,j]]=Normal@Series[Sum[Rie[[s,i,s,j]],{s,n}],{pert,0,perti}],{i,n},{j,i,n}]];
If[compile, Do[res[[i,j]]=If[NumberQ[res[[i,j]]],res[[i,j]],Compile[Evaluate@({#,_Real}&/@xx),Evaluate[res[[i,j]]],CompilationTarget->"C",CompilationOptions->"InlineExternalDefinitions"->True]],{i,n},{j,i,n}];] ;
(* Applying symmetries *)
Do[res[[i+1,j]]=res[[j,i+1]];,{i,n-1},{j,i}];
simpl@res]
Options[ElectricTensor3p1Dev]=Options[ChristoffelSymbolDev];
ElectricTensor3p1Dev[xx_,g_,Kt_,pert_:0,OptionsPattern[]]:=Block[{compile,gup,Ks,Ktup,Ric,res,n,perti,simpl},
perti=OptionValue["PerturbationIndex"];
compile=OptionValue["Compile"];
simpl=OptionValue["SimplifyFunction"];
n=Length@xx;
Ric=RicciTensorDev[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
gup=InverseMetric[g];
Ktup=(gup.gup.Kt);
Ks=(Sum[gup[[a,b]]Kt[[a,b]],{a,n},{b,n}]);
res=ConstantArray[0,{n,n}];
If[NumericQ[pert], Do[res[[i,j]]=Ric[[i,j]]+Ks Kt[[i,j]]-Sum[Kt[[i,m]]Ktup[[m,j]],{m,n}],{i,n},{j,i,n}],
Do[res[[i,j]]=Normal@Series[Ric[[i,j]]+Ks Kt[[i,j]]-Sum[Kt[[i,m]]Ktup[[m,j]],{m,n}],{pert,0,perti}],{i,n},{j,i,n}]];
If[compile, Do[res[[i,j]]=If[NumberQ[res[[i,j]]],res[[i,j]],Compile[Evaluate@({#,_Real}&/@xx),Evaluate[res[[i,j]]],CompilationTarget->"C",CompilationOptions->"InlineExternalDefinitions"->True]],{i,n},{j,i,n}];];
(* Applying symmetries *)
Do[res[[i+1,j]]=res[[j,i+1]];,{i,n-1},{j,i}];
simpl@res]
Options[MagneticTensor3p1Dev]=Options[ChristoffelSymbolDev];
MagneticTensor3p1Dev[xx_,g_,Kt_,pert_:0,OptionsPattern[]]:=Block[{compile,gup,Ks,Ktup,Ric,res,rescov,n,perti,simpl,\[Epsilon]},
perti=OptionValue["PerturbationIndex"];
compile=OptionValue["Compile"];
simpl=OptionValue["SimplifyFunction"];
n=Length@xx;
Ric=RicciTensorDev[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
gup=InverseMetric[g];
\[Epsilon]=Det[gup]^(-1/2)LeviCivitaTensor[n];\[Epsilon]=Table[Sum[g[[s,i]]\[Epsilon][[s,j,k]],{s,n}],{i,n},{j,n},{k,n}];
rescov=ConstantArray[0,{n,n,n}];
Do[rescov[[k,i,j]]=CovDer[xx,g,Kt,{k,i,j},"Valence"->"Covariant"],{k,n},{i,n},{j,i,n}];
Do[rescov[[k,i+1,j]]=rescov[[k,j,i+1]];,{k,n},{i,n-1},{j,i}];
res=ConstantArray[0,{n,n}];
If[NumericQ[pert], Do[res[[i,j]]=Sum[\[Epsilon][[i,m,k]]rescov[[m,k,j]],{m,n},{k,n}],{i,n},{j,i,n}],
Do[res[[i,j]]=Normal@Series[Sum[\[Epsilon][[i,m,k]]rescov[[m,k,j]],{m,n},{k,n}],{pert,0,perti}],{i,n},{j,i,n}]];
If[compile, Do[res[[i,j]]=If[NumberQ[res[[i,j]]],res[[i,j]],Compile[Evaluate@({#,_Real}&/@xx),Evaluate[res[[i,j]]],CompilationTarget->"C",CompilationOptions->"InlineExternalDefinitions"->True]],{i,n},{j,i,n}];];
(* Applying symmetries *)
Do[res[[i+1,j]]=res[[j,i+1]];,{i,n-1},{j,i}];
simpl@res]
Options[RicciScalarDev]=Options[ChristoffelSymbolDev];
RicciScalarDev[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{Ricc,ig,res,n,perti,simpl},
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
n=Length@xx;
Ricc=RicciTensorDev[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
ig=InverseMetric[g,"SimplifyFunction"->simpl];
res=If[NumericQ[pert],Sum[If[i==s,ig[[s,i]] Ricc[[s,i]],2 ig[[s,i]] Ricc[[s,i]]],{s,n},{i,s}],
Normal@Series[Sum[ig[[s,i]] Ricc[[s,i]],{s,1,n},{i,1,n}],{pert,0,perti}]];
simpl@res]
Options[WeylTensorDev]=Options[ChristoffelSymbolDev];
WeylTensorDev[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{n,Chr,riemann,riemanndown,ricciS,ricciT,res,perti,simpl,verbose},
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
verbose=OptionValue["Verbose"];
n=Length@xx;
If[verbose,Print["Starting with RicciScalar..."]];
ricciS=RicciScalarDev[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
If[verbose,Print["Following with RicciTensor..."]];
ricciT=RicciTensorDev[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
If[verbose,Print["Following with Riemann..."]];
riemann=RiemannTensorDev[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl,"IndexDown"->False];
res=Table[If[NumericQ[pert],riemann[[i,k,l,m]]+1/(n-2)(ricciT[[i,m]]g[[k,l]]-ricciT[[i,l]]g[[k,m]]+ricciT[[k,l]]g[[i,m]]-ricciT[[k,m]]g[[i,l]])+1/((n-1)(n-2))ricciS(g[[i,l]]g[[k,m]]-g[[i,m]]g[[k,l]]),
Normal@Series[riemann[[i,k,l,m]]+1/(n-2)(ricciT[[i,m]]g[[k,l]]-ricciT[[i,l]]g[[k,m]]+ricciT[[k,l]]g[[i,m]]-ricciT[[k,m]]g[[i,l]])+1/((n-1)(n-2))ricciS(g[[i,l]]g[[k,m]]-g[[i,m]]g[[k,l]]),{pert,0,perti}]],{i,1,n},{k,1,n},{l,1,n},{m,1,n}];
(*Simplify[res]*)
simpl@res];
Options[\[Omega]eff]={"ST"->False,"SimplifyFunction"->Identity};
\[Omega]eff[fR_,OptionsPattern[]]:=Block[{dR,dR2,expr,Global`x,f,fRterm1,fRc,fRterm2,dalem,dfRc,dfRc2,dfRc3,Global`H,Hc,Global`R,Rc,Rc\[Phi],res,riscal,riscalvars,simpl,st,Global`t,tc,Global`\[Phi],\[Phi]c},
st=OptionValue["ST"];
simpl=OptionValue["SimplifyFunction"];
Rc=Global`R;
tc=Global`t;
Hc=Global`H;
\[Phi]c=Global`\[Phi];
If[st,Rc\[Phi]=fR2\[Phi][fR]/.\[Phi]c->\[Phi]c[t]];
dR=D[Rc[tc],tc];
dR2=D[Rc[tc],{tc,2}];
fRc=fR/.Global`R->Rc[tc];
dfRc=D[fRc,Rc[tc]];
dfRc2=D[fRc,{Rc[tc],2}];
dfRc3=D[fRc,{Rc[tc],3}];
If[fRc==Rc[tc],Return[-1]];
res=(dR^2 dfRc3 + 2Hc[tc]dfRc2 dR + dR2 dfRc2 + 1/2(fRc-Rc[tc]dfRc))/((Rc[tc]dfRc-fRc)/2-3Hc[tc]dR dfRc2);
If[st,res=Block[{R},Rc[t_]:=Evaluate@Rc\[Phi];res]];
simpl@res
];
Options[EinsteinfR]={{"Metric"->True},Join[Options[ChristoffelSymbol]]};
EinsteinfR[xx_,g_,fR_,OptionsPattern[]]:=Block[{res,fRterm1,fRterm2,dalem,Global`R,riscal,riscalvars,simpl,Rc,fRc,dfRc,covterm1,covterm2,metric},
metric=OptionValue["Metric"];
simpl=OptionValue["SimplifyFunction"];
riscal=simpl@RicciScalar[xx,g];
riscalvars=Complement[xx,Complement[xx,AtomsList[riscal]]];
Rc=Global`R@@riscalvars;
fRc=fR/.Global`R->Rc;
dfRc=D[fRc,Rc];
dalem=DAlembert[xx,g,dfRc];
fRterm1=Simplify[dfRc RicciTensor[xx,g]-1/2 g fRc];
If[metric,
covterm1=Table[CovDer[xx,g,dfRc,{i},"Valence"->"Covariant","Verbose"->False],{i,Length@xx}];
covterm2=Table[CovDer[xx,g,covterm1,{i,j},"Valence"->"Covariant","Verbose"->False],{i,Length@xx},{j,Length@xx}];
fRterm2=covterm2-g*dalem;
res=(fRterm1-fRterm2);
Return[{covterm1}];
,
res=(fRterm1);
];
Simplify@res
]
STensorT\[Psi][coor_,met_,\[Psi]potl_]:=Module[{a,b,l,m,g,xx,\[Psi],pot,gup},
xx=coor;
g=met;
pot=\[Psi]potl[[2]];
\[Psi]=\[Psi]potl[[1]];
gup=Inverse@g;
Table[D[\[Psi],xx[[a]]]D[\[Psi],xx[[b]]]-1/2g[[a,b]]Sum[gup[[l,m]]D[\[Psi],xx[[l]]]D[\[Psi],xx[[m]]],{l,4},{m,4}]-g[[a,b]] pot ,{a,Length@xx},{b,Length@xx}]
]
Options[TeffFR]=Join[{"Metric"->True},Options[ChristoffelSymbol]];
TeffFR[xx_,g_,fR_,pert_:0,OptionsPattern[]]:=Block[{res,fRterm1,fRterm2,dalem,Global`R,riscal,riscalvars,Rc,fRc,dfRc,simpl,covterm1,covterm2,metric,perti},
metric=OptionValue["Metric"];
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
riscal=simpl@RicciScalar[xx,g,pert,"PerturbationIndex"->perti];
riscalvars=Complement[xx,Complement[xx,AtomsList[riscal]]];
Rc=Global`R@@riscalvars;
fRc=fR/.Global`R->Rc;
dfRc=D[fRc,Rc];
dalem=DAlembert[xx,g,dfRc];
If[metric,
covterm1=Table[CovDer[xx,g,dfRc,{i},"Valence"->"Covariant","Verbose"->False],{i,Length@xx}];
covterm2=Table[CovDer[xx,g,covterm1,{i,j},"Valence"->"Covariant","Verbose"->False],{i,Length@xx},{j,Length@xx}];
fRterm2=covterm2-g*dalem;
res=If[Not@NumericQ@pert,Normal@Series[1/(8\[Pi])(fRterm2-dfRc/2 Rc g + 1/2 fRc g),{pert,0,perti}],1/(8\[Pi])(fRterm2-dfRc/2 Rc g + 1/2 fRc g)];
,
res=If[Not@NumericQ@pert,Normal@Series[(fRterm1),{pert,0,perti}],(fRterm1)];
];
Simplify@res
]
Options[TeffST]=Join[{"Frame"->"Einstein"},Options[ChristoffelSymbol]];
TeffST[xx_,g_,{V\[CurlyPhi]_,var\[CurlyPhi]_},pert_:0,OptionsPattern[]]:=Block[{covterm1,covterm2,dalem,der,gup,frame,perti,res,riscal,riscalvars,simpl,sumder,V\[CurlyPhi]c,\[CurlyPhi]c,\[Phi]term2},
perti=OptionValue["PerturbationIndex"];
frame=OptionValue["Frame"];
simpl=OptionValue["SimplifyFunction"];
riscal=simpl@RicciScalar[xx,g,pert,"PerturbationIndex"->perti];
riscalvars=Complement[xx,Complement[xx,AtomsList[riscal]]];
gup=Inverse[g];
\[CurlyPhi]c=var\[CurlyPhi]@@riscalvars;
V\[CurlyPhi]c=V\[CurlyPhi]/.var\[CurlyPhi]->\[CurlyPhi]c;
Which[frame=="Einstein",
der=Table[D[\[CurlyPhi]c,xx[[a]]]D[\[CurlyPhi]c,xx[[b]]],{a,Length@xx},{b,Length@xx}];
sumder=Sum[ gup[[a,b]](der[[a,b]]),{a,Length@xx},{b,Length@xx}];
res=1/(8\[Pi]) (2 der-g sumder -1/2 g V\[CurlyPhi]c);,
frame=="Jordan",
dalem=DAlembert[xx,g,\[CurlyPhi]c];
covterm1=Table[CovDer[xx,g,\[CurlyPhi]c,{i},"Valence"->"Covariant","Verbose"->False],{i,Length@xx}];
covterm2=Table[CovDer[xx,g,covterm1,{i,j},"Valence"->"Covariant","Verbose"->False],{i,Length@xx},{j,Length@xx}];
\[Phi]term2=covterm2-g*dalem;
res=If[Not@NumericQ@pert,Normal@Series[1/(8\[Pi])(\[Phi]term2-V\[CurlyPhi]c/2 g),{pert,0,perti}],1/(8\[Pi])(\[Phi]term2-V\[CurlyPhi]c/2 g)];,
True,
Print["Wrong option for Frame"];Return[];
];
If[NumberQ@pert,Simplify@res,Simplify@Normal@Series[res,{pert,0,perti}]]
]
(* Are you sure this is right??? *)
EinsteinST[xx_,g_,V\[Phi]_]:=Block[{res,fRterm1,fRterm2,dalem,R,riscal,riscalvars,\[Phi]c,fRc,dfRc,covterm1,covterm2},
riscal=Simplify@RicciScalar[xx,g];
riscalvars=Complement[xx,Complement[xx,AtomsList[V\[Phi]]]];
\[Phi]c=\[Phi]@@riscalvars;
dalem=DAlembert[xx,g,\[Phi]c];
fRterm1=Simplify[RicciTensor[xx,g]-1/2 g D[V\[Phi],\[Phi]c]+ 1/2 g V\[Phi]/(2 \[Phi]c )];
covterm1=Table[CovDer[xx,g,dfRc,{i},"Valence"->"Covariant","Verbose"->False],{i,Length@xx}];
covterm2=Table[CovDer[xx,g,covterm1,{i,j},"Valence"->"Covariant","Verbose"->False],{i,Length@xx},{j,Length@xx}];
fRterm2=-(1/\[Phi]c) (covterm2-g*dalem);
res=(fRterm1-fRterm2);
Simplify@res
]
Options[STTOV]=Join[{"Signature"->1},Options[TeffST]];
STTOV[coors_,met_,{V\[CurlyPhi]_,var\[CurlyPhi]_},metvars_,pert_:0,OptionsPattern[]]:=Block[{Global`x,Aa,\[Alpha]a,dV\[CurlyPhi]c,dpc,d\[Phi]c,d\[Phi]\[CurlyPhi]c,eq,eqaux,eqKG,eqp,eval,ET,eq2,frame,g,gup,perti,riscal,riscalvars,sign,simpl,T,TEF,Teff,Ttot,T\[Mu]\[Nu],T\[Mu]\[Nu]up,V\[CurlyPhi]c,\[CurlyPhi]c,\[CurlyPhi]2c,Global`r,Global`\[Rho],\[Rho]c,Global`p,pc,Global`\[Phi],\[Phi]c,\[Phi]\[CurlyPhi]c,Global`\[CurlyPhi]},
sign=OptionValue["Signature"];
frame=OptionValue["Frame"];
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
Which[frame!= "Einstein" && frame!="Jordan",Return[];];
Print["Variables must be given as: {p,Var_gtt,Var_grr,\[CurlyPhi]}"];
g=met;
gup=Inverse[g];
riscal=simpl@RicciScalar[coors,g,pert,"PerturbationIndex"->perti];
riscalvars=Complement[coors,Complement[coors,AtomsList[riscal]]];
\[CurlyPhi]c=var\[CurlyPhi]@@riscalvars;
V\[CurlyPhi]c=V\[CurlyPhi]/.var\[CurlyPhi]->\[CurlyPhi]c;
dV\[CurlyPhi]c=D[V\[CurlyPhi]c,\[CurlyPhi]c];
\[CurlyPhi]2c=(var\[CurlyPhi]'')@@riscalvars;
pc=Global`p@@riscalvars;
\[Rho]c=Global`\[Rho]@@riscalvars;\[Phi]c=Global`\[Phi]@@riscalvars;
dpc=D[pc,riscalvars];
d\[Phi]c=D[\[Phi]c,riscalvars];
(* Define the coupling of the scalar field with matter *)
\[Phi]\[CurlyPhi]c[\[CurlyPhi]_]:=Exp[2/Sqrt[3]\[CurlyPhi]];
d\[Phi]\[CurlyPhi]c=D[\[Phi]\[CurlyPhi]c,Global`r];
eval[x_]:=Evaluate[x];
(* Computation of the matter energy-momentum tensor *)
T\[Mu]\[Nu]=ETensor[coors,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
T\[Mu]\[Nu]up=Normal@Series[(gup.T\[Mu]\[Nu])gup,{pert,0,perti}];
T=Normal@Series[Sum[(g.T\[Mu]\[Nu]up)[[i,i]],{i,Length@coors}],{pert,0,perti}];
If[frame=="Einstein",
TEF=(T* Exp[-4/Sqrt[3]\[CurlyPhi]c]);,
TEF=(T)];
(* Computation of the effective energy-momentum tensor. Aa converts matter quantities to the Jordan-Frame. *)
Teff=TeffST[coors,g,{V\[CurlyPhi],var\[CurlyPhi]},pert,"Frame"->frame,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
If[frame=="Einstein",
Aa[\[CurlyPhi]_]:=Exp[-\[CurlyPhi]/Sqrt[3]];
\[Alpha]a[\[CurlyPhi]_]:=D[Log[Aa[\[CurlyPhi]]],\[CurlyPhi]];
Ttot=FullSimplify[8\[Pi](T\[Mu]\[Nu] Aa[\[CurlyPhi]c]^4 +Teff)];
;,
Aa[\[CurlyPhi]_]:=1;
\[Alpha]a[\[CurlyPhi]_]:=1;
Ttot=FullSimplify[8\[Pi]/\[CurlyPhi]c(T\[Mu]\[Nu] +Teff)];];
(* Einstein tensor *)
ET=Einstein[coors,g,pert,"PerturbationIndex"->perti];
(* Solving the equations *)
eq=ET-Ttot;
eq2=Simplify[{Solve[eq[[1]]==0,metvars[[3]]'@@riscalvars][[1,1]],Solve[eq[[2]]==0,metvars[[2]]'@@riscalvars][[1,1]]}];
(* Continuity equation and conversion to EF if frame\[Rule]Einstein *)
eqp=Flatten[Solve[Table[Sum[CovDer[coors,g,T\[Mu]\[Nu]up,{i,i,j},"Valence"->"Contravariant"],{i,Length@coors}],{j,Length@coors}]=={0,0,0,0},metvars[[1]]'@@riscalvars]];
If[frame=="Einstein",
eqaux=Table[T*\[Alpha]a[\[CurlyPhi]c] Sum[gup[[a,b]]D[\[CurlyPhi]c,coors[[a]]],{a,4}],{b,4}];
eqp=Equal@@@Flatten[Solve[Table[Sum[CovDer[coors,g,T\[Mu]\[Nu]up,{i,i,j},"Valence"->"Contravariant"],{i,Length@coors}],{j,Length@coors}]==eqaux,dpc]];
eqp=Simplify@Solve[eqp/.metvars[[1]]->Function[r,Global`\[Phi][r]^(-2)*metvars[[1]][r]]/.\[Rho]c->(\[Rho]c \[Phi]c^(-2))/.Global`\[Phi]->Function[r,Evaluate[\[Phi]\[CurlyPhi]c[Global`\[CurlyPhi][r]]]],dpc];
eqp=Simplify@Normal@Series[Flatten[eqp],{pert,0,perti}];,
eqp=Simplify[eqp];];
(* KG equation *)
If[frame=="Einstein",
eqKG=Solve[DAlembert[coors,g,\[CurlyPhi]c]==-4\[Pi] \[Alpha]a[\[CurlyPhi]c]TEF+1/4 dV\[CurlyPhi]c, \[CurlyPhi]2c];,
eqKG=Solve[3DAlembert[coors,g,\[CurlyPhi]c]==8\[Pi] \[Alpha]a[\[CurlyPhi]c]TEF+ \[CurlyPhi]c dV\[CurlyPhi]c -2 V\[CurlyPhi]c, \[CurlyPhi]2c];];
Normal@Series[Equal@@@Flatten[Join[eqp,eq2,eqKG]],{pert,0,perti}]
]
Options[FRTOV]=Join[{"Signature"->1},Options[TeffST]];
FRTOV[coors_,met_,fR_,metvars_,pert_,OptionsPattern[]]:=Block[{Global`R,teff,perti,sign,u,udown,g,gup,simpl,T\[Mu]\[Nu]up,T\[Mu]\[Nu],T,Teff,Ttot,ET,riscal,riscalvars,Rc,fRc,dfRc,eq,eq2,eqp,eqKG,R2c},
sign=OptionValue["Signature"];
perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
Print["Variables must be given as: {p,Var_gtt,Var_grr,R}"];
g=met;
gup=Inverse[g];
riscal=RicciScalar[coors,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
riscalvars=Complement[coors,Complement[coors,AtomsList[riscal]]];
Rc=Global`R@@riscalvars;
R2c=(Global`R'')@@riscalvars;
fRc=fR/.Global`R->Rc;
dfRc=D[fRc,Rc];
(* Computation of the matter energy-momentum tensor *)T\[Mu]\[Nu]=ETensor[coors,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
T\[Mu]\[Nu]up=(gup.T\[Mu]\[Nu])gup;
T=Sum[(g.T\[Mu]\[Nu]up)[[i,i]],{i,Length@coors}];
(* Computation of the effective energy-momentum tensor *)
Teff=TeffFR[coors,g,fR,pert,"PerturbationIndex"->perti];
Ttot=FullSimplify[8\[Pi]/dfRc(T\[Mu]\[Nu]+Teff)];
(* Einstein tensor *)
ET=Einstein[coors,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
(* Solving the equations *)
eq=ET-Ttot;
eq2=Simplify[{Solve[eq[[1]]==0,metvars[[3]]'@@riscalvars][[1,1]],Solve[eq[[2]]==0,metvars[[2]]'@@riscalvars][[1,1]]}];
(* Continuity equation *)
eqp=Flatten[Solve[Table[Sum[CovDer[coors,g,T\[Mu]\[Nu]up,{i,i,j},"Valence"->"Contravariant"],{i,Length@coors}],{j,Length@coors}]=={0,0,0,0},metvars[[1]]'@@riscalvars]];
(* KG equation *)
eqKG=Solve[(3DAlembert[coors,g,dfRc]+dfRc Rc-2fRc)==8\[Pi] T,R2c];
Equal@@@Flatten[Join[eqp,eq2,eqKG]]
]
fR2Pot[fR_,\[Phi]\[CurlyPhi]_:1]:=Block[{dfRc,fRc,fun,Rc,uc,Global`R,Global`\[Phi],\[Phi]c,Rc\[Phi]},
Rc=Global`R;
\[Phi]c=Global`\[Phi];
fRc=fR/.Global`R->Rc;
dfRc=D[fRc,Rc];
Rc\[Phi]=(Rc/.Solve[\[Phi]c==dfRc,Rc]);
uc=Simplify[Rc\[Phi] \[Phi]c - (fRc/.Rc->Rc\[Phi])];
(* pot=1 gives the JF potential else gives the SF *)
$Assumptions = _ \[Element] Reals;
If[NumericQ@\[Phi]\[CurlyPhi],fun=uc,fun=FullSimplify[(uc/\[Phi]c^2)/.\[Phi]c->\[Phi]\[CurlyPhi]];];
fun
];
fR2\[Phi][fR_,\[Phi]\[CurlyPhi]_:1]:=Block[{dfRc,fRc,fun,Rc,uc,Global`R,Global`\[Phi],\[Phi]c,Rc\[Phi]},
Rc=Global`R;
\[Phi]c=Global`\[Phi];
fRc=fR/.Global`R->Rc;
dfRc=D[fRc,Rc];
Rc\[Phi]=(Rc/.Solve[\[Phi]c==dfRc,Rc]);
$Assumptions = _ \[Element] Reals;
Rc\[Phi][[1]]
];
Pot2fR[V\[CurlyPhi]_,\[Phi]\[CurlyPhi]_]:=Block[{dfRc,fRc,fun,Rc,uc,Global`R,Global`\[Phi],Global`\[CurlyPhi],\[Phi]c,Rc\[Phi],\[CurlyPhi]c},
Rc=Global`R;
\[Phi]c=Global`\[Phi];
\[CurlyPhi]c=Global`\[CurlyPhi];
fRc=fR/.Global`R->Rc;
dfRc=D[fRc,Rc];
Return[Solve[\[Phi]\[CurlyPhi]==\[Phi]c,\[CurlyPhi]c]];
Rc\[Phi]=(Rc/.Solve[\[Phi]c==dfRc,Rc])[[1]];
uc=Simplify[Rc\[Phi] \[Phi]c - (fRc/.Rc->Rc\[Phi])];
(* pot=1 gives the JF potential else gives the SF *)
$Assumptions = _ \[Element] Reals;If[NumericQ@\[Phi]\[CurlyPhi],fun=uc,fun=FullSimplify[(uc/\[Phi]c^2)/.\[Phi]c->\[Phi]\[CurlyPhi]];];
fun
];
EoSCallParsBSks[model_]:=
Block[{a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,table,pars},
table={{a1,3.916`,4.078`,4.857`},{a2,7.701`,7.587`,6.981`},{a3,0.00858`,0.00839`,0.00706`},{a4,0.22114`,0.21695`,0.19351`},{a5,3.269`,3.614`,4.085`},{a6,11.964`,11.942`,12.065`},{a7,13.349`,13.751`,10.521`},{a8,1.3683`,1.3373`,1.5905`},{a9,3.254`,3.606`,4.104`},{a10,-12.953`,-22.996`,-28.726`},{a11,0.9237`,1.6229`,2.0845`},{a12,6.2`,4.88`,4.89`},{a13,14.383`,14.274`,14.302`},{a14,16.693`,23.56`,22.881`},{a15,-1.0514`,-1.5564`,-1.769`},{a16,2.486`,2.095`,0.989`},{a17,15.362`,15.294`,15.313`},{a18,0.085`,0.084`,0.091`},{a19,6.23`,6.36`,4.68`},{a20,11.68`,11.67`,11.65`},{a21,-0.029`,-0.042`,-0.086`},{a22,20.1`,14.8`,10.`},{a23,14.19`,14.18`,14.15`}};
pars=Table[table[[i,1]]->table[[i,model+1]],{i,Length@table}]
]
EoSCallParsSly[model_]:=
Block[{a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,table,pars},
table={{a1, 6.22},{a10,11.4950},
{a2,6.121},{a11,-22.775},
{a3,0.005925},{a12,1.5707},
{a4,0.16326},{a13,4.3 },
{a5,6.48},{a14,14.08 },
{a6,11.4971},{a15,27.80 },
{a7,19.105},{a16,-1.653},
{a8,0.8938},{a17,1.50},
{a9,6.54},{a18,14.67}};
pars=Table[table[[i,1]]->table[[i,model+1]],{i,Length@table}]
]
EoSBSks[model_]:=Block[{\[Zeta],\[Rho],pars},
pars=EoSCallParsBSks[model];
\[Zeta]=(a1+a2 \[Rho] +a3 \[Rho]^3)/(1+a4 \[Rho]) (Exp[a5(\[Rho]-a6)]+1)^(-1)+(a7+a8 \[Rho])(Exp[a9(a6-\[Rho])]+1)^(-1)
+(a10+a11 \[Rho])(Exp[a12(a13-\[Rho])]+1)^(-1)+(a14+a15 \[Rho])(Exp[a16(a17-\[Rho])]+1)^(-1)
+a18/(1+(a19(\[Rho]-a20))^2)+a21/(1+(a22(\[Rho]-a20))^2);
\[Zeta]/.pars
]
EoSSly[model_]:=Block[{\[Zeta],\[Rho],pars,x,f0},
pars=EoSCallParsSly[model];
f0[x_]:=1/(Exp[x]+1);
\[Zeta]=(a1+a2 \[Rho] +a3 \[Rho]^3)/(1+a4 \[Rho]) f0[a5(\[Rho]-a6)]+
(a7+a8 \[Rho])f0[(a9(a10-\[Rho]))]+
(a11+a12 \[Rho])f0[(a13(a14-\[Rho]))]+
(a15+a16 \[Rho])f0[(a17(a18-\[Rho]))];
\[Zeta]/.pars
]
SlyInner={{0.00020905`,3.4951`*^11,6.214999999999999`*^29,1.177`,0.0099795`,1.6774`*^13,3.072`*^31,1.342`},{0.00022059000000000003`,3.6883`*^11,6.4304`*^29,0.527`,0.012513000000000002`,2.1042`*^13,4.1574`*^31,1.332`},{0.00023114`,3.865`*^11,6.5813`*^29,0.476`,0.016547`,2.7844000000000004`*^13,6.023399999999999`*^31,1.322`},{0.00026426`,4.4199`*^11,6.9945`*^29,0.447`,0.021405`,3.6043`*^13,8.461299999999999`*^31,1.32`},{0.00030533000000000003`,5.1079999999999994`*^11,7.4685`*^29,0.466`,0.024157`,4.068800000000001`*^13,9.928599999999999`*^31,1.325`},{0.00035331`,5.9119`*^11,8.0149`*^29,0.504`,0.027894000000000002`,4.7001`*^13,1.2022999999999999`*^32,1.338`},{0.00040763999999999997`,6.8224`*^11,8.644299999999999`*^29,0.554`,0.031941000000000004`,5.3843`*^13,1.4430000000000001`*^32,1.358`},{0.000468`,7.8339`*^11,9.3667`*^29,0.61`,0.036264`,6.115300000000001`*^13,1.7175`*^32,1.387`},{0.00053414`,8.9426`*^11,1.0190999999999999`*^30,0.668`,0.039888`,6.7284`*^13,1.9626`*^32,1.416`},{0.00060594`,1.0146`*^12,1.1128`*^30,0.726`,0.044578`,7.5224`*^13,2.3024`*^32,1.458`},{0.00076608`,1.2831`*^12,1.337`*^30,0.84`,0.048425`,8.1738`*^13,2.6018`*^32,1.496`},{0.0010471`,1.7543`*^12,1.7792`*^30,0.987`,0.052327000000000005`,8.835000000000002`*^13,2.9261000000000002`*^32,1.536`},{0.0012616`,2.1141`*^12,2.1547000000000002`*^30,1.067`,0.056264`,9.5022`*^13,3.2756`*^32,1.576`},{0.0016246000000000001`,2.7232`*^12,2.8565000000000003`*^30,1.16`,0.060218999999999995`,1.0173000000000002`*^14,3.6505000000000004`*^32,1.615`},{0.0020384`,3.4178`*^12,3.7461000000000004`*^30,1.227`,0.064183`,1.0845`*^14,4.0509000000000005`*^32,1.65`},{0.0026726000000000002`,4.4827`*^12,5.2679`*^30,1.286`,0.067163`,1.1351`*^14,4.3681`*^32,1.672`},{0.0034064`,5.7153`*^12,7.230400000000001`*^30,1.322`,0.070154`,1.1859`*^14,4.6998`*^32,1.686`},{0.0044746`,7.5106`*^12,1.0404999999999999`*^31,1.344`,0.073174`,1.2372`*^14,5.0462`*^32,1.685`},{0.005726`,9.6148`*^12,1.4513`*^31,1.353`,0.075226`,1.272`*^14,5.2856`*^32,1.662`},{0.0074963`,1.2593`*^13,2.0894`*^31,1.351`,0.075959`,1.2845`*^14,5.3739`*^32,1.644`}};
SlyLCore={{0.0771`,1.3038`*^14,5.3739`*^32,2.159`,0.49`,8.850899999999999`*^14,1.0315`*^35,2.953`},{0.08`,1.3531`*^14,5.8259999999999996`*^32,2.217`,0.52`,9.4695`*^14,1.2289`*^35,2.943`},{0.085`,1.4381`*^14,6.6828`*^32,2.309`,0.55`,1.0102`*^15,1.4491`*^35,2.933`},{0.09`,1.5232`*^14,7.6443`*^32,2.394`,0.58`,1.0748`*^15,1.693`*^35,2.924`},{0.1`,1.6935`*^14,9.9146`*^32,2.539`,0.61`,1.1408`*^15,1.9616`*^35,2.916`},{0.11`,1.8641`*^14,1.2700999999999999`*^33,2.655`,0.64`,1.2085`*^15,2.2558999999999998`*^35,2.908`},{0.12`,2.035`*^14,1.6063`*^33,2.708`,0.67`,1.2777`*^15,2.5769`*^35,2.9`},{0.13`,2.2063`*^14,1.9971`*^33,2.746`,0.7`,1.3486`*^15,2.9255`*^35,2.893`},{0.16`,2.7223000000000003`*^14,3.5926999999999995`*^33,2.905`,0.75`,1.4706`*^15,3.5702`*^35,2.881`},{0.19`,3.2424`*^14,5.9667`*^33,2.99`,0.8`,1.5977`*^15,4.2981`*^35,2.869`},{0.22`,3.7675`*^14,9.2766`*^33,3.025`,0.85`,1.7302`*^15,5.1128999999999996`*^35,2.858`},{0.25`,4.2983`*^14,1.3668`*^34,3.035`,0.9`,1.8683`*^15,6.0183`*^35,2.847`},{0.28`,4.8358`*^14,1.9276999999999999`*^34,3.032`,0.95`,2.0123000000000002`*^15,7.0176`*^35,2.836`},{0.31`,5.3808`*^14,2.6234999999999996`*^34,3.023`,1.`,2.1624`*^15,8.113899999999998`*^35,2.824`},{0.34`,5.934`*^14,3.467`*^34,3.012`,1.1`,2.482`*^15,1.0609`*^36,2.801`},{0.37`,6.4963`*^14,4.4702`*^34,2.999`,1.2`,2.8289`*^15,1.3524`*^36,2.778`},{0.4`,7.0684`*^14,5.6451`*^34,2.987`,1.3`,3.2048`*^15,1.6876`*^36,2.754`},{0.43`,7.651`*^14,7.0033`*^34,2.975`,1.4`,3.6113`*^15,2.0678999999999998`*^36,2.731`},{0.46`,8.244999999999999`*^14,8.5561`*^34,2.964`,1.5`,4.0498000000000005`*^15,2.4947`*^36,2.708`}};
SlyLCoreAll=Join[TakeColumn[SlyLCore,{2,3}],TakeColumn[SlyLCore,{6,7}]]/.{zz_,yy_}->{Log[10,zz],Log[10,yy]};
SlyInnerAll=Join[TakeColumn[SlyInner,{2,3}],TakeColumn[SlyInner,{6,7}]]/.{zz_,yy_}->{Log[10,zz],Log[10,yy]};
EoSFitsPars[model_,verbose_:False]:=Module[{eostable},
eostable={{"PAL6",34.38`,2.227`,2.189`,2.159`,0.0011`,0.693`,1.37`,1.477`,-0.47`,0.374`,-0.51`,1660,-0.97`,1.051`,-2.03`,10.547`,-0.54`},{"SLy",34.384`,3.005`,2.988`,2.851`,0.002`,0.989`,1.41`,2.049`,0.02`,0.592`,0.81`,1810,0.1`,1.288`,-0.08`,11.736`,-0.21`},{"APR1",33.943`,2.442`,3.256`,2.908`,0.019`,0.924`,9.94`,1.683`,-1.6`,0.581`,2.79`,2240,1.05`,0.908`,-2.57`,9.361`,-1.85`},{"APR2",34.126`,2.643`,3.014`,2.945`,0.0089`,1.032`,0.42`,1.808`,-1.5`,0.605`,0.33`,2110,-0.02`,1.024`,-2.34`,10.179`,-1.57`},{"APR3",34.392`,3.166`,3.573`,3.281`,0.0091`,1.134`,2.72`,2.39`,-1.`,0.704`,0.57`,1810,-0.14`,1.375`,-1.59`,12.094`,-0.96`},{"APR4",34.269`,2.83`,3.445`,3.348`,0.0068`,1.16`,1.45`,2.213`,-1.08`,0.696`,0.22`,1940,0.05`,1.243`,-1.36`,11.428`,-0.9`},{"FPS",34.283`,2.985`,2.863`,2.6`,0.005`,0.883`,2.29`,1.799`,-0.03`,0.53`,0.67`,1880,0.11`,1.137`,0.03`,10.85`,0.12`},{"WFF1",34.031`,2.519`,3.791`,3.66`,0.018`,1.185`,7.86`,2.133`,-0.29`,0.739`,2.21`,2040,0.3`,1.085`,0.1`,10.414`,0.02`},{"WFF2",34.233`,2.888`,3.475`,3.517`,0.017`,1.139`,7.93`,2.198`,-0.14`,0.717`,0.71`,1990,0.03`,1.204`,-0.59`,11.159`,-0.28`},{"WFF3",34.283`,3.329`,2.952`,2.589`,0.017`,0.835`,8.11`,1.844`,-0.48`,0.53`,2.26`,1860,0.59`,1.16`,-0.25`,10.926`,-0.12`},{"BBB2",34.331`,3.418`,2.835`,2.832`,0.0055`,0.914`,7.75`,1.918`,0.1`,0.574`,0.97`,1900,0.47`,1.188`,0.17`,11.139`,-0.29`},{"BPAL12",34.358`,2.209`,2.201`,2.176`,0.001`,0.708`,1.03`,1.452`,-0.18`,0.382`,-0.29`,1700,-1.03`,0.974`,0.2`,10.024`,0.67`},{"ENG",34.437`,3.514`,3.13`,3.168`,0.015`,1.`,10.71`,2.24`,-0.05`,0.654`,0.39`,1820,-0.44`,1.372`,-0.97`,12.059`,-0.69`},{"MPA1",34.495`,3.446`,3.572`,2.887`,0.0081`,0.994`,4.91`,2.461`,-0.16`,0.67`,-0.05`,1700,-0.18`,1.455`,-0.41`,12.473`,-0.26`},{"MS1",34.858`,3.224`,3.033`,1.325`,0.019`,0.888`,12.44`,2.767`,-0.54`,0.606`,-0.52`,1400,1.67`,1.944`,-0.09`,14.918`,0.06`},{"MS2",34.605`,2.447`,2.184`,1.855`,0.003`,0.582`,3.96`,1.806`,-0.42`,0.343`,2.57`,1250,2.25`,1.658`,0.46`,14.464`,-2.69`},{"MS1b",34.855`,3.456`,3.011`,1.425`,0.015`,0.889`,11.38`,2.776`,-1.03`,0.614`,-0.56`,1420,1.38`,1.888`,-0.64`,14.583`,-0.32`},{"PS",34.671`,2.216`,1.64`,2.365`,0.028`,0.691`,7.36`,1.755`,-1.53`,0.355`,-1.45`,1300,-2.39`,2.067`,-3.06`,15.472`,3.72`},{"GS1",34.504`,2.35`,1.267`,2.421`,0.018`,0.695`,0.49`,1.382`,-1.`,0.395`,-0.64`,1660,9.05`,0.766`,-3.13`,Null},{"GS2",34.642`,2.519`,1.571`,2.314`,0.026`,0.592`,16.1`,1.653`,-0.3`,0.339`,7.71`,1340,3.77`,1.795`,-3.33`,14.299`,0.07`},{"BGN1H1",34.623`,3.258`,1.472`,2.464`,0.029`,0.878`,-7.42`,1.628`,0.39`,0.437`,-3.55`,1670,-2.08`,1.504`,0.56`,12.901`,-1.96`},{"GNH3",34.648`,2.664`,2.194`,2.304`,0.0045`,0.75`,2.04`,1.962`,0.13`,0.427`,0.37`,1410,-0.04`,1.713`,0.38`,14.203`,-0.28`},{"H1",34.564`,2.595`,1.845`,1.897`,0.0019`,0.561`,2.81`,1.555`,-0.92`,0.311`,-1.47`,1320,-1.46`,1.488`,-1.45`,12.861`,-0.03`},{"H2",34.617`,2.775`,1.855`,1.858`,0.0028`,0.565`,1.38`,1.666`,-0.77`,0.322`,-0.55`,1280,-1.29`,1.623`,-0.82`,13.479`,0.29`},{"H3",34.646`,2.787`,1.951`,1.901`,0.007`,0.564`,7.05`,1.788`,-0.79`,0.343`,1.07`,1290,-0.88`,1.702`,-1.18`,13.84`,0.31`},{"H4",34.669`,2.909`,2.246`,2.144`,0.0028`,0.685`,4.52`,2.032`,-0.85`,0.428`,-1.01`,1400,-1.28`,1.729`,-1.18`,13.774`,1.34`},{"H5",34.609`,2.793`,1.974`,1.915`,0.005`,0.596`,1.65`,1.727`,-1.`,0.347`,-0.82`,1340,-1.55`,1.615`,-1.31`,13.348`,0.68`},{"H6",34.593`,2.637`,2.121`,2.064`,0.0087`,0.598`,11.71`,1.778`,0.07`,0.346`,8.65`,1310,5.33`,1.623`,-2.19`,13.463`,0.37`},{"H7",34.559`,2.621`,2.048`,2.006`,0.0046`,0.63`,1.82`,1.683`,-1.12`,0.357`,-0.57`,1410,-1.52`,1.527`,-2.33`,12.992`,0.23`},{"PCL2",34.507`,2.554`,1.88`,1.977`,0.0069`,0.6`,1.74`,1.482`,-0.79`,0.326`,-2.25`,1440,-1.87`,1.291`,-3.27`,11.761`,-1.15`},{"ALF1",34.055`,2.013`,3.389`,2.033`,0.04`,0.565`,18.59`,1.495`,-0.53`,0.386`,3.52`,1730,2.44`,0.987`,-0.4`,9.896`,-0.22`},{"ALF2",34.616`,4.07`,2.411`,1.89`,0.043`,0.642`,1.5`,2.086`,-5.26`,0.436`,-0.62`,1440,1.01`,1.638`,-6.94`,13.188`,-3.66`},{"ALF3",34.283`,2.883`,2.653`,1.952`,0.017`,0.565`,11.29`,1.473`,-0.06`,0.358`,2.46`,1620,1.79`,1.041`,0.76`,10.314`,-0.25`},{"ALF4",34.314`,3.009`,3.438`,1.803`,0.023`,0.685`,14.78`,1.943`,-0.93`,0.454`,0.59`,1590,0.52`,1.297`,-2.38`,11.667`,-1.2`}};If[verbose,eostable[[All,1]],Flatten@Select[eostable,#[[1]]==model&]]
]
EoSSlyCrust[\[Rho]_]:=Module[{pol1,pol2,pol3,pol4,\[Rho]l1,\[Rho]l2,\[Rho]l3,k1,k2,k3,k4,\[CapitalGamma]1,\[CapitalGamma]2,\[CapitalGamma]3,\[CapitalGamma]4,c},
c=2.99792458 10^10;
\[Rho]l1=Log[10,2.44034 10^(07)];
\[Rho]l2=Log[10,3.78358 10^(11)];
\[Rho]l3=Log[10,2.62780 10^(12)];
k1=Log[10,(6.80110 10^(-9))*c^2];
k2=Log[10,(1.06186 10^(-6))*c^2];
k3=Log[10,(5.32697 10)*c^2];
k4=Log[10,(3.99874 10^(-8))*c^2];
\[CapitalGamma]1=1.58425;
\[CapitalGamma]2=1.28733;
\[CapitalGamma]3=0.62223;
\[CapitalGamma]4=1.35692;
pol1=k1+ \[CapitalGamma]1 \[Rho];
pol2=k2+ \[CapitalGamma]2 \[Rho];
pol3=k3+ \[CapitalGamma]3 \[Rho];
pol4=k4+ \[CapitalGamma]4 \[Rho];
Piecewise[{{pol1,\[Rho]<=\[Rho]l1},{pol2,\[Rho]l1<\[Rho]<=\[Rho]l2},{pol3,\[Rho]l2<\[Rho]<=\[Rho]l3},{pol4,\[Rho]>\[Rho]l3}}]
];
EoSPol[model_:"PolR"]:=Block[{k,\[CapitalGamma],pol,Global`\[Rho],c,pol1,pol2,k2},
If[ListQ[model],k=Log10[model[[1]]];\[CapitalGamma]=model[[2]];pol=k+ \[CapitalGamma] Global`\[Rho],
Which[model=="PolNR",k=Log[10,(3.3)];\[CapitalGamma]=2; pol=k+ \[CapitalGamma] Global`\[Rho];,
model=="PolR",k=Log[10,(1.98183*10^-6)];\[CapitalGamma]=2.75; pol=k+ \[CapitalGamma] Global`\[Rho];,
model=="PolMS",k=Log[10,3.849119840037`*^14];\[CapitalGamma]=4/3; pol= k +\[CapitalGamma] Global`\[Rho];,
model=="PolMSMix",k=Log[10,3.849119840037`*^14]; pol1=k+ 4/3 Global`\[Rho];k2=Log[10,1.2392481667219052`*^15];pol2= k2+ 5/3 Global`\[Rho]; pol=Piecewise[{{pol1,Global`\[Rho]>Log10[0.029964]},{pol2,Global`\[Rho]<=Log10[0.029964]}}];,
True,Return[]];
];
pol
]
EoSPol\[Epsilon][model_:"PolR"]:=Block[{k,\[CapitalGamma],pol,Global`\[Rho],c,pol1,pol2,k2},
c=2.99792458 10^10;
If[ListQ[model],k=Log10[model[[1]]];\[CapitalGamma]=model[[2]];pol=k+ \[CapitalGamma] Global`\[Rho],
Which[model=="PolNR",k=Log[10,(3.3)];\[CapitalGamma]=2; pol=k+ \[CapitalGamma] Global`\[Rho],
model=="PolR",k=Log[10,(1.98183*10^-6)];\[CapitalGamma]=2.75; pol=k+ \[CapitalGamma] Global`\[Rho],
model=="PolMS",k=Log[10,3.849119840037`*^14];\[CapitalGamma]=4/3; pol=k+ \[CapitalGamma] Global`\[Rho];];
];
Global`\[Rho] +10^k/(c^2(\[CapitalGamma]-1))Global`\[Rho]^(\[CapitalGamma])
]
Options[EoSFits]:={"PhysUnits"->False,"Verbose"->False}
EoSFits[model_,OptionsPattern[]]:=Block[{P0Sly=Log[10,5.37*10^(32)],\[Rho]0Sly=Log[10,1.2845 10^14],\[CapitalGamma]0=4/3,\[Rho]1=Log[10,10^(14.7)],\[Rho]2=Log[10,10^(15.)],fit0,fit1,fit2,
fit3,Global`\[Rho],p,pars,p1, \[CapitalGamma]1,\[CapitalGamma]2,\[CapitalGamma]3,K0,K1,K2,K3,\[Rho]0c,p2,verbose,picore,res,\[Rho]1c,\[Rho]2c,\[Rho]3c,k1c,k2c,k3c,k4c,\[CapitalGamma]1c,\[CapitalGamma]2c,\[CapitalGamma]3c,\[CapitalGamma]4c,fit1c,fit2c,fit3c,fit4c,physuns,c},
verbose=OptionValue["Verbose"];
physuns=OptionValue["PhysUnits"];
If[verbose,Print[" EoS available "];Return[EoSFitsPars[model,True]]];
c=2.99792458 10^10;
(* values for the crust taken from SLy crust model *)
\[Rho]1c=Log[10,2.44034 10^(07)];
\[Rho]2c=Log[10,3.78358 10^(11)];
\[Rho]3c=Log[10,2.62780 10^(12)];
k1c=Log[10,(6.80110 10^(-9))*c^2];
k2c=Log[10,(1.06186 10^(-6))*c^2];
k3c=Log[10,(5.32697 10)*c^2];
k4c=Log[10,(3.99874 10^(-8))*c^2];
\[CapitalGamma]1c=1.58425;
\[CapitalGamma]2c=1.28733;
\[CapitalGamma]3c=0.62223;
\[CapitalGamma]4c=1.35692;
fit1c=k1c+ \[CapitalGamma]1c Global`\[Rho];
fit2c=k2c+ \[CapitalGamma]2c Global`\[Rho];
fit3c=k3c+ \[CapitalGamma]3c Global`\[Rho];
fit4c=k4c+ \[CapitalGamma]4c Global`\[Rho];
(* p1, \[CapitalGamma]1,\[CapitalGamma]2,\[CapitalGamma]3*)
{p1, \[CapitalGamma]1,\[CapitalGamma]2,\[CapitalGamma]3}=EoSFitsPars[model][[2;;5]];
(* Using a polytropic SLy for the crust *)
K1=(p1-(\[CapitalGamma]1 \[Rho]1));
fit0=Piecewise[{{fit1c,Global`\[Rho]<=\[Rho]1c},{fit2c,\[Rho]1c<Global`\[Rho]<=\[Rho]2c},{fit3c,\[Rho]2c<Global`\[Rho]<=\[Rho]3c},{fit4c,Global`\[Rho]>\[Rho]3c}}];
fit1=K1+(\[CapitalGamma]1 Global`\[Rho]);
\[Rho]0c=(Global`\[Rho]/.Solve[fit0==K1+\[CapitalGamma]1 Global`\[Rho],Global`\[Rho]])[[1]];
K2=p1-(\[CapitalGamma]2 \[Rho]1);
fit2=K2+(\[CapitalGamma]2 Global`\[Rho]);
p2=K2+(\[CapitalGamma]2 \[Rho]2);
K3=p2-(\[CapitalGamma]3 \[Rho]2);
fit3=K3+(\[CapitalGamma]3 Global`\[Rho]);
If[Not@physuns,
res=Piecewise[{{fit1c,Global`\[Rho]<=\[Rho]1c},{fit2c,\[Rho]1c<Global`\[Rho]<=\[Rho]2c},{fit3c,\[Rho]2c<Global`\[Rho]<=\[Rho]3c},{fit4c,\[Rho]3c<Global`\[Rho]<=\[Rho]0c},{fit1,\[Rho]0c<Global`\[Rho]<=\[Rho]1},{fit2,\[Rho]1<Global`\[Rho]<=\[Rho]2},{fit3,Global`\[Rho]>\[Rho]2}}];,
res=Piecewise[{{10^(fit1c/.Global`\[Rho]->Log[10,Global`\[Rho]]),Global`\[Rho]<=10^\[Rho]1c},{10^(fit2c/.Global`\[Rho]->Log[10,Global`\[Rho]]),10^\[Rho]1c<Global`\[Rho]<=10^\[Rho]2c},{10^(fit3c/.Global`\[Rho]->Log[10,Global`\[Rho]]),10^\[Rho]2c<Global`\[Rho]<=10^\[Rho]3c},{10^(fit4c/.Global`\[Rho]->Log[10,Global`\[Rho]]),10^\[Rho]3c<Global`\[Rho]<=10^\[Rho]0c},{10^(fit1/.Global`\[Rho]->Log[10,Global`\[Rho]]),10^\[Rho]0c<Global`\[Rho]<=10^\[Rho]1},{10^(fit2/.Global`\[Rho]->Log[10,Global`\[Rho]]),10^\[Rho]1<Global`\[Rho]<=10^\[Rho]2},{10^(fit3/.Global`\[Rho]->Log[10,Global`\[Rho]]),Global`\[Rho]>10^\[Rho]2}}];
];
res
]
From\[Rho]To\[Epsilon]Fits[eos_]:=Block[{Global`\[Rho],ks,\[Gamma]s,pols,polsd,a,\[Rho]s,\[Epsilon],tab,k,k1,\[Gamma],\[Gamma]1,\[Rho]v,c},
c=2.99792458 10^10;
If[eos=="PolMSMix",pols=EoSPol[eos][[1]],pols=EoSFits[eos][[1]]];
{ks,\[Gamma]s}=Transpose[CoefficientList[#,Global`\[Rho]]&/@pols[[All,1]]];
polsd=pols[[All,2]];
\[Rho]s=DeleteDuplicates@polsd[[All,-1]];
a=0;
\[Epsilon][Global`\[Rho]_,a_,k_,\[Gamma]_]:=(1+a)Global`\[Rho] + k/(c^2(\[Gamma]-1))Global`\[Rho]^(\[Gamma]);
polsd=Table[10^\[Rho]s[[i]]<Global`\[Rho]<=10^\[Rho]s[[i+1]],{i,Length@\[Rho]s-1}];
polsd=Join[{Global`\[Rho]<=10^\[Rho]s[[1]]},polsd,{Global`\[Rho]>10^\[Rho]s[[-1]]}];
tab=Table[k=10^ks[[i]];\[Gamma]=\[Gamma]s[[i]];\[Rho]v=10^\[Rho]s[[i]];k1=10^ks[[i+1]];\[Gamma]1=\[Gamma]s[[i+1]];
a=\[Epsilon][\[Rho]v,a,k,\[Gamma]]/\[Rho]v-1- k1/(c^2(\[Gamma]1-1))\[Rho]v^(\[Gamma]1-1);
{\[Epsilon][Global`\[Rho],a,k1,\[Gamma]1],polsd[[i+1]]}
,{i,1,Length@\[Rho]s}];
Simplify@Piecewise[Join[{{\[Epsilon][Global`\[Rho],0,10^ks[[1]],\[Gamma]s[[1]]],polsd[[1]]}},tab]]
]
RK4[func_?ListQ,vars_?ListQ,ivals_?ListQ,pars_?ListQ,step_]:=Module[{k1,k2,k3,k4,x2,f1,f2,x3,f3,x4,f4,dx,x,x0,sol},
dx=step;
{x,x0}={pars[[1]],pars[[2]]};
k1=dx ((func/.x->(x0))/.MapThread[Rule, {vars,ivals}]);
k2=dx ((func/.x->(x0+dx/2))/.MapThread[Rule, {vars,ivals +1/2 k1}]);
k3=dx ((func/.x->(x0+dx/2))/.MapThread[Rule, {vars,ivals +1/2 k2}]);
k4=dx ((func/.x->(x0+dx))/.MapThread[Rule, {vars,ivals + k3}]);
ivals+1/6(k1+2k2+2k3+k4)
]
TestCode[eqs_,vars_,icond_,rlst_,drlst_]:=Module[{R1,R,v1,v,w1,w,\[Lambda]1,\[Lambda],p1,p,eqsrules,atomlst,atomlstaux,solvevars,varsp1,rvar,rval,drvar,drval},
{rvar,rval}=rlst;
{drvar,drval}=drlst;
eqsrules=(Flatten[eqs/.MapThread[Rule, {vars,icond}]])/.drvar->drval/.rvar->rval;
varsp1=Flatten[eqs][[All,1]];
eqsrules=Flatten@Solve[eqsrules,varsp1];
Return[varsp1/.eqsrules];
]
Options[BracketingSTNStars]={"Tolerance"->10^(-8),"Verbose"->False,"MaxIteraton"->100,"NPoints"->1000,"AssymptoticMatch"->None,"AssymptoticValue"->10^-8};
BracketingSTNStars[eqs_,eqsRg_,Global`r_,vars_,shtdInd_,varshtdRg_,OptionsPattern[]]:=Module[{dom,eqsht,
Sh0,Sh0m,posref,mean,tol,Sh\[Infinity],Sh\[Infinity]m,Sh\[Infinity]m2,Sh0m2,Sh0ref,
a,posreftest,verbose,begin,eqsRga,brack,varshtdRga,varsa,varshta,varshtalw,dvarshta,out,ShtStr,raux,np,i,amax,assymptotic,
Rs,r,datab,datfit,a0,m,y1,y2,A3,assval,threshold,Sh\[Infinity]maux},
(* Loading options *)
tol=OptionValue["Tolerance"];
verbose=OptionValue["Verbose"];
amax=OptionValue["MaxIteraton"];
assymptotic=OptionValue["AssymptoticMatch"];
np=OptionValue["NPoints"];
assval=OptionValue["AssymptoticValue"];
(* Some auxiliary variables *)
eqsRga=eqsRg;
varshtdRga=varshtdRg;
varsa=ToExpression[(ToString[#]<>"a")&/@vars];
varshta=Join[varshtdRga,{Mean[varshtdRga]}];
ShtStr=ToString[varsa[[shtdInd]]];
a=0 (* auxiliary counter *) ;
i=0 (* auxiliary counter *) ;
dvarshta=(varshta[[2]]-varshta[[1]])/np;
(* In case 'bracketing' is activated, we split [varshta[[1]], varshta[[2]]] in np points *)
Sh\[Infinity]m=1.1;
If[assval>=1,threshold=assval,threshold=0];
While[Sh\[Infinity]m>threshold&&i< np,
i=i+1;
Sh\[Infinity]maux=Sh\[Infinity]m;
varshtalw=varshtdRg[[2]] - i*dvarshta;
If[verbose, Print["n: ",i," Redefining upper limit as: ",varshtalw]];
eqsht={vars[[shtdInd]][eqsRga[[1]]]==varshtalw};
eqsht=Join[eqs,eqsht];
varsa=vars/.Flatten[NDSolve[eqsht,vars,{Global`r,eqsRga[[1]],eqsRga[[2]]},AccuracyGoal->15,PrecisionGoal->15,WorkingPrecision->30,MaxSteps->Infinity]];
dom=InterpolationDomain[varsa[[1]]];
Sh\[Infinity]m=varsa[[shtdInd]]@dom[[2]];
];
Return[{varshtalw,varshtalw + dvarshta,i}]
(* The - is to avoid varshta[[1]]=varshta[[2]] *)
(*varshta[[2]]=varshtalw + dvarshta;
If[verbose, Print[" New upper limit: ",varshta[[2]]]];
Sh\[Infinity]m=-0.1;
i=0;
varshtalw=varshta[[1]];
If[assval\[GreaterEqual]1,threshold=assval,threshold=0];
While[Sh\[Infinity]m<threshold && i<np,
i=i+1;
varshtalw=varshtdRg[[1]] +i*dvarshta;
eqsht={vars[[shtdInd]][eqsRga[[1]]]\[Equal]varshtalw};
eqsht=Join[eqs,eqsht];
varsa=vars/.Flatten[NDSolve[eqsht,vars,{Global`r,eqsRga[[1]],eqsRga[[2]]},AccuracyGoal\[Rule]15,PrecisionGoal\[Rule]15,WorkingPrecision\[Rule]30,MaxSteps->Infinity]];
dom=InterpolationDomain[varsa[[1]]];
Sh\[Infinity]m=varsa[[shtdInd]]@dom[[2]]; varshta[[1]]=varshtalw;
If[verbose, Print["n: ",i,". Redefining lower limit as: ",varshta[[1]]]];
];
varshta[[1]]=varshtalw - dvarshta;
Return[Flatten[Join[varshta[[1;;2]],{i}]]];
*)];
Options[ShootingNStars]={"Tolerance"->10^(-8),"Verbose"->True,"Bracketing"->False,"MaxIteraton"->100,"NPoints"->1000,"AssymptoticMatch"->None,"AssymptoticValue"->10^-8};
ShootingNStars[eqs_,eqsRg_,Global`r_,vars_,shtdInd_,varshtdRg_,optNDS__,OptionsPattern[]]:=Module[{dom,eqsht,
Sh0,Sh0m,posref,mean,tol,Sh\[Infinity],Sh\[Infinity]m,Sh\[Infinity]m2,Sh0m2,Sh0ref,
a,posreftest,verbose,begin,eqsRga,brack,varshtdRga,varsa,varshta,varshtalw,dvarshta,out,ShtStr,raux,np,i,amax,assymptotic,
Rs,r,datab,datfit,a0,m,y1,y2,A3,assval,threshold},
(* Loading options *)
tol=OptionValue["Tolerance"];
verbose=OptionValue["Verbose"];
brack=OptionValue["Bracketing"];
amax=OptionValue["MaxIteraton"];
assymptotic=OptionValue["AssymptoticMatch"];
np=OptionValue["NPoints"];
assval=OptionValue["AssymptoticValue"];
(* Some auxiliary variables *)
eqsRga=eqsRg;
varshtdRga=varshtdRg;
varsa=ToExpression[(ToString[#]<>"a")&/@vars];
varshta=Join[varshtdRga,{Mean[varshtdRga]}];
ShtStr=ToString[varsa[[shtdInd]]];
a=0 (* auxiliary counter *) ;
i=0 (* auxiliary counter *) ;
dvarshta=(varshta[[2]]-varshta[[1]])/np;(* In case 'bracketing' is activated, we split [varshta[[1]], varshta[[2]]] in np points *)
(* In case shooting is not required *)
If[Length@varshtdRga==1,
eqsht={vars[[shtdInd]][eqsRga[[1]]]==(varshta[[1]])};
eqsht=Join[eqs,eqsht];
varsa=vars/.Flatten[NDSolve[Chop[eqsht],vars,{Global`r,eqsRga[[1]],eqsRga[[2]]},AccuracyGoal->15,PrecisionGoal->15,WorkingPrecision->30,MaxSteps->Infinity]];
out=Join[varsa,{a,Rationalize[varsa[[shtdInd]]@eqsRga[[1]]],0}];
Return[out];
];
(* We also include provide a 'bracketing' option in case the bracketing on the shooting variable is required *)
If[brack,
varshtalw=Chop[varshta[[2]]];
eqsht={vars[[shtdInd]][eqsRga[[1]]]==varshtalw};
eqsht=Join[eqs,eqsht];
varsa=vars/.Flatten[NDSolve[Chop[eqsht],vars,{Global`r,eqsRga[[1]],eqsRga[[2]]},AccuracyGoal->15,PrecisionGoal->15,WorkingPrecision->30,MaxSteps->Infinity]];
dom=InterpolationDomain[varsa[[1]]][[1]];
(*If[dom[[2]]<eqsRga[[2]],eqsRga[[2]]=Rationalize[0.9*eqsRga[[2]],1];If[verbose,Print[Style[" Breakup found. Changing rfin to ",Red],eqsRga[[2]]] ];Goto[begin];];*)
Which[Sh\[Infinity]m<0&&assval<= 1, Print[Style[" Upper value of the shooted variable is not positive at rmax: Redefine the brackets !",Red]]; Return[];,
Sh\[Infinity]m<1&&assval>= 1, Print[Style[" Upper value of the shooted variable is <1 at rmax: Redefine the brackets !",Red]]; Return[];];
Sh\[Infinity]m=-0.1;
i=0;
If[assval>=1,threshold=assval,threshold=0];
While[Sh\[Infinity]m<threshold && i< np,
i=i+1;
varshtalw=varshtdRg[[1]] +i*dvarshta;
eqsht={vars[[shtdInd]][eqsRga[[1]]]==varshtalw};
eqsht=Join[eqs,eqsht];
varsa=vars/.Flatten[NDSolve[eqsht,vars,{Global`r,eqsRga[[1]],eqsRga[[2]]},AccuracyGoal->15,PrecisionGoal->15,WorkingPrecision->30,MaxSteps->Infinity]];
dom=InterpolationDomain[varsa[[1]]][[1]];
Sh\[Infinity]m=Chop[varsa[[shtdInd]]@dom[[2]]];
varshta[[1]]=Chop[varshtalw];
If[verbose, Print["n: ",i,". Redefining lower limit as: ",varshta[[1]]]];
];
Return[Flatten[Join[varshta[[1;;2]],{i}]]];
];
(* First estimates on the value of the shooted variable at eqsRg1[[2]]. There is a goto to correct possible breakups on the equations. We also include provide a 'bracketing' option in case the bracketing on the unknown variable is needed*)
Label[begin];
Sh\[Infinity]=Table[
eqsht={vars[[shtdInd]][eqsRga[[1]]]==(varshta[[i]])};
eqsht=Join[eqs,eqsht];
varsa=vars/.Flatten[NDSolve[eqsht,vars,{Global`r,eqsRga[[1]],eqsRga[[2]]},AccuracyGoal->15,PrecisionGoal->15,WorkingPrecision->30,MaxSteps->Infinity]];
dom=InterpolationDomain[varsa[[1]]][[1]];
(*If[dom[[2]]<eqsRga[[2]],eqsRga[[2]]=Rationalize[0.9*eqsRga[[2]],1];If[verbose,Print[Style[" Breakup found. Changing rfin to ",Red],eqsRga[[2]]] ];Goto[begin];];*)
Sh\[Infinity]m=varsa[[shtdInd]]@dom[[2]]
,{i,3}];
Sh\[Infinity]m=Sh\[Infinity][[3]];
Sh0m=varshta[[3]];
posref=Position[{Abs[Sh\[Infinity][[1]]],Abs[Sh\[Infinity][[2]]]},_?(#==Min[{Abs[Sh\[Infinity][[1]]],Abs[Sh\[Infinity][[2]]]}] &)][[1]];
Sh\[Infinity]m2=Sh\[Infinity][[posref]][[1]];
Sh0m2=varshtdRga[[posref]];
Sh0ref=Flatten[{Sh0m,Sh0m2}];
mean=Mean[Sh0ref];
(* Starts the shooting loop. *)
If[verbose,Print["Output: vars, a, "<>ShtStr<>" variable at \[Infinity], Error"]];
If[verbose,Print["Intermediate prints: iteration, "<>ShtStr<>" at \[Infinity], Error"]];
While[ Abs[Sh\[Infinity]m]>assval &&Abs[1-Sh\[Infinity]m/Sh\[Infinity]m2]>tol && Length@posref>0 && a<= amax,
a=a+1;
If[verbose,Print[{a,Round[mean,10^-16],varsa[[shtdInd]]@dom[[2]],Abs[1-Sh\[Infinity]m/Sh\[Infinity]m2]}]];
eqsht={vars[[shtdInd]][eqsRga[[1]]]==mean};
eqsht=Join[eqs,eqsht];
varsa=vars/.Flatten[NDSolve[eqsht,vars,{Global`r,eqsRga[[1]],eqsRga[[2]]},
AccuracyGoal->15,PrecisionGoal->15,WorkingPrecision->30,MaxSteps->Infinity]];
(*If[dom[[2]]<eqsRga[[2]],eqsRga[[2]]=Rationalize[0.9*eqsRga[[2]],1];If[verbose,Print[Style[" Breakup found. Changing rfin to ",Red],eqsRga[[2]] ]];Goto[begin];];*)
dom=InterpolationDomain[varsa[[1]]][[1]];
Which[assymptotic=="Exponential"&&dom[[2]]==eqsRga[[2]],
Rs=0.95(r/.FindRoot[varsa[[1]]@r,{r,6}]);datab=Table[{r,varsa[[shtdInd]]@r},{r,Rs,dom[[2]],0.01}];
y1=varsa[[shtdInd]]@Rs;y2=(D[varsa[[shtdInd]]@r,r]/.r->Rs);m=(-y1-Rs^2 y2)/(Rs y1);a0= y1;
datfit=NonlinearModelFit[datab,A3 + a0 Exp[-r m],{A3},r];Print[{a0,m,datfit["BestFitParameters"]}]];
posreftest=Quiet@Position[{Abs[Sh\[Infinity]m],Abs[varsa[[shtdInd]]@dom[[2]]]},_?(#==Min[{Abs@Sh\[Infinity]m,Abs[varsa[[shtdInd]]@dom[[2]]]}] &)];
If[Length@posreftest>0,posref=posreftest[[1]],posref={};];
Sh0ref={mean,Sh0ref[[posref]][[1]]};
mean=Mean[Sh0ref];
Sh\[Infinity]m2={Sh\[Infinity]m,Sh\[Infinity]m2}[[posref[[1]]]];
Sh\[Infinity]m=varsa[[shtdInd]]@dom[[2]];
];
If[verbose,Print["Output: vars, a, Shooted variable at \[Infinity]"]];
out=Join[varsa,{a,Round[mean,10^-8],Sh\[Infinity]m}];
Return[out]
]
ComputeEdges[pts_]:=Module[{ptsx,auxvar,auxvar2,nears,i},
ptsx=SortBy[pts,First];
auxvar={};
i=1;
AppendTo[auxvar,{ptsx[[i]]}];
While[i<= Length@ptsx-1,
If[ptsx[[i+1,1]]==ptsx[[i,1]],i=i+1,AppendTo[auxvar,{ptsx[[i]]}];i=i+1]
];
AppendTo[auxvar,{ptsx[[i]]}];
auxvar=Flatten[auxvar,1];
i=Length@ptsx-1;
While[i> 1,
If[ptsx[[i+1,1]]==ptsx[[i,1]],i=i-1,AppendTo[auxvar,ptsx[[i+1]]];i=i-1]
];
AppendTo[auxvar,ptsx[[i]]];
Do[auxvar=AppendTo[auxvar,0.5auxvar[[1]]+0.5auxvar[[-1]]],{i,3}];
auxvar
]
CredibleRegion[data_,level_]:=Module[{datasrt,prob,cumprob,pbound},
(* Last column must be the PDF *)
datasrt=SortBy[data,Last];
prob=TakeColumn[datasrt,-1];
cumprob=Accumulate[prob];
pbound=Quiet@Position[cumprob,_?(#>= (1-level) cumprob[[-1]]& ),1][[1,1]];
ComputeEdges[datasrt[[pbound-1;;-1]]][[All,1;;-2]]
]
LoveNumber[eos_,mtot_]:=Module[{y,r,p,\[Epsilon],m,eqy,\[CapitalGamma],eqsGR,listeos,yy,zz,eqEoS,\[Rho]int,pint,rin,\[Rho]c,pc,eqsIC,rMax,alleqs,alleqsnd,Pr,mr,yr,Rm,Mm,Cc,G,c,m0,rg,P0,\[Rho]0,R0},
G=6.67428 10^-8;
c=2.99792458 10^10;
m0=1.989 10^33;
rg= G m0/c^2;
P0=m0 c^2/rg^3;
\[Rho]0=m0/rg^3;
R0=1/rg^2;
eqEoS={p[r]==(10^EoSFits[eos]/.\[Rho]->Log[10,\[Rho]0 \[Rho][r]])/P0};
listeos=Table[{(EoSFits[eos]/.\[Rho]->x),(From\[Rho]To\[Epsilon]Fits[eos]/.\[Rho]->10^x)},{x,1,16,0.01}];
\[Rho]int=Interpolation[listeos/.{yy_,zz_}->{(10^yy)/P0,(zz)/( \[Rho]0)}];
pint=Interpolation[listeos/.{yy_,zz_}->{(zz)/( \[Rho]0),(10^yy)/P0}];
(* Equations 1-2 GR ToVs. Equation 3 k2 eq. *)
eqsGR={((m[r]+4 \[Pi] r^3 p[r]) (p[r]+\[Rho][r]))/(r^2-2 r m[r])+Derivative[1][p][r]==0,4 \[Pi] r^2 \[Rho][r]-Derivative[1][m][r]==0};
eqy={y'[r]==-(y[r]^2/r)-(r+4\[Pi] r^3 (p[r]-\[Rho][r]))/(r(r-2m[r])) y[r]+(4(m[r] +4\[Pi] r^3 p[r])^2)/(r(r-2m[r]))+6/(r-2m[r])-(4\[Pi] r^2)/(r-2m[r]) (5\[Rho][r]+9p[r]+(\[Rho][r] +p[r])/D[pint@\[Rho][r],\[Rho][r]])};
alleqs=Join[eqsGR,eqy];
(* --------*)
(* Solve the system *)
\[Rho]c=\[Rho]max[[1]]0.76;
rin=10^-5;
pc=eqEoS[[1,2]]/.\[Rho][r]->\[Rho]c;
eqsIC={p[rin]==pc,m[rin]==0,y[rin]==2,WhenEvent[p[r]/pc<10^(-12),rMax=r;"StopIntegration"]};
alleqsnd=alleqs/.\[Rho][r]->\[Rho]int[p[r]]/.p[r]->Max[p[r],0];
{Pr,mr,yr}={p,m,y}/.Flatten[NDSolve[Flatten@Join[alleqsnd,Join[eqsIC]],{p,m,y},{r,rin,100},Method->{"ExplicitRungeKutta","DifferenceOrder"->8},AccuracyGoal->16,PrecisionGoal->13]];
Rm=InterpolationDomain[Pr][[2]];
Mm=(mr@Rm);
Cc=Mm/(Rm);
(2/3((c^2/G)(rg Rm)/(Mm*m0))^5)k2[Cc,yr[Rm]]
];
k2[c_,y_]:=(8 c^5)/5 (1-2c)^2(2+2c(y-1)-y)(2c(6-3y +3c(5y-8))+4c^3(13-11y+c(3y-2)+2c^2(1+y))
+3(1-2c)^2(2-y+2c(y-1))Log[1-2c])^(-1)
(* ::Code::Initialization:: *)
AtomsList[expr_]:=Union@Select[Level[expr,{0,Infinity}],AtomQ];
InterpolationDomain[fun_]:=Module[{min,max},fun[[1]]];
TakeColumn[list1_?ListQ,list2_?ListQ]:=Map[Part[#,list2]&,list1];
TakeColumn[list1_?ListQ,n_?IntegerQ]:=(list1//Transpose)[[n]];
(* ::Code::Initialization:: *)
End[];
EndPackage[];
Options[RiemannTensorDev2]=Join[Options[ChristoffelSymbolDev],{"IndexDown"->False}];
RiemannTensorDev2[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{n,Chr,compile,index,res,perti,simpl,verbose},
index=OptionValue["IndexDown"];perti=OptionValue["PerturbationIndex"];
simpl=OptionValue["SimplifyFunction"];
compile=OptionValue["Compile"];
verbose=OptionValue["Verbose"];
n=Length@xx;
If[verbose,Print["Starting with Christoffel symbols..."]];
Chr=ChristoffelSymbolDev[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl,"Compile"->False];
If[verbose,Print["Christoffel symbols computed. Starting with Riemann..."]];
res=ConstantArray[0,{n,n,n,n}];
If[index,
If[NumericQ[pert], Do[res[[i,k,l,m]]=Sum[g[[i,p]](D[Chr[[p,k,m]],xx[[l]]]-D[Chr[[p,k,l]],xx[[m]]]+Sum[Chr[[p,s,l]]*Chr[[s,k,m]]-Chr[[p,s,m]]*Chr[[s,k,l]],{s,n}]),{p,n}],{i,n},{k,n},{l,n},{m,n}],
Do[res[[i,k,l,m]]=Sum[g[[i,p]](Normal@Series[D[Chr[[p,k,m]],xx[[l]]]-D[Chr[[p,k,l]],xx[[m]]]+Sum[Chr[[p,s,l]]*Chr[[s,k,m]],{s,n}]-Sum[Chr[[p,s,m]]*Chr[[s,k,l]],{s,n}],{pert,0,perti}]),{p,n}],{i,n},{k,n},{l,n},{m,n}]];
If[compile, Do[res[[i,k,l,m]]=If[NumberQ[res[[i,k,l,m]]],res[[i,k,l,m]],Compile[Evaluate@({#,_Real}&/@xx),Evaluate[res[[i,k,l,m]]],CompilationTarget->"C",CompilationOptions->"InlineExternalDefinitions"->True]],{i,n},{k,n},{l,n},{m,n}];];
(* Applying simmetries *)
,
If[NumericQ[pert], Do[res[[i,k,l,m]]=(D[Chr[[i,k,m]],xx[[l]]]-D[Chr[[i,k,l]],xx[[m]]]+Sum[Chr[[i,s,l]]*Chr[[s,k,m]]-Chr[[i,s,m]]*Chr[[s,k,l]],{s,n}]),{i,n},{k,n},{l,n},{m,n}],
Do[res[[i,k,l,m]]=(Normal@Series[D[Chr[[i,k,m]],xx[[l]]]-D[Chr[[i,k,l]],xx[[m]]]+Sum[Chr[[i,s,l]]*Chr[[s,k,m]],{s,n}]-Sum[Chr[[i,s,m]]*Chr[[s,k,l]],{s,n}],{pert,0,perti}]),{i,n},{k,n},{l,n},{m,n}]];
];
If[verbose,Print["...Riemann computed"]];
simpl@res];
Options[RicciTensorDev]=Options[ChristoffelSymbolDev];
RicciTensorDev[xx_,g_,pert_:0,OptionsPattern[]]:=Block[{compile,Rie,res,n,perti,simpl},
perti=OptionValue["PerturbationIndex"];
compile=OptionValue["Compile"];
simpl=OptionValue["SimplifyFunction"];
n=Length@xx;
Rie=RiemannTensorDev[xx,g,pert,"PerturbationIndex"->perti,"SimplifyFunction"->simpl];
res=ConstantArray[0,{n,n}];
If[NumericQ[pert], Do[res[[i,j]]=Sum[Rie[[s,i,s,j]],{s,n}],{i,n},{j,n}],
Do[res[[i,j]]=Normal@Series[Sum[Rie[[s,i,s,j]],{s,n}],{pert,0,perti}],{i,n},{j,i,n}]];
If[compile, Do[res[[i,j]]=If[NumberQ[res[[i,j]]],res[[i,j]],Compile[Evaluate@({#,_Real}&/@xx),Evaluate[res[[i,j]]],CompilationTarget->"C",CompilationOptions->"InlineExternalDefinitions"->True]],{i,n},{j,i,n}];] ;
(* Applying symmetries *)
Do[res[[i+1,j]]=res[[j,i+1]];,{i,n-1},{j,i}];
simpl@res]