Notatnik porównuje Standard Solar Model (SSM) Bahcala i Soneiry z 2005 roku z modelem Eddingtona (politropowy z n=3).

Import["http://th.if.uj.edu.pl/~odrzywolek/aa/Lane-Emden.m"]

zero[3]

slope[3]

zero[1]/Pi-1

bs=Import["http://www.sns.ias.edu/~jnb/SNdata/Export/BS2005/bs05op.dat","Table"];

bs[[1;;24]]

UnitConvert[Quantity[1,"dyn/cm^2"],"Pa"]

bs[[-5;;-1]]

bs[[23]]

temp=bs[[24;;-4,{2,3}]];

temp[[1]]

temp[[-1]]

w2=ListLogPlot[Legended[temp,"BS2005-OP"],Joined->True,PlotStyle->Black]

Tc=temp[[1,-1]]

Teff=SetPrecision[QuantityMagnitude[%,"Kelvins"],$MachinePrecision]

dens=bs[[24;;-4,{2,4}]];

w1=ListLogPlot[Legended[dens,"BS2005"],Joined->True,PlotStyle->Black]

ρC=dens[[1,2]]

R=dens[[-1,1]]

R=1

Aczkolwiek nie jest to zupełnie oczywiste, rachunek pokazuje, że rozkład temperatury wewnątrz gwiazdy w modelu Eddingtona jest dany po prostu funkcją Lane-Emdena z n=3.

w2Eddi=LogPlot[{Legended[Tc w[3][r zero[3]],"Model\nEddingtona"],Legended[Teff,"(Subscript[T, eff])=5772 K"]},{r,0,1},PlotStyle->{Red,Yellow},PlotRange->{Teff,All}]

Show[w2,w2Eddi,PlotRange->All,Frame->True,Axes->None,PlotLabel->"Model Eddingtona vs Bahcall-Soneira (2005)",FrameLabel->{"r/(Subscript[R, ⊙])","T(r) [K]"},ImageSize->600,BaseStyle->FontSize->18]

Manipulate[Show[w1,LogPlot[ρC w[n][zero[n]r]^n,{r,0,1}]],{{n,3},0,5,0.1}]

Show[LogPlot[Evaluate@Table[ρC w[n][zero[n]r]^n,{n,3,3+1/2,1/4}],{r,0,1},PlotStyle->{Red,Green,Blue},PlotLegends->{"n=3","n=3.25","n=3.5","n=0.6"}],w1,AxesLabel->{"r/(Subscript[R, ⊙]) [km]","ρ [g/(cm^3)]"},BaseStyle->FontSize->18,ImageSize->600,PlotLabel->"Modele politropowe vs Bahcall-Soneira (2005)"]

Poniżej przeanalizujmy równanie stanu materii użyte w Standarowym Modelu Słońca (SSM BP-2005_.

EOS={bs[[24;;-4,4]],0.1bs[[24;;-4,5]]}//Transpose;

Na wykresie podwójnie logarytmicznym równanie stanu P(ρ) ma postać liniii łamanej. Każdy prostoliniowy odcinek to prawo potęgowe, P = K (ρ^γ).

ListLogLogPlot[EOS]

Length[EOS]

Clear[K]

Ma początek spróbujemy dopasować jedno prawo potęgowe do całego równania stanu.

fitLog10=FindFit[Log10@EOS,lgK+4/3 lgρ,{{lgK,14}},lgρ]

fitAll=FindFit[EOS,K ρ^γ,{{K,10^14},{γ,1.4}},ρ]

fitAllloglog=FindFit[Log10@EOS,lgK +γ lgρ,{{lgK,14},{γ,1.4}},lgρ]

Jak widać, fit w postaci logarytmicznej rzeczywiście daje równanie stanu z wykładnikiem bliskim 4/3 !

{{"","K","γ"},{"lgP-lgρ", 10^lgK,γ}/.fitAllloglog,{"P-ρ",K,γ}/.fitAll}//TableForm (* UWAGA: fitowanie praw potęgowych daje bardzo różne (!) wyniki w postali P-ρ vs lgP - lgρ !   *)

Teraz przystępujemy do zbadania równania stanu dla niskich i dużych gęstości osobno.

fitLOWdens=FindFit[Log10@EOS[[900;;]],lgK +γ lgρ,{{lgK,14},{γ,1.4}},lgρ]

1/(γ-1)/.% (* Wynik bliski n=3/2, γ=5/3 *)

fitHIGHdens=FindFit[Log10@EOS[[1;;900]],lgK +γ lgρ,{{lgK,14},{γ,1.4}},lgρ]

1/(γ-1)/.% (* Wynik bliski n=4, γ=1.25 *)

1+1/(3/2)

1/(5/4-1)

Show[
ListLogLogPlot[Legended[EOS[[1;;]],"True EOS"]],
LogLogPlot[
{
Legended[10^(lgK +4/3Log10[ρ])/.fitLog10,"P = K (ρ^(4/3))"],
Legended[10^(lgK +γ Log10[ρ])/.fitAllloglog,"P = K (ρ^γ)"],
Legended[10^(lgK +γ Log10[ρ])/.fitLOWdens,"Low density,\n γ=1.66≃5/3\nn=3/2"],
Legended[10^(lgK +γ Log10[ρ])/.fitHIGHdens,"High density,\n γ=1.23≃5/4\nn=4"]
},
{ρ,10^-3,200},PlotStyle->{Red,{Blue,Dashed},{Black,Dotted},{Black,DotDashed}}]
,PlotRange->All,BaseStyle->FontSize->18,ImageSize->600,AxesLabel->{"ρ [g/(cm^3)]","P [Pa]"},PlotLabel->"Równanie stanu P(ρ) wewnątrz Słońca"]

Skuteczność modelu Eddingtona można częściowo zrozumieć patrząc przez pryzmat wynikającej z niego średniej drogi swobodnej fotonu.

Lsun=SetPrecision[QuantityMagnitude[%,"W"],$MachinePrecision]

Rsun=SetPrecision[QuantityMagnitude[%,"meters"],$MachinePrecision]

Msun=SetPrecision[QuantityMagnitude[%,"kg"],$MachinePrecision]

β=10^-3

c=299792458

G=6.67*10^-11

mass=bs[[24;;-4,1]];

ρ=bs[[24;;-4,4]];

Wzór na średnią drogę swobodną w modelu Eddingtona można wyprowadzić w poniższy sposób:

L(r)/(4 Pi r^2) =- 1/3 Lγ c D[a T[r]^4, r] =- 1/3 Lγ c D[Prad,r] =- 1/3 Lγ c D[(1-β) P,r] = -1/3 (1-β) Lγ c D[P,r] = 1/3 (1-β) Lγ c G m(r) ρ(r)/r^2

Lγ=3Lsun/(4 Pi β c G)/(ρ*1000)/(mass*Msun);

radius=bs[[24;;-4,2]];

QuantityMagnitude[Quantity[1,"ThomsonCrossSection"],"m^2"]

Lthompson=QuantityMagnitude[1/Quantity[1,"ThomsonCrossSection"]/(Quantity[dens[[1,2]],"g/cm^3"]*0.7/Quantity[1,"AtomicMassUnit"]),"Meters"]

lgamma3=ListLogPlot[Lγ*1000,Frame->True,FrameLabel->{"Numer komórki","Średnia droga swobodna fotonu [mm]"},BaseStyle->FontSize->14,ImageSize->600,GridLines->{None,{Lthompson*1000}}]

lgamma2=ListLogPlot[Transpose[{mass,Lγ*1000}],Frame->True,FrameLabel->{"m/(Subscript[M, ⊙])","Średnia droga swobodna fotonu [mm]"},BaseStyle->FontSize->14,ImageSize->600,GridLines->{None,{Lthompson*1000}}]

lgamma1=ListLogPlot[Transpose[{radius,Lγ*1000}],Frame->True,FrameLabel->{"r/(Subscript[R, ⊙])","Średnia droga swobodna fotonu [mm]"},BaseStyle->FontSize->18,ImageSize->600,GridLines->{None,{Lthompson*1000}}]