Implementation of the 2-D MOI algorithm

Definitions of user interface symbols and their default values

(* rectangular cavity dimensions (meters) *)
Lx=0.152247; (*  long dimension *)
Ly = 0.114417; (* short dimension *)

(* this is the linear Q(f) model used in the paper *)
(* change the definition if desired *)
Qfunc[f_]:=Evaluate[Simplify[((220.-80.)/(4500.-1100.))(f -1100.)+80.]]

(* source position offsets from its corner (meters); all offsets should be ≥ 0 *)
Sx=0.;
Sy = 0.;

(* response position offsets from diagonally opposite corner (meters) *)
(* the (x,y) response position must be provided as arguments 
	to the response function, so these values need not apply. *)
x = 0.;
y = 0.;

(* speed of sound (m/s); this is the value used in paper's model *)
Cs = 344.794;

(* number of image clusters away from source corner along Lx direction *)
(* will go ±2\[Cross]Nmax\[Cross]Lx away from source corner in both Lx and Ly directions*)
Nmax = 300;

Functions for the MOI Algorithm

Rtable: Create the list of distances from image sources to detector(x,y)

Rtable[x,y] uses global values set for Lx, Ly, Sx, Sy, Nmax. 
Response offsets supplied as arguments. 
Returns a sorted list (largest first) of {count, distance} pairs to source and images.

(* Rtable[x,y]: create the list of distances from image sources to response (x,y) *)
(* uses global values set for Lx, Ly, Sx, Sy, Nmax *)
(* returns list of {count, distance} pairs *) 
Clear[Rtable];
Rtable[x_,y_]:=
Block[{Rfunc, X=1.x, Y=1.y, Ny = Ceiling[(Lx/Ly)Nmax]},

(* Rfunc[i,j]: return the distances from {x,y} to the source images 
clustered around an indexed lattice location *) 
If[Sx==0.&&Sy==0., (* look at source position in cavity *)
(* Sx = Sy = 0, so source is in a corner, and only one image 
at each lattice point *)
Rfunc = 
Compile[{{i, _Integer},{j,_Integer}}, (* compile for speed*)
Evaluate[√(((2i -1)Lx+X)^2+((2j -1)Ly+Y)^2)],
RuntimeOptions->"Speed"
],
(* else we calculate 4 images around each lattice point, ±Sx, ±Sy *)
Rfunc = 
Compile[{{i, _Integer},{j,_Integer}}, 
Evaluate[
{
√(((2i -1)Lx-Sx+X)^2+((2j -1)Ly-Sy+Y)^2),
√(((2i -1)Lx-Sx+X)^2+((2j -1)Ly+Sy+Y)^2),
√(((2i -1)Lx+Sx+X)^2+((2j -1)Ly-Sy+Y)^2),
√(((2i -1)Lx+Sx+X)^2+((2j -1)Ly+Sy+Y)^2)
}
],
RuntimeOptions->"Speed"
]
]; (* end If *)

(* Here's what the next, quite powerful, expression does:  
1. Table[]: iterates over the lattice, calculating distances using Rfunc[i,j];
  2. //Flatten: gets rid of all inner braces, leaving a simple list of distances;
  3. //Gather[]&: organizes the list into sublists of equal-valued distances;
  4. //Sort[]&: sorts the sublists from longest to shortest distances to improve precision;
  5. //Map[]: replaces each sublist with a {count, distance} pair.
*)
Table[Rfunc[i,j], {i, -Nmax, Nmax}, {j, -Ny, Ny}]
//Flatten 
 //Gather[#,Equal]&
   //Sort[#, First[#1]>First[#2]&] &
     //Map[{1.Length[#],First[#]}&] (* note that count will be type Real *)
]

(* Execute Rtable[0,0] and return the length of the resulting list *)
(* Timing returns {CPU time (sec), result} *) 
Length[distances=Rtable[0,0]]//Timing  (* output saved in distances *)

Term: calculate one term of the sum of complex amplitudes from the sources

(* process a single element of Rtable output to get its complex phasor contribution *)
Term = Compile[{
	{k,_Real}, (* wave number in cycles/meter (f/Cs) *)
	{Q,_Real}, (* Q at this k (or freq) *)
	{m,_Real}, (* count from Rtable entry *)
	{d,_Real}     (* distance from Rtable entry *)
	},
	Evaluate[(m/(Sqrt[d])) Exp[I 2 Pi * FractionalPart[k * d]]Exp[-Pi*k*d/Q]]]

Response: Total  Term[ ] calls over the list generated by Rtable

(* If a single frequency in kHz is its second argument returns the response phasor at that frequency *) 
Response[
	RtableList:{{_Real,_Real}..} , (* list of results from call to Rtable *)
	f_ (* frequency (Hz) to be evaluated *)
	]:= 
	Block[
	{k=f/Cs,q=Qfunc[1.f]},
	Total[Term[k,q,First[#],Last[#]]&/@RtableList]
	]

Response[distances,1500]

(* If a list of frequencies is input as the second argument, it returns a list of
{frequency,complex response} pairs *) 
Response[
	RtableList:{{_Real,_Real}..} , (* list of results from call to Rtable *)
	fList:{__} (* list of frequencies (Hz) to be evaluated *)
	]:= Map[{#,Response[RtableList,#]}& ,fList]

Length[model=Response[distances,Evaluate[Range[900.,5000.,10.]]]]//Timing

ListLogPlot[Abs[model],Joined->True,Frame->True]

ListPlot[Map[{First[#],Arg[Last[#]]/Degree}&,model],Joined->True,Frame->True]