SetDirectory[NotebookDirectory[]]

Get["Figures_style.wl"]

u=Import["U.wl"] (* Elastic function *)

(* HERE we select central overdensity for plot !!!!  *)

cases={{"a",1.001},{"b",1.01},{"c",1.1},{"d",2}};

case="c" (* Select cases: a,b,c,d *)

(* Data for Fe from CHRISTINE J.WU et al.PHYSICAL REVIEW B 108,014102 (2023) Table V, faza bcc (α) *)

ρ0= Quantity[7.927,"g/cc"]//UnitConvert//QuantityMagnitude(* g/cc *)

K = Quantity[164.8,"GPa"]//UnitConvert//QuantityMagnitude (* GPa *)

Kprim = 5.55

(* Liquid iron *)

(* 

ρ0 = Quantity[8.027,"g/cc"]//UnitConvert//QuantityMagnitude
K=Quantity[147.6,"GPa"]//UnitConvert//QuantityMagnitude
Kprim=7.82;

*)

μ=ρ/ρ0-1

P[ρ_]:=Evaluate[K μ-1/2 K (1-Kprim)μ^2] (* Quadratic reference EOS *)

Pgeneralized[ρ_]:=Evaluate[K/Kprim((ρ/ρ0)^Kprim-1)] (* Generalized elastic EOS *)

Plinearized[ρ_]:=Evaluate[K μ] (* Linearized EOS = Quadratic/Generalized for Kprim=1 *)

Pelastic[ρ_]:=Evaluate[K Log[ρ/ρ0]]

P[8500]

LogPlot[{P[ρ],Pgeneralized[ρ],Plinearized[ρ],Pelastic[ρ]},{ρ,ρ0,2ρ0},PlotLegends->"Expressions"]

G=Quantity[1,"GravitationalConstant"]//UnitConvert//QuantityMagnitude

ρC=Switch[case,"a",cases[[1]],"b",cases[[2]],"c",cases[[3]],"d",cases[[4]]][[-1]]*ρ0

ηC=Switch[case,"a",cases[[1]],"b",cases[[2]],"c",cases[[3]],"d",cases[[4]]][[-1]]

(* Fig5, panel a  *)

(*ρC=1.01*ρ0//UnitConvert//QuantityMagnitude *)
(* Fig5, panel b *)

(*ρC=1.1*ρ0//UnitConvert//QuantityMagnitude *)
(* Fig5, panel c *)

(*ρC=2*ρ0//UnitConvert//QuantityMagnitude *)
(* Fig5, panel d *)

ϵ=Quantity[2^-13,"Meters"]//UnitConvert//QuantityMagnitude

Rmax=Quantity[10000,"Kilometers"]//UnitConvert//QuantityMagnitude

rozwNUMref=NDSolve[{D[P[ρ[r]],r]==-G m[r]ρ[r]/r^2,m'[r]==4 Pi r^2 ρ[r],ρ[ϵ]==ρC,m[ϵ]==4/3Pi ρC ϵ^3,WhenEvent[ρ[r]<ρ0,"StopIntegration"]},{ρ,m},{r,ϵ/128,Rmax}][[1]]

rozwNUMgeneralized=NDSolve[{D[Pgeneralized[ρ[r]],r]==-G m[r]ρ[r]/r^2,m'[r]==4 Pi r^2 ρ[r],ρ[ϵ]==ρC,m[ϵ]==4/3Pi ρC ϵ^3,WhenEvent[ρ[r]<ρ0,"StopIntegration"]},{ρ,m},{r,ϵ/128,Rmax}][[1]]

rozwNUMlinearized=NDSolve[{D[Plinearized[ρ[r]],r]==-G m[r]ρ[r]/r^2,m'[r]==4 Pi r^2 ρ[r],ρ[ϵ]==ρC,m[ϵ]==4/3Pi ρC ϵ^3,WhenEvent[ρ[r]<ρ0,"StopIntegration"]},{ρ,m},{r,ϵ/128,Rmax}][[1]]

{dens,mass}={ρ,m}/.rozwNUMref;

{densGen,massGen}={ρ,m}/.rozwNUMgeneralized;

{densLin,massLin}={ρ,m}/.rozwNUMlinearized;

dens["Properties"]

R=(dens["Region"]//RegionBounds)[[1,-1]]

mass[R]

mass[R]/(Quantity[1,"MoonMass"]//UnitConvert//QuantityMagnitude)

mass[R]/(Quantity[1,"CeresMass"]//UnitConvert//QuantityMagnitude)

λ=Sqrt[K/(4Pi G ρC^2)]

R0=Sqrt[K/(4Pi G ρ0^2)]

densGen[R]

R/λ

kmTBL=FindDivisions[{0,R/1000},10]

kmTicks=Table[{1000*km,ToString[km]},{km,kmTBL}]

xTBL=FindDivisions[{0,R/R0},10]

λTicks=Table[{R0 x,N@x},{x,xTBL}]//Rest

ρTBL=FindDivisions[{ρ0,ρC},10]

digits=Switch[case,"a",4,"b",3,"c",2,"d",1]-1

ρTicks=Table[{ρ,NumberForm[ρ/1000.,{4,digits}]},{ρ,ρTBL}]

Log10[μ/.ρ->ρC]

μTicks=Join[{{ρ0,"-∞"}},Table[{ρ0(1+10^lgμ),NumberForm[lgμ,{3,2}]},{lgμ,Log10[μ/.ρ->ρC]-1.25,0,0.25}]]

profiles={
dens[r],
ρC(1-r^2/R^2)+ρ0 r^2/R^2,
ρC u[r/λ],
densGen[r],
densLin[r]
};

rTBL={R/2};

dataForPlot=Table[profiles,{r,rTBL}]//Transpose

curves=Table[{rTBL,dataForPlot[[k]]}//Transpose,{k,1,Length[dataForPlot]}]

markers={"▫","○", "□", "◇", "△", "▽"};

styles={
{Black,AbsoluteThickness[thickness]},
{Blue,Dashed,AbsoluteThickness[thickness]},
{Red,DotDashed,AbsoluteThickness[thickness]},
{Purple,Dashed,AbsoluteThickness[thickness]}, 
{Green,Dashed,AbsoluteThickness[thickness]}
};

(*legend=If[case=="a",Placed[LineLegend[{"Reference [Quad EOS]","Quadratic profile","Elastic model","Generalized elastic","Linear EOS"}], {Left,Bottom}],False,False]*)

legend=Placed[LineLegend[{"Reference [Quad EOS]","Quadratic profile","Elastic model","Generalized elastic","Linear EOS"}], {Left,Bottom}];

fig5symbols=ListPlot[
curves,
PlotMarkers->markers,
Joined->True,
PlotStyle->styles,
Frame->True,
FrameTicks->{{ρTicks,μTicks},{kmTicks,λTicks}},
PlotLegends->legend]

fig5lines=Plot[profiles,
{r,0,3/2R},
PlotRange->{{0,9/8R},{ρ0,ρC}},
AxesOrigin->{0,0},
PlotStyle->styles,
Frame->True,
FrameLabel->{{"ρ [g/(cm^3)]", "(Subscript[log, 10])μ"}, {"r [km]", "r/(Subscript[R, 0])"}},
Epilog->{Text["(Subscript[ρ, c])="<>ToString[ηC]<>" (Subscript[ρ, 0])\n"<>"R="<>ToString[EngineeringForm[R/1000,{4,1}],StandardForm]<>" km\nM="<>ToString[ScientificForm[mass[R],{3,2}],StandardForm]<>" kg",{8/9R,7(ρC-ρ0)/8+ρ0}]},
ImageSize->450,
FrameTicks->{{ρTicks,μTicks},{kmTicks,λTicks}},
FrameStyle->{{Black,AbsoluteThickness[1]},{Black,AbsoluteThickness[1]},{Black,AbsoluteThickness[1]},{Black,AbsoluteThickness[1]}}
]

fig5=Show[fig5lines,fig5symbols]

case

Export["Fig5_"<>case<>".pdf",fig5]