MathGroup Archive 2004

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

Search the Archive

Re: Do-loop conversion

  • To: mathgroup at
  • Subject: [mg50265] Re: Do-loop conversion
  • From: Paul Abbott <paul at>
  • Date: Mon, 23 Aug 2004 06:34:54 -0400 (EDT)
  • Organization: The University of Western Australia
  • References: <cg6s9i$o9l$>
  • Sender: owner-wri-mathgroup at

In article <cg6s9i$o9l$1 at>,
 "Rich Matyi" <rjmatyi at> wrote:

> 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],

You have k(z,j) on the left-hand side but wavelength and phi on the 
right-hand side. To turn this into a Mathematica function the 
relationship between (implicit) variables needs to be made clear.

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

Shouldn't these be k(z,j) and k(z,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. 

Entering the recurrence formula for X(j),

 x[j_] := (x[j + 1] p[j + 1]^2 + r[j])/(r[j] x[j + 1] p[j + 1]^2 + 1)

and the bottom condition,

 X[5] = 0; 

you can compute the results for intermediate layers, say


or for the top layer (depending on your numbering system),


in terms of r[j] and p[j]. These results are valid for _any_ angle phi 
and are not all that complicated. You could use dynamic programming
(x[j_] := x[j] = ...) to make this more efficient, if required. 

> 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

z, I assume?

> interfacial roughness) by reading them out of a list {params} that is input
> at the start of the program. 

Both r[j] and p[j] depend (implicitly) on z, lambda, and phi. Computing 
these functions should be very fast for specific parameter values.

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

If I understand what you're trying to do, I think that you can leave phi 
as a symbolic parameter. It appears that all other parameters will be 
numerical so that the final result, although complicated, should be 


Paul Abbott                                   Phone: +61 8 9380 2734
School of Physics, M013                         Fax: +61 8 9380 1014
The University of Western Australia      (CRICOS Provider No 00126G)         
35 Stirling Highway
Crawley WA 6009                      mailto:paul at 

  • Prev by Date: Re: Re: Beware of adding 0.0
  • Next by Date: Re: Do-loop conversion
  • Previous by thread: Re: Do-loop conversion
  • Next by thread: Re: Do-loop conversion