Implementation of the 2-D MOI algorithm

This is a heavily annotated version of the MOI model parameters and functions. If you are very 
familiar with Mathematica, then you may want to look at MOI.nb, which is more streamlined.

A more complete, updated version is available on the Wolfram website:
https://community.wolfram.com/groups/-/m/t/3256594

Definitions of user interface symbols and their default values

This Mathematica code is commented and not very compactly written to help make it more understandable.
The default values used for various parameters match those used to generate the results presented in the paper.

Rectangular cavity characteristics (2 - D)

Lx - Long dimension of cavity (will be in x direction) (meters)
Ly - Short dimension (y direction, meters)

Qfunc[f] - cavity Q as a function of frequency f (kilohertz).

Predefined, default values (what were used for paper’s model):

Lx=0.152247;
Ly = 0.114417;

(* this is the linear Q(f) model used in the paper *)
(* change to   Qfunc[f_]:= Q   to use the value you assign to symbol Q, independent of f *)
Qfunc[f_]:=Evaluate[Simplify[((220.-80.)/(4500.-1100.))(f -1100.)+80.]]

Source and Receiver position offsets - all offsets should be ≥ 0 - away
from their respective, diagonally opposite corners

Sx - Source offset away from its corner in x direction (meters)
Sy - in y direction (meters)

x - Detector offset away from its corner in x direction (meters)
y - in y direction (meters)
Note: [x, y] values will usually be supplied as function arguments

Default values (what were used for paper’s model):

Sx=0.;
Sy = 0.;

x = 0.;
y = 0.;

Cs - speed of Sound for freq ↔ wave number conversions (m/sec)

Cs = 344.794;

Nmax - number of image clusters away from source corner along long cavity dimension Lx
(will go ± 2×Nmax×Lx away from source corner in both x and y directions, forming a 
square-shaped area containing source + images)

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.

Clear[Rtable];
Rtable[x_,y_]:=
Block[{Rfunc, X=1.x, Y=1.y, Ny = Ceiling[(Lx/Ly)Nmax]},

(* 
   Rfunc[i,j] is a local function which returns the distances from {x,y} 
  to the source images clustered around an indexed lattice location. 
 By using Block[], it is given a local redefinition with each 
  call to Rtable.
*)

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

A few examples of Rtable using various parameter values

Call with default source and detector locations in opposite corners. The suffix //Timing 
will return the CPU execution time (in seconds) required along with the result of the expression.

(* Execute Rtable[0,0] and return the length of the resulting list *)
Length[distances=Rtable[0,0]]//Timing

(* Execute Rtable[0,0] and return the value of the 20,000th member (count,distance) *)
Rtable[0,0][[20000]]//Timing

Now call with the detector located in the center of the cavity. The list will be longer because
of the reduced symmetry of the images around the detector, taking longer to gather and sort.

(* Return the length of the resulting list *)
Length[Rtable[Lx/2,Ly/2]]//Timing

(* The 20,000th member (count,distance) *)
Rtable[Lx/2,Ly/2][[20000]]//Timing

With both the source and detector positioned away from their corners, the symmetry is broken, and
the list of distances is approximately 16 times that for both in their corners. So is the time it takes to
generate the results. Using Block[] gives temporary, local values to Sx and Sy.

(* Return the length of the resulting list *)
Block[{Sx=Lx/4,Sy=Ly/4},Length[Rtable[Lx/2,Ly/2]]]//Timing

(* The 20,000th member (count,distance) *)
Block[{Sx=Lx/4,Sy=Ly/4},Rtable[Lx/2,Ly/2][[20000]] ]//Timing

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

The wave number argument k must have units of cycles/meter, not radians/meter. This makes it 
quicker to get the phase of kd, which will be 2π times its fractional part, rather that use modulo 2π. 
Thus: k=f/Cs.

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 the Terms over the list generated by Rtable

Response takes a list returned by Rtable as its first argument. If a single frequency in kHz is its
second argument, Response returns the sum of Term[] calls over the elements in Rtable, generating
the complex-valued response phasor at the supplied 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 Response’s second argument, then it returns a list of
{frequency, complex response} pairs suitable for generating frequency vs. gain or phase plots, etc.

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]