Polytropic Equation of State implementation
In[]:=
ClearAll[EoS];​​​​EoS[asoc_][key_]:=asoc[key]/;MemberQ[Keys[asoc],key]​​
In[]:=
EoSPoly[k_,Γ_]:=Module{ρ},​​ ρ[p_?NumericQ]:=
1
Γ
p
k
+
p
Γ-1
;(*Returnsenergy-densitygivenpressure*)​​ EoS<|​​ "name"->"PolyEos("<>ToString[k]<>"|"<>ToString[Γ]<>")",​​ "k"->k,"γ"->Γ,​​ "ρ"->ρ​​ |>​​​​
TOVsolver for TOV equations in usual form
In[]:=
ClearAll[TOVSolver`solution];​​​​TOVSolver`solution[asoc_][key_]:=asoc[key]/;MemberQ[Keys[asoc],key]​​
In[]:=
Options[TOVSolver]=​​"rmin"->1.*
-4
10
, ​​ (*NDSolveParams*)​​ AccuracyGoal12,​​ PrecisionGoal12,​​ MaxSteps->
6
10
​​;​​​​TOVSolver[pc_,eos_]:=Module{eqs,p,ρ,ρc,m,mc,r,rc,sol},​​rc=OptionValue[TOVSolver,"rmin"];​​ρ=eos["ρ"];​​ρc=ρ[pc];​​mc=
4
3
πρc
2
rc
;​​​​eqs=m'[r]==4πρ[r]
2
r
,p'[r]==
(ρ[r]+p[r])(4πp[r]
3
r
+m[r])
(r+0.001)(r+0.001-2m[r])
,m[0]==mc,p[0]==pc;​​​​sol=First@NDSolve[eqs,{m,p},{r,rc,20},​​Method->{"StiffnessSwitching"},​​AccuracyGoalOptionValue[TOVSolver,AccuracyGoal],​​PrecisionGoalOptionValue[TOVSolver,PrecisionGoal],​​MaxStepsOptionValue[TOVSolver,MaxSteps]​​];​​{m,p}={m,p}/.sol;​​​​Module[{ri,pi,ρi,mi},​​ri={rc}~Join~m["Coordinates"][[1]];​​pi={pc}~Join~p["ValuesOnGrid"];​​mi={mc}~Join~m["ValuesOnGrid"];​​ρi=Map[ρ[#]&,pi];​​TOVSolver`solution[<|​​"eos"->eos,​​"pc"->pc,"ρc"->ρc,"mc"->mc,​​"R"->Last@ri,"M"->Last@mi​​|>]​​]​​​​
Mass-Radius curves for polytropes
In[]:=
EOSValues=EoSPoly[0.05,N@#]&/@Table
8
6
+
i
6
,{i,0,4};
In[]:=
data=ParallelTable[Table[Quiet@TOVSolver[
i
10
,e],{i,-5,1,1/20}],{e,EOSValues},DistributedContextsAll,Method"FinestGrained"];
In[]:=
mri=Map[d|->Map[{#["R"],#["M"]}&,d],data];
In[]:=
ListLogLinearPlot[mri]
Out[]=