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