Re: Do-loop conversion

• To: mathgroup at smc.vnet.net
• Subject: [mg50265] Re: Do-loop conversion
• From: Paul Abbott <paul at physics.uwa.edu.au>
• Date: Mon, 23 Aug 2004 06:34:54 -0400 (EDT)
• Organization: The University of Western Australia
• References: <cg6s9i\$o9l\$1@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```In article <cg6s9i\$o9l\$1 at smc.vnet.net>,
"Rich Matyi" <rjmatyi at comcast.net> 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

Factor[X[3]]

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

Factor[X[1]]

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

Cheers,
Paul

--
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 physics.uwa.edu.au
AUSTRALIA                            http://physics.uwa.edu.au/~paul

```

• 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