Celem notatnika jest zademonstrowanie, że równanie transportu promienistego energii w gwiazdach jest de facto równaniem dyfuzji energii gazu fotonowego w sferycznej symetrii. Pozwala to na uzyskanie bardzo prostej, łatwej do zapamiętania i fizycznie umotywowanej formy tegoż równania:

Fγ = L[r]/(4 Pi r^2) = -Dc D[aT[r]^4,r]

gdzie  Fγ to strumień fotonów przez sferę o promieniu r, L[r]  - strumień energii (jasność gwiazdy), Dc - współczynnik dyfuzji, ε=aT^4 - gęstość energii gazu fotonowego.

Kluczowy i trudny problem to powiązanie współczynnika Dc z mikroskopowymi wielkościami (obliczanymi z teorii lub mierzonymi w laboratoriach) jak średnia droga swobodna fotonu Lγ, powiązana z przekrojem czynnym σ i gęstością cząstek n poprzez Lγ=1/(σ n) LUB z nieprzeźroczystością materii κ oraz gęstością ρ za pomocą Lγ=1/(κ ρ). Poniżej problem ten rozwiązany jest poprzez porównanie procesu błądzenia przypadkowego fotonów z analitycznym rozwiązaniem r. dyfuzji.

Zaczniemy od sprawdzenia/przypomnienia analitycznego rozwiązania równania dyfuzji we współrzędnych sferycznych.  Rozwiązanie to ma postać rozpływającego się impulsu gaussowskiego:

((exp(-((x^2)/((4 Dc) t))))/(((4 π Dc t))^(3/2)));

(* Rozwiązanie analityczne równania dyfuzji; współczynnik dyfuzji oznaczam jako Dc ponieważ D w Mathematice to różniczkowanie. *)

D[u[x,t],t]-Dc D[u[x,t],x,x]-Dc 2/x D[u[x,t],x]/.u->Function[{x,t},1/(4 Pi Dc t)^(3/2) Exp[-x^2/4/Dc/t]] (*We wsp. sferycznych *)

Simplify[%]

D[u[x,y,z,t],t]- Dc D[u[x,y,z,t],x,x]- Dc D[u[x,y,z,t],y,y]- Dc D[u[x,y,z,t],z,z]/.u->Function[{x,y,z,t},1/Sqrt[4Pi Dc  t]^3Exp[-(x^2+y^2+z^2)/4/Dc/t]] (* We wsp. kartezjańskich *)

Simplify[%]

Dc=1

Manipulate[Plot[1/(4 Pi Dc t)^(3/2) Exp[-x^2/4/Dc/t],{x,-3,3},PlotRange->{0,32}],{{t,0.01},0.005,1}]

(*DSolve[-1+σ[t] σ'[t]==0,σ,t]*)

Przystępujemy teraz do symulacji Monte Carlo błądzenia przypadkowego w 3D. Potrzebujemu m.in. losowego kierunku (losowego wersora) .

(* Losowy kierunek/wektor jednostkowy *)

RandomDirection=Function[{},Module[{tmp},
tmp=RandomReal[{-1,1},3];
While[Norm[tmp]>1,tmp=RandomReal[{-1,1},3];];
tmp/Norm[tmp]
]
];

ListPointPlot3D[Table[RandomDirection[],{i,10000}],BoxRatios->{1,1}]

Proces wyobrażamy sobie następująco: foton jest emitowany z r=0, porusza się prostoliniowo aż do momentu zderzenia np. z elektronem, po czym zmienia kierunek i dokonuje kolejnego przeskoku. Proces powtarzamy wielokrotnie.

(* Przykład błądzenia przypadkowego pojedynczego fotonu *)

Graphics3D[Line[Accumulate[Table[RandomVariate[ExponentialDistribution[1]]RandomDirection[],{182}]]],Axes->True,Boxed->False,AxesOrigin->{0,0,0}]

Ponieważ naszym celem jest oszacowanie gęstości fotonów, potrzebujemy ,,dużego’’ ich zbioru.

(* Symulacja błądzenia przypadkowego *)

num=128 (* Liczba fotonów *)

steps=256 (* Liczba skoków *)

Mean@ExponentialDistribution[1/L]

StandardDeviation@ExponentialDistribution[1/L]

L=1

Symulację powtarzamy w celach ilustracyjnych przyjmując różne sposoby losowania długości skoków. Fizycznie motywowany jest rozkład wykładniczy, jako że prawdopodobieństwo, że foton w nic nie trafi spada w stały sposób. Inne przykłady to: skoki o stałej długości, skoki po siatce sześciennej, skoki o jednorodnym rozkładzie i jako ciekawostka rozkład ,,ciężkoogonowy’’.

rnd2=Table[Prepend[{0,0,0}]@Accumulate[Table[L RandomDirection[],{steps}]],{num}];

rnd1=Table[Prepend[{0,0,0}]@Accumulate[Table[RandomChoice[L{{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}}],{steps}]],{num}];

rnd4=Table[Prepend[{0,0,0}]@Accumulate[Table[RandomVariate[ExponentialDistribution[1/L]]RandomDirection[],{steps}]],{num}];

rnd3=Table[Prepend[{0,0,0}]@Accumulate[Table[RandomVariate[UniformDistribution[{0,2L}]]RandomDirection[],{steps}]],{num}];

rnd5=Table[Prepend[{0,0,0}]@Accumulate[Table[RandomVariate[LevyDistribution[L,1]]RandomDirection[],{steps}]],{num}];

w1=Graphics3D[{Thin,Line[rnd1[[;;32,;;-1]]]},Axes->True,Boxed->False,AxesOrigin->{0,0,0},PlotRange->{{-1,1},{-1,1},{-1,1}}8];

w2=Graphics3D[{Thin,Line[rnd2[[;;32,;;-1]]]},Axes->True,Boxed->False,AxesOrigin->{0,0,0},PlotRange->{{-1,1},{-1,1},{-1,1}}8];

w3=Graphics3D[{Thin,Line[rnd3[[;;32,;;-1]]]},Axes->True,Boxed->False,AxesOrigin->{0,0,0},PlotRange->{{-1,1},{-1,1},{-1,1}}8];

w4=Graphics3D[{Thin,Line[rnd4[[;;32,;;-1]]]},Axes->True,Boxed->False,AxesOrigin->{0,0,0},PlotRange->{{-1,1},{-1,1},{-1,1}}8];

GraphicsGrid[{{w1,w2},{w3,w4}}]

rnd={rnd1,rnd2,rnd3,rnd4};

Manipulate[Graphics3D[{Thin,Line[rnd[[type]][[;;part,;;k]]]},Axes->True,Boxed->False,AxesOrigin->{0,0,0},PlotRange->{{-1,1},{-1,1},{-1,1}}12],{{k,32},1,steps,1},{{part,32},1,num,1},{type,{1,2,3,4}}]

Manipulate[
SetOptions[Graphics3D,Axes->True,Boxed->False,AxesOrigin->{0,0,0},PlotRange->{{-1,1},{-1,1},{-1,1}}60,ImageSize->300];
Graphics3D[{PointSize[0.01],Point[rnd3[[;;part,k]]]}]
,{{k,8},1,steps,1},
{{part,32},1,num,1}]

(*Export["RandomWalk_delta_step8.pdf",%]*)

(*Export["RandomWalk_grid_step8.pdf",%]*)

(*Export["RandomWalk_exp_step8.pdf",%]*)

(*Export["RandomWalk_uniform_step8.pdf",%]*)

Teraz przystępujemy do przeliczenia dyskretnego procesu krokowego na proces zachodzący w czasie. Jest to o tyle łatwe, iż fotony poruszają się ze stałą prędkością równą c.

c=1

diff=1/3 L c

(* Porówanie dyfuzji typu błądzenia przypadkowego ze wzorem analitycznym *)

Show[Plot[4Pi x^2/Sqrt[4 Pi diff steps L]^3 Exp[-x^2/(4diff steps L)],{x,0,4Sqrt@steps L},PlotStyle->{Black,DotDashed},PlotRange->{0,All}],
SmoothHistogram[{Norm/@rnd1[[;;,-1]],Norm/@rnd2[[;;,-1]],Norm/@rnd3[[;;,-1]],Norm/@rnd4[[;;,-1]]},Automatic,"PDF",PlotRange->All,PlotLegends->{"Po krawędziach sześcianu","Stała długość skoku","Jednorodny","Wykładniczy"}]]

(Legended[, (LineLegend["Po krawędziach sześcianu", "Stała długość skoku", "Jednorodny", "Wykładniczy"])])