MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

Do-loop conversion

  • To: mathgroup at
  • Subject: [mg50211] Do-loop conversion
  • From: "Rich Matyi" <rjmatyi at>
  • Date: Sat, 21 Aug 2004 03:04:12 -0400 (EDT)
  • Sender: owner-wri-mathgroup at

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:

     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 +
     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

  • Prev by Date: Re: Re: Re: FindMinimum and the minimum-radius circle
  • Next by Date: Re: Re: Plotting a contour plot with cylindrical co-ordinates
  • Previous by thread: PDE with boundary condition ODE
  • Next by thread: Re: Do-loop conversion