Do-loop conversion
- To: mathgroup at smc.vnet.net
- Subject: [mg50211] Do-loop conversion
- From: "Rich Matyi" <rjmatyi at comcast.net>
- Date: Sat, 21 Aug 2004 03:04:12 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
Mathematica colleagues,
I have a problem that involves converting a Do-loop structure into
faster-executing version using functional programming. This is an optics
problem -- specifically, calculating the angle dependence of the X-ray
reflectivity of an n-layered medium on an infinitely thick substrate. The
calculation relies on a well-known recursion relation:
X(j) = [ r(j) + X(j+1) p^2(j+1) ]/[1+ r(j) X(j+1) p^2(j+1)]
where X(j) is the ratio of the reflected and transmitted X-ray amplitudes
at the bottom of layer j, p(j+1) is a phase factor (given below), and r(j)
is the Fresnel coefficient for reflection at the interface between layers j
and j+1:
r(j) = [k(z,j) - k(z,j+1)]/[k(z,j) + k(z,j+1)]
where the wavevector k(z,j) = (2*Pi/wavelength)Sqrt[n^2(j) - Cos(phi)^2],
n(j) is the complex index of refraction for X-rays for the layer n(j) = 1 -
delta(j) + I*beta(j), and phi is incident angle of the X-rays. The phase
factor mentioned above is given by
p(j) = r(j) exp[-2 (k(j) k(j+1)) s(j+1)] with s being the roughness at
interface j+1.
The recursion layer works because with an infinitely thick substrate, the
reflected amplitude coming up from the substrate = 0, so at a given angle
of incidence, you work from the bottom up starting with X(bottom) = 0 and
add the amplitudes at each interface until you get to the top surface. The
Mathematica code that works just fine as a Do-loop is:
Do[
wavevectmap := wavetabc[revdelta, revbeta, phi];
leftvectmap := leftvectmapc[wavevectmap];
thtest := thtestcomp[revrough, leftvectmap, wavevectmap];
maproughrev := roughrevcomp[thtest];
coeffmap := coeffcomp[wavevectmap, leftvectmap, maproughrev];
phasemap := phasecomp[wavevectmap, revthick];
coeffs := fcoeffcomp[coeffmap, phasemap];
recurr[prior_, {fc_, pf_}] := (fc + (prior*pf^2))/(1 +
(fc*prior*pf^2));
intensity := intensitycomp[Fold[recurr, 0, coeffs ]];
newdata = newdatacomp[phi, intensity, params];
newc = AppendTo[newc, newdata],
{phi, phistart, phiend, phiinc}]
The various functions above -- wavetabc, leftvectmapc, thtestcomp,
roughrevcomp, coeffcomp, phasecomp, fcoeffcomp, intensitycomp, and
newdatacomp -- are complied functions that take structural parameters for
each layer (index of refraction parameters delta and beta, layer thickness,
interfacial roughness) by reading them out of a list {params} that is input
at the start of the program. I've left out these functions for brevity; if
needed, I'll post them as well.
As I said, the above Do-loop works just fine by calculating at each angle
phi the final intensity at that angle, normalizing it with a background and
maximum peak factor (also in the {params} list and applied by the
newdatacomp line) and appending it to the intensity datafile
newc. Executing this loop for a five-layer system and a couple of hundred
angular points takes around 0.2 seconds -- not bad, except that I need to
use this in a genetic algorithm minimization routine which needs many
hundreds of these curves to be calculated with GA- dictated changes to the
initial parameter matrix. So: if anyone can give me some guidance on how to
eliminate the Do-loop over the angular range so I can speed things up, I
would be very grateful.
Thanks in advance,
Rich Matyi
rjmatyi at comcast.net
- Follow-Ups:
- Re: Do-loop conversion
- From: DrBob <drbob@bigfoot.com>
- Re: Do-loop conversion