(* Rozwiązanie analityczne równania dyfuzji *)

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]]

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]]

Simplify[%]

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

(* 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,1000}],BoxRatios->{1,1}]

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

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

(* Symulacja błądzenia przypadkowego *)

num=256 (* Liczba fotonów *)

steps=512 (* Liczba skoków *)

Mean@ExponentialDistribution[1/L]

StandardDeviation@ExponentialDistribution[1/L]

L=1

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}];

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}}20],{{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",%]*)

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->Red,PlotRange->{0,All}],SmoothHistogram[Norm/@rnd1[[;;,-1]],Automatic,"PDF",PlotRange->All],SmoothHistogram[Norm/@rnd2[[;;,-1]],Automatic,"PDF"],SmoothHistogram[Norm/@rnd3[[;;,-1]],Automatic,"PDF",PlotStyle->Blue],SmoothHistogram[Norm/@rnd4[[;;,-1]],Automatic,"PDF",PlotStyle->Green]]

2Sqrt[2]+1+Sqrt[17]//N