Visualizing Plane Wave Partial Wave Expansions

Mathematica code as it appeard in Mathematica in Education,
Volume 2 Number 1, Fall 1992

William J. Thompson
University of North Carolina at Chapel Hill
Department of Physics and Astronomy
Chapel Hill, NC  27599-3255

(Mathematica version 2.1)

Print["Plane waves in 3 D; wave along z\n"]
(* Cartesian -> plane polar functions *)
rho := Sqrt[z^2 + x^2]
theta := ArcTan[z,x]

(*	Spherical Bessel function *)
jl[l_,rho_] := Sqrt[Pi/(2 rho)]*
	 BesselJ[l+(1/2),rho]

L = Input["Maximum l value ( >=0 )"];

rL = Floor[L/2];
iL = Floor[(L-1)/2];
rhoMx = Input["Max rho (=kr) ( >0 )"];

Print["l <= ",L,", & rho <= ",rhoMx]
NumZX = 20; (* points in z & x *)
stepZX = rhoMx/(NumZX-1);

(* Real part of complex exponential *)
ReExp3D = Table[N[Cos[rho Cos[theta]]],
 {z,0.00001,rhoMx,stepZX},
 {x,0.00001,rhoMx,stepZX}];

(* Real part of sum to L *)
RePW3D = 
Table[NSum[
		(-1)^lp (4*lp + 1)*
		jl[2*lp,rho]*
		LegendreP[2*lp,Cos[theta]],
	  {lp,0,rL}],
 {z,0.00001,rhoMx,stepZX},
 {x,0.00001,rhoMx,stepZX}];

$DefaultFont = {"Times-Roman",8};

Figure 1

ReExpGraph = ListPlot3D[ReExp3D,
	DisplayFunction->Identity,
  	Ticks->None,
  	AxesLabel->{"x","z","Re"},
	PlotLabel->"Re exp[i rho cos(theta)]"];

RePWGraph = ListPlot3D[RePW3D,
 	DisplayFunction -> Identity,
 	Ticks -> None,
 	AxesLabel -> {"x","z","Re"},
 	PlotLabel -> "Re Partial Wave Expansion"];

(* Show Re side by side *)
Show[GraphicsArray[
 {ReExpGraph,RePWGraph},GraphicsSpacing->{.25,0}],
DisplayFunction -> $DisplayFunction];

Figure 2

(* Imaginary part of exponential *)
ImExp3D = Table[N[Sin[rho Cos[theta]]],
{z,0.00001,rhoMx,stepZX},
{x,0.00001,rhoMx,stepZX}];

(* Imaginary part of sum to L *)
ImPW3D = 
 Table[NSum[(-1)^lp (4*lp+3)*jl[2*lp+1,rho]*LegendreP[2*lp+1,Cos[theta]],
       {lp,0,iL}],
 {z,0.00001,rhoMx,stepZX},
 {x,0.00001,rhoMx,stepZX}];

ImExpGraph = 
 ListPlot3D[ImExp3D,
    DisplayFunction->Identity,
    Ticks->None,
    AxesLabel->{"x","z","Im"},
    PlotLabel->"Im exp[i rho cos(theta)]"];

ImPWGraph = 
 ListPlot3D[ImPW3D,
    DisplayFunction->Identity,
    Ticks->None,
    AxesLabel->{"x","z","Im"},
    PlotLabel->"Im Partial Wave Expansion"];

(* Show Im side by side *)
Show[GraphicsArray[{ImExpGraph,ImPWGraph}],
  DisplayFunction->$DisplayFunction];