This file generates the equation of state and MR plots from  tabulated particle numbers (use the ones from tabluated.nb or create your own using EoS.nb). The parameter muMin gives the minimal  quark number chemical potential (in GeV) for which the absolutely stable quark matter is chosen to exist.

-Pick a value for minimal quark chemical potential, where quark matter exist :

muMin=0.9/3.

-Choose quark equation of state :

Central values

No CSC

Ntable=LambdaMSbar0378ms0010Δ0X1;

Ntable=LambdaMSbar0378ms0010Δ0X2;

Ntable=LambdaMSbar0378ms0010Δ0X4;

With CSC

Ntable=LambdaMSbar0378ms0010Δ01X4;

Ntable=LambdaMSbar0378ms0010Δ01X2;

Ntable=LambdaMSbar0378ms0010Δ01X1

Favoured values

No CSC

Ntable=LambdaMSbar0346ms0007Δ0X1;

Ntable=LambdaMSbar0346ms0007Δ0X2;

Ntable=LambdaMSbar0346ms0007Δ0X4;

With CSC

Ntable=LambdaMSbar0346ms0007Δ01X4;

Ntable=LambdaMSbar0346ms0007Δ01X2;

Ntable=LambdaMSbar0346ms0007Δ01X1

Run This :

Quark EoS

gevtoinvfm=5.06842372022301;

n0f3[mu_]:=(((mu^3) Nf)/(π^2))/.Nf-> 3

Neutraln2total=Transpose[{Ntable[[All,1]],(Ntable[[All,2]]+Ntable[[All,3]]+Ntable[[All,4]]+Ntable[[All,5]])/n0f3[Ntable[[All,1]]]}];
Neutraln2quark=Transpose[{Ntable[[All,1]],(Ntable[[All,2]]+Ntable[[All,3]]+Ntable[[All,4]])/n0f3[Ntable[[All,1]]]}];

neutraln2totalI:=Interpolation[Neutraln2total]
neutraln2quark:=Interpolation[Neutraln2quark]
Ns:=Interpolation[Transpose[{Ntable[[All,1]],(Ntable[[All,2]])/n0f3[Ntable[[All,1]]]}]]
Nu:=Interpolation[Transpose[{Ntable[[All,1]],(Ntable[[All,3]])/n0f3[Ntable[[All,1]]]}]]
Nd:=Interpolation[Transpose[{Ntable[[All,1]],(Ntable[[All,4]])/n0f3[Ntable[[All,1]]]}]]
Ne:=Interpolation[Transpose[{Ntable[[All,1]],(Ntable[[All,5]])/n0f3[Ntable[[All,1]]]}]]
Mue[mu_]:=((3π^2)Ne[mu]*n0f3[mu])^(1/3)
dMuedmu[mu_]:=Evaluate[D[Mue[mu],mu]]
MuQ[mu_]:=mu-1/3 Mue[mu]
MuS[muq_]:=mu/.FindRoot[mu-1/3Mue[mu]==muq,{mu,muq}]

neutralp2total[mu_,B_]:=If[mu<MuS[muMin],-B,NIntegrate[Ns[mup] n0f3[mup]+Nd[mup] n0f3[mup]+(1-dMuedmu[mup]) Nu[mup] n0f3[mup]+dMuedmu[mup] Ne[mup] n0f3[mup],{mup,MuS[muMin],mu}]-B]

neutrale2total[mu_,B_]:=-neutralp2total[mu,B]+(mu Ns[mu]+mu Nd[mu]+(mu-Mue[mu])Nu[mu]+Mue[mu]Ne[mu])*n0f3[mu]

strange = Table[{neutralp2total[mu, 0], neutrale2total[mu,0]}, {mu,muMin, .6, 0.004}];
strangeEoS={};
Do[If[strange[[i,1]]>10^(-7),strangeEoS=Append[strangeEoS,strange[[i]]]],{i,1,Length[strange]}]

Mass - Radius curves from Tolmann - Oppenheimer - Volkov

G=1/(1.2211* 10^(19))^2;

GeVtokg=1.7827 *10^(-27);

ooGeVtokm=1.9733 10^(-16)/1000;

newNUQSPress[EOST_,LAST_]:=Module[{rho,rhou,M,dp,p,r,dr,dr2,cnt},
(*LAST=Length[EOST];*)
(*LAST=30;*)
p=EOST[[LAST,1]];
rho=EOST[[LAST,2]];
dp=(EOST[[LAST-1,1]]-EOST[[LAST,1]]);
dr=Sqrt[-dp/4/Pi/G/(p+rho)/(p+rho/3)];
r=dr;
M=4*Pi*dr^3/3*rho;
cnt=0;
While[LAST-cnt>2,cnt=cnt+1;
p=EOST[[LAST-cnt,1]];
rho=EOST[[LAST-cnt,2]];
rhou=EOST[[LAST-cnt-1,2]];
dp=(EOST[[LAST-1-cnt,1]]-EOST[[LAST-cnt,1]]);
dr=-dp*1/N[G*(p+rho)*(M+4*Pi*p*r^3)/(r-2*G*M)/r];p=EOST[[LAST-cnt-1,1]];
dr2=-dp*1/N[G*(p+rhou)*(M+4 Pi rhou (r+dr)^2 dr+4*Pi*p*(r+dr)^3)/(r+dr-2*G*(M+4 rhou Pi (r+dr)^2 dr))/(r+dr)];
dr=(dr+dr2)/2;
M=N[M+4*Pi*(rho+rhou)/2*(r+dr/2)^2*dr];
r=N[r+dr];];{r*ooGeVtokm,M*GeVtokg/1.989/10^(30)}]

MassvsRadius=Table[newNUQSPress[Sort[strangeEoS],n],{n,2,Length[strangeEoS]-1,2}];

Results

Some information:

stop=1;imin=0;Do[If[stop>0,If[Ns[Ntable[[i,1]]]>0,stop=0;imin=i;]],{i,1,Length[Ntable]}];mumMin=Ntable[[imin,1]];Print["Minimal mu_B where 3flavor SQM exists:",mumMin*3," GeV compared to your choice: ",muMin*3," GeV"];If[mumMin<muMin,Print["You may consider decreasing muMin"],Print["WARNING: Your choice contradicts experiment (no 2 flavor quark matter at zero pressure). It is recommended that you increase muMin!!!"]]

Mass Radius relationship

ListPlot[{MassvsRadius}, 
AxesLabel-> {"R[km]","M / (Subscript[M, sol])"},PlotRange-> {{0,20},{0,3}}]

Tabulate EoS

Gives the EoS in tabulated format EoSlist[[n]]={Baryon mu[GeV], Pressure[GeV/fm^3], Energy Dens.[GeV/fm^3], Baryon Dens.[1/fm^3], Electron fraction}

EoSlist={};
Do[Quiet[If[neutralp2total[Ntable[[i,1]],0]>0,EoSlist=Append[EoSlist,{3*MuQ[Ntable[[i,1]]],neutralp2total[Ntable[[i,1]],0]*gevtoinvfm^3,neutrale2total[Ntable[[i,1]],0]*gevtoinvfm^3,neutraln2quark[Ntable[[i,1]]*gevtoinvfm^3*n0f3[Ntable[[i,1]]]]/3,Ne[Ntable[[i,1]]]/(neutraln2quark[Ntable[[i,1]]]/3)}]]],{i,1,Length[Ntable]-1}]

EoSlist