Skip to content
Snippets Groups Projects
Select Git revision
  • 4966dc2c23db379cae6df76f427147ad048ad371
  • master default protected
  • Binary
  • add-version-information
  • os-path-join
  • develop-GA
  • timeFstatmap
  • add-higher-spindown-components
  • develop-DK
  • adds-header-to-grid-search
  • v1.3
  • v1.2
  • v1.1.2
  • v1.1.0
  • v1.0.1
15 results

plot_data.py

Blame
  • Forked from Gregory Ashton / PyFstat
    Source project has a limited visibility.
    GRTensor.m 69.95 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])";
    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.";
    
    
    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[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
    ];
    
    
    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]