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