Free particle

Afree=Compile[{{l,_Integer},{kred,_Real},{r,_Real}},Module[{kr,ans},
kr=2 π kred r;
SphericalBesselJ[l,kr]
]]

ψfreel=Compile[{{l,_Integer},{kred,_Real},{x,_Real},{y,_Real}},Module[{r,cosθ,kr},
r=Sqrt[(x^2)+(y^2)];
If[r==0.0,cosθ=1.0,cosθ=x/r];
kr=2 π kred r;
(2l+1)I^l SphericalBesselJ[l,kr] LegendreP[l,cosθ]
]]

DensityPlot[Re[ψfreel[2,1.0,x,y]],{x,-5,5},{y,-5,5},PlotPoints->100,PlotRange->All]

ψfree=Compile[{{kred,_Real},{x,_Real},{y,_Real},{lmax,_Integer}},Module[{r,cosθ,kr},
r=Sqrt[(x^2)+(y^2)];
If[r==0.0,cosθ=1.0,cosθ=x/r];
kr=2 π kred r;
Sum[(2l+1)I^l SphericalBesselJ[l,kr] LegendreP[l,cosθ],{l,0,lmax}]
]]

DensityPlot[Re[ψfree[1.0,x,y,40]],{x,-5,5},{y,-5,5},PlotPoints->100,PlotRange->All]

Hard sphere scattering

Scattering of scalar waves by an impenetrable sphere of radius a.

Lengths in units of a, wave numbers in units of 1/a.

kred=k a/(2π)    k a = 2 π  kred

(* not needed for the plots of the waves *)(* δ[l_,ka_]:=ArcTan[-SphericalBesselY[l,ka],-SphericalBesselJ[l,ka]] *)

A=Compile[{{l,_Integer},{kred,_Real},{r,_Real}},Module[{ka,kr,H1,ans},
ans=0.0+I 0.0;
If[r>1.0,
ka=2 π kred;
kr=ka r;
(* H1 = SphericalHankelH1[l,ka] *)
H1=SphericalBesselJ[l,ka]+I SphericalBesselY[l,ka]; ans=I(SphericalBesselY[l,ka]SphericalBesselJ[l,kr]-SphericalBesselJ[l,ka]SphericalBesselY[l,kr])/H1;
];
ans
]]

Table[Plot[Re[A[l,1.0,kr]],{kr,0,5},PlotRange->{-0.14,0.14}],{l,0,7}]

Do[anipartial1d[l]=Table[Plot[Re[A[l,1.0,kr]Exp[-I 2 π t]],{kr,0,5},PlotRange->{-0.14,0.14}],{t,0,0.99,0.01}],{l,0,7}];

Save["anipartial1d.mathematica",anipartial1d]

ψplusl=Compile[{{l,_Integer},{kred,_Real},{x,_Real},{y,_Real}},Module[{r,cosθ,kr},
r=Sqrt[(x^2)+(y^2)];
If[r==0.0,cosθ=1.0,cosθ=x/r];
(2l+1)I^l A[l,kred,r] LegendreP[l,cosθ]
]]

partialwaves=Table[DensityPlot[Re[ψplusl[l,1.0,x,y]],{x,-5,5},{y,-5,5},PlotPoints->100,PlotRange->All,Epilog->{Disk[]}],{l,0,7}]

SetDirectory[NotebookDirectory[]];Do[Export["p7230_hard_sphere_scattering_l="<>ToString[l]<>".jpg",partialwaves[[l+1]]],{l,0,7}]

With[{aa=0.14,l=0},ToExpression["anipartial"<>ToString[l],InputForm,Function[name,
name=Table[DensityPlot[(1+Re[ψplusl[l,1.0,x,y]Exp[-I 2 π t]]/aa)/2,{x,-5,5},{y,-5,5},PlotPoints->50,PlotRange->{-1,1},ColorFunctionScaling->False, Epilog->{Disk[]}],{t,0,0.9,0.1}],HoldAll]]];

Export["anipartial0.mov",anipartial[0]]

ψplus=Compile[{{kred,_Real},{x,_Real},{y,_Real},{lmax,_Integer}},Module[{r,cosθ,kr},
r=Sqrt[(x^2)+(y^2)];
If[r==0.0,cosθ=1.0,cosθ=x/r];
Sum[(2l+1)I^l A[l,kred,r] LegendreP[l,cosθ],{l,0,lmax}]
]]

DensityPlot[Re[ψplus[1.0,x,y,40]],{x,-5,5},{y,-5,5},PlotPoints->100,PlotRange->All,Epilog->{Disk[]}]

σscattering=Compile[{{kred,_Real},{lmax,_Integer}},Module[{ka},
ka=2 π kred;
((4π)/(ka^2))Sum[(2l+1)(((SphericalBesselJ[l,ka])^2)/(((SphericalBesselJ[l,ka])^2)+((SphericalBesselY[l,ka])^2))),{l,0,lmax}]
]]

Plot[σscattering[ka,20]/π,{ka,0.001,3.0},PlotRange->{0,4.2},AxesLabel->{"k a","σ/(π (a^2))"}]

Animations

Create and save the frames

This part is very slow to compute. It could take hours.

anifree=Table[DensityPlot[Cos[2π (1.0 x-t)],{x,-5,5},{y,-5,5},PlotPoints->50,PlotRange->All],{t,0,0.99,0.01}];

SetDirectory[NotebookDirectory[]];

Export["anifree.swf",anifree]

Save["anifree.mathematica",anifree]

anihardsphere=Table[DensityPlot[Re[ψplus[1.0,x,y,40]Exp[-2π I t]],{x,-5,5},{y,-5,5},PlotPoints->50,PlotRange->All,Epilog->{Disk[]}],{t,0,0.99,0.01}];

SetDirectory[NotebookDirectory[]];

Export["anihardsphere.swf",anihardsphere]

Save["anihardsphere.mathematica",anihardsphere]

Get the frames

SetDirectory[NotebookDirectory[]];

Get["anihardsphere.mathematica"];

open created mov file with Quicktime and Export as 480p

Export["anihardsphere.mov",anihardsphere]

Get["anifree.mathematica"];

Export["anifree.mov",anifree]