SetSystemOptions["CatchMachineUnderflow" -> False];

prepend this to my previous code.

I translated MATLAB code to Mathematica code directly.

A bottleneck is in MagneticSlabTransmission.

If xLim=54, underflow occurs in MagneticSlabTransmission,

and Mathematica will use arbitrary precision arithmetic, not FPU.

trapz[x_, y_] := (x[[2 ;;]] – x[[;; -2]]).(y[[2 ;;]] + y[[;; -2]])/2;

fftshift[x_] := RotateRight[x, Quotient[Length[x], 2]];

PerformDft1[fh_, xLim_, smpRateX_] :=
Module[{step, Nw, kLim, kstep, xvec, kvec, data, dataTransform},
step = smpRateX^-1;
Nw = 2*xLim*smpRateX;
kLim = smpRateX;
kstep = (2*xLim)^-1;
xvec = Range[-xLim, xLim - step, step];
kvec = Range[-kLim/2, kLim/2 - kstep, kstep];
data = fftshift[fh[xvec]];
dataTransform = step*Fourier[data, FourierParameters -> {1, -1}];
{dataTransform, kvec}
];

TESlabTransmissionNice[kx_, mu11_, d_] :=
Module[{eps, m0, epar, mperp, mu1, kz0, kz1, kz2},
eps = $MachineEpsilon;
m0 = 1. + eps*I;
epar = 1. + eps*I;
mperp = 1.;
mu1 = mu11 + eps*I;
kz0 = Sqrt[m0 - kx^2];
kz1 = Sqrt[mu1*(epar - kx^2/mperp)];
kz2 = kz1*m0/(kz0*mu1);
(Cos[d*kz1] – (kz2 + kz2^(-1))*Sin[d*kz1]*I/2)^(-1)
];

MagneticSlabTransmission[kx_, mu1_, d_] :=
With[{h = 0.5},
Exp[-4*Pi*h*Sqrt[kx^2 - 1 + 0.*I]]*
TESlabTransmissionNice[kx, mu1, 2*Pi*d](*
unpacking occurs . bottleneck *)
];

CostFunctionIntegrand[fh_] :=
Module[{xLim, xSmpRate, dftStruct, absxform},
xLim = 100.;
xSmpRate = 100.;
dftStruct = PerformDft1[fh, xLim, xSmpRate];
absxform = Abs[fftshift[dftStruct[[1]]]];
{dftStruct[[2]], absxform/Max[absxform]}
];

Timing[
outfile = OpenWrite["costfn1.dat"];
Do[
Write[outfile, {d, mpp,
trapz @@
CostFunctionIntegrand[MagneticSlabTransmission[#, mpp, d] &]}],
{d, 0.05, 1, 0.05}, {mpp, 0.05, 1, 0.05}];
Close[outfile];
]

Timing[n = 100;
zlim = 2*Pi*{-1, 2};
zpts = 5*n;
z = N[Range[zlim[[1]], zlim[[2]], (zlim[[2]] – zlim[[1]])/zpts]];
xlim = 10.;
xpts = 3*n;
x = Range[-xlim, xlim, 2*xlim/xpts];
klim = 1.;
kpts = 1*n;
kx = Range[-klim, klim, 2*klim/kpts];
kxmask = Exp[-(Pi/2*kx/klim)^2];
kz0 = Sqrt[1 - kx^2];
fwdPropFieldK =
DiagonalMatrix[kxmask].Exp[I*Transpose[Outer[Times, z, kz0]]];
xxkxx = Exp[I*Transpose[Outer[Times, x, kx]]];
beamMtx = Re[Transpose[xxkxx].fwdPropFieldK];]

This runs in around 0.17 seconds on my machine. Rendering the plot (code below) took around 1.1 seconds.

Timing[ListDensityPlot[beamMtx, ColorFunction -> "Rainbow",MaxPlotPoints -> 200]]

Daniel Lichtblau
Wolfram Research

Matlab? 48 seconds. still.
Now, I have to disqualify the NIntegrate results, because for several values of parameters, they seem to be numerically different from the Integrate[] or from matlab’s trapz by about 5%. (Matlab agrees with the Integrate[] approach to within 0.0003%)

(1) Use of exact inputs that will slow some functions (trigs, Sqrt[], and so on).

(2) Repeated omputation of a function that is mapped over a list. Much faster to have the function take a list argument and compute a list result. This is in part because we can now work with the vector as a packed array, and do not need to break it up element-by-element.

Here is the faster code. It runs in about 80 seconds on my machine.

FftShift1D[x_] := Module[{n = Ceiling[Length[x]/2]},
RotateRight[x, n]]

Out[8]= {50.1584, Null}

Locating the bottlenecks was in part knowing to look for (1) above, coupled with use of Print[] statements to isolate (2).

Daniel Lichtblau
Wolfram Research

It’s still nowhere near 47 seconds, but I thought you might be interested anyway!

Jony